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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
"""
==========================================
Reprojecting Images to Different Observers
==========================================
This example demonstrates how you can reproject images to the view from
different observers. We use data from these two instruments:
* AIA on SDO, which is in orbit around Earth
* EUVI on STEREO A, which is in orbit around the Sun away from the Earth
You will need `reproject <https://reproject.readthedocs.io/en/stable/>`__ v0.6 or higher installed.
"""
# sphinx_gallery_thumbnail_number = 2
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord
import sunpy.map
from sunpy.coordinates import get_body_heliographic_stonyhurst
from sunpy.data.sample import AIA_193_JUN2012, STEREO_A_195_JUN2012
######################################################################
# In this example we are going to make a lot of side by side figures, so
# let's change the default figure size.
plt.rcParams['figure.figsize'] = (16, 8)
######################################################################
# Create a map for each image, after making sure to sort by the
# appropriate name attribute (i.e., "AIA" and "EUVI") so that the
# order is reliable.
map_list = sunpy.map.Map([AIA_193_JUN2012, STEREO_A_195_JUN2012])
map_list.sort(key=lambda m: m.detector)
map_aia, map_euvi = map_list
# We downsample these maps to reduce memory consumption, but you can
# comment this out.
out_shape = (512, 512)
map_aia = map_aia.resample(out_shape * u.pix)
map_euvi = map_euvi.resample(out_shape * u.pix)
######################################################################
# Plot the two maps, with the solar limb as seen by each observatory
# overlaid on both plots.
fig = plt.figure()
ax1 = fig.add_subplot(121, projection=map_aia)
map_aia.plot(axes=ax1)
map_aia.draw_limb(axes=ax1, color='white')
map_euvi.draw_limb(axes=ax1, color='red')
ax2 = fig.add_subplot(122, projection=map_euvi)
map_euvi.plot(axes=ax2)
limb_aia = map_aia.draw_limb(axes=ax2, color='white')
limb_euvi = map_euvi.draw_limb(axes=ax2, color='red')
plt.legend([limb_aia[0], limb_euvi[0]],
['Limb as seen by AIA', 'Limb as seen by EUVI A'])
######################################################################
# Data providers can set the radius at which emission in the map is assumed
# to have come from. Most maps use a default value for photospheric
# radius (including EUVI maps), but some maps (including AIA maps) are set
# to a slightly different value. A mismatch in solar radius means a reprojection
# will not work correctly on pixels near the limb. This can be prevented by
# modifying the values for rsun on one map to match the other.
map_euvi.meta['rsun_ref'] = map_aia.meta['rsun_ref']
######################################################################
# We can reproject the EUVI map to the AIA observer wcs using
# :meth:`~sunpy.map.GenericMap.reproject_to`. This method defaults to using
# the fast :func:`reproject.reproject_interp` algorithm, but a different
# algorithm can be specified (e.g., :func:`reproject.reproject_adaptive`).
outmap = map_euvi.reproject_to(map_aia.wcs)
######################################################################
# We can now plot the STEREO/EUVI image as seen from the position of
# SDO, next to the AIA image.
fig = plt.figure()
ax1 = fig.add_subplot(121, projection=map_aia)
map_aia.plot(axes=ax1)
ax2 = fig.add_subplot(122, projection=outmap)
outmap.plot(axes=ax2, title='EUVI image as seen from SDO')
map_euvi.draw_limb(color='blue')
# Set the HPC grid color to black as the background is white
ax2.grid(color='k')
######################################################################
# AIA as Seen from Mars
# =====================
#
# The new observer coordinate doesn't have to be associated with an
# existing Map. sunpy provides a function which can get the location
# coordinate for any known body. In this example, we use Mars.
mars = get_body_heliographic_stonyhurst('mars', map_aia.date)
######################################################################
# Without a target Map wcs, we can generate our own for an arbitrary observer.
# First, we need an appropriate reference coordinate. This will be similar to
# the one contained in ``map_aia``, except with the observer placed at Mars.
mars_ref_coord = SkyCoord(0*u.arcsec, 0*u.arcsec,
obstime=map_aia.reference_coordinate.obstime,
observer=mars,
rsun=map_aia.reference_coordinate.rsun,
frame="helioprojective")
######################################################################
# We now need to construct our output WCS; we build a custom header using
# :func:`sunpy.map.header_helper.make_fitswcs_header` using the ``map_aia``
# properties and our new, mars-based reference coordinate.
mars_header = sunpy.map.make_fitswcs_header(
out_shape,
mars_ref_coord,
scale=u.Quantity(map_aia.scale),
instrument="AIA",
wavelength=map_aia.wavelength
)
######################################################################
# We generate the output map and plot it next to the original image.
outmap = map_aia.reproject_to(mars_header)
fig = plt.figure()
ax1 = fig.add_subplot(121, projection=map_aia)
map_aia.plot(axes=ax1)
map_aia.draw_grid(color='w')
ax2 = fig.add_subplot(122, projection=outmap)
outmap.plot(axes=ax2, title='AIA observation as seen from Mars')
map_aia.draw_grid(color='w')
map_aia.draw_limb(color='blue')
plt.show()
|
