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
|
"""
=============================
Differentially rotating a map
=============================
How to apply differential rotation to a Map.
.. note::
This example requires `reproject` 0.6 or later to be installed.
The example uses the :func:`~sunpy.coordinates.propagate_with_solar_surface`
context manager to apply differential rotation during coordinate
transformations.
"""
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import sunpy.map
from sunpy.coordinates import Helioprojective, propagate_with_solar_surface
from sunpy.data.sample import AIA_171_IMAGE
##############################################################################
# First, load an AIA observation.
aiamap = sunpy.map.Map(AIA_171_IMAGE)
in_time = aiamap.date
##############################################################################
# Let's define the output frame to be five days in the future for an observer
# at Earth (i.e., about five degrees offset in heliographic longitude compared
# to the location of AIA in the original observation).
out_time = in_time + 5*u.day
out_frame = Helioprojective(observer='earth', obstime=out_time,
rsun=aiamap.coordinate_frame.rsun)
##############################################################################
# Construct a WCS object for the output map. If one has an actual ``Map``
# object at the desired output time (e.g., the actual AIA observation at the
# output time), one can use the WCS object from that ``Map`` object (e.g.,
# ``mymap.wcs``) instead of constructing a custom WCS.
out_center = SkyCoord(0*u.arcsec, 0*u.arcsec, frame=out_frame)
header = sunpy.map.make_fitswcs_header(aiamap.data.shape,
out_center,
scale=u.Quantity(aiamap.scale))
out_wcs = WCS(header)
##############################################################################
# Reproject the map from the input frame to the output frame. We use the
# :func:`~sunpy.coordinates.propagate_with_solar_surface` context manager so
# that coordinates are treated as points that evolve in time with the
# rotation of the solar surface rather than as inertial points in space.
with propagate_with_solar_surface():
out_warp = aiamap.reproject_to(out_wcs)
##############################################################################
# Let's plot the differentially rotated Map next to the original Map.
fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(121, projection=aiamap)
aiamap.plot(axes=ax1, vmin=0, vmax=20000, title='Original map')
plt.colorbar()
ax2 = fig.add_subplot(122, projection=out_warp)
out_warp.plot(axes=ax2, vmin=0, vmax=20000,
title=f"Reprojected to an Earth observer {(out_time - in_time).to('day')} later")
plt.colorbar()
plt.show()
|
