Reprojecting a Map

How to apply differential rotation to a Map.

The example uses the RotatedSunFrame coordinate metaframe in sunpy.coordinates to apply differential rotation to a Map. See Differential rotation using coordinate frames for more details on using RotatedSunFrame.

Note

This example requires reproject 0.6 or later to be installed.

import matplotlib.pyplot as plt
from reproject import reproject_interp

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

import sunpy.map
from sunpy.coordinates import Helioprojective, RotatedSunFrame, transform_with_sun_center
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)

For the reprojection, the definition of the target frame can be counter-intuitive. We will be transforming from the original frame to the RotatedSunFrame version of the output frame with rotated_time set to the time of the original frame.

rot_frame = RotatedSunFrame(base=out_frame, rotated_time=in_time)
print(rot_frame)

Out:

<RotatedSunHelioprojective Frame (base=<Helioprojective Frame (obstime=2011-06-12T06:33:02.770, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>)>, duration=-5.0 d, rotation_model=howard)>

Construct a WCS object for the output map with the target RotatedSunHelioprojective frame specified instead of the regular Helioprojective frame. 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_shape = aiamap.data.shape
out_center = SkyCoord(0*u.arcsec, 0*u.arcsec, frame=out_frame)
header = sunpy.map.make_fitswcs_header(out_shape,
                                       out_center,
                                       scale=u.Quantity(aiamap.scale))
out_wcs = WCS(header)
out_wcs.coordinate_frame = rot_frame

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/units/quantity.py:1058: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError.
  warnings.warn('The truth value of a Quantity is ambiguous. '

Reproject the map from the input frame to the output frame. We use the transform_with_sun_center() context manager so that the coordinates stay defined relative to Sun center as the Sun moves slowly over time.

with transform_with_sun_center():
    arr, _ = reproject_interp(aiamap, out_wcs, out_shape)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/wcs_utils.py:224: RuntimeWarning: invalid value encountered in greater
  reset |= (np.abs(inputs_check[ipix] - inputs[ipix]) > 1)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/wcs_utils.py:224: RuntimeWarning: invalid value encountered in greater
  reset |= (np.abs(inputs_check[ipix] - inputs[ipix]) > 1)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/array_utils.py:30: RuntimeWarning: invalid value encountered in less
  reset |= (coords[i] < -0.5)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/array_utils.py:31: RuntimeWarning: invalid value encountered in greater
  reset |= (coords[i] > original_shape[i] - 0.5)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/array_utils.py:30: RuntimeWarning: invalid value encountered in less
  reset |= (coords[i] < -0.5)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/reproject/array_utils.py:31: RuntimeWarning: invalid value encountered in greater
  reset |= (coords[i] > original_shape[i] - 0.5)

Finally, we create the output map and preserve the original map’s plot settings.

Let’s plot the differentially rotated Map next to the original Map.

fig = plt.figure(figsize=(12, 4))

ax1 = fig.add_subplot(1, 2, 1, projection=aiamap)
aiamap.plot(vmin=0, vmax=20000, title='Original map')
plt.colorbar()

ax2 = fig.add_subplot(1, 2, 2, projection=out_warp)
out_warp.plot(vmin=0, vmax=20000,
              title=f"Reprojected to an Earth observer {(out_time - in_time).to('day')} later")
plt.colorbar()

plt.show()
Original map, Reprojected to an Earth observer 5.0 d later

Total running time of the script: ( 0 minutes 11.215 seconds)

Gallery generated by Sphinx-Gallery