1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
"""
===========================
Plotting a difference image
===========================
This example shows how to compute and plot a difference image.
"""
# sphinx_gallery_thumbnail_number = 3
import matplotlib.colors as colors
import matplotlib.pyplot as plt
from astropy.visualization import ImageNormalize, SqrtStretch
import sunpy.data.sample
import sunpy.map
###########################################################################
# When analyzing solar imaging data it is often useful to look at the
# difference from one time step to the next (running difference) or the
# difference from the start of the sequence (base difference). In this
# example, we'll use a sequence of AIA 193 cutout images taken during a
# flare.
#
# First, load a series of images into a `~sunpy.map.MapSequence`.
m_seq = sunpy.map.Map([
sunpy.data.sample.AIA_193_CUTOUT01_IMAGE,
sunpy.data.sample.AIA_193_CUTOUT02_IMAGE,
sunpy.data.sample.AIA_193_CUTOUT03_IMAGE,
sunpy.data.sample.AIA_193_CUTOUT04_IMAGE,
sunpy.data.sample.AIA_193_CUTOUT05_IMAGE,
], sequence=True)
###########################################################################
# Let's take a look at each image in the sequence. Note that these images
# are sampled at a 6.4 minute cadence, much lower than the actual 12 s
# AIA cadence. We adjust the plot setting to ensure the colorbar is
# the same at each time step.
plt.figure()
ani = m_seq.plot(norm=ImageNormalize(vmin=0, vmax=5e3, stretch=SqrtStretch()))
plt.show()
###########################################################################
# And now we can take the actual difference. We will compute the running
# difference and the base difference for each image in the sequence. For the
# case of the first entry in the running difference, we just subtract it
# from itself.
#
# But we have to decide what to do with the metadata. For example, what
# time does this difference image correspond to? The time of the first or
# second image? The mean time? You'll have to decide what makes most sense
# for your application. Here, by subtracting the data of the first (and
# previous) data array from each map in the sequence, the resulting
# difference maps have the same metadata as each corresponding map in the
# sequence.
#
# Note that, because arithmetic operations between `~sunpy.map.GenericMap`
# objects are not supported, we subtract just the array data (with units
# attached) of the second map from the first map. The ``quantity`` attribute
# returns the image data as an `~astropy.units.Quantity`, where the resulting
# units are those returned by the ``unit`` attribute of the map.
m_seq_base = sunpy.map.Map([m - m_seq[0].quantity for m in m_seq[1:]], sequence=True)
m_seq_running = sunpy.map.Map(
[m - prev_m.quantity for m, prev_m in zip(m_seq[1:], m_seq[:-1])],
sequence=True
)
###########################################################################
# Finally, let's plot the difference maps. We'll apply a colormap and
# re-normalize the intensity so that it shows up well.
# First, we show the base difference map.
plt.figure()
ani = m_seq_base.plot(title='Base Difference', norm=colors.Normalize(vmin=-200, vmax=200), cmap='Greys_r')
plt.colorbar(extend='both', label=m_seq_base[0].unit.to_string())
plt.show()
###########################################################################
# Then, we show the running difference map.
plt.figure()
ani = m_seq_running.plot(title='Running Difference', norm=colors.Normalize(vmin=-200, vmax=200), cmap='Greys_r')
plt.colorbar(extend='both', label=m_seq_running[0].unit.to_string())
plt.show()
|
