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 v0.6 or higher 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 get_body_heliographic_stonyhurst
from sunpy.net import Fido
from sunpy.net import attrs as a

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)

Let’s download an EUV image from both AIA and EUVI A, when the two spacecraft were separated by approximately 120 degrees.

euvi = (a.vso.Source('STEREO_A') &
        a.Instrument("EUVI") &
        a.Time('2011-11-01', '2011-11-01T00:10:00'))

aia = (a.Instrument.aia &
       a.Sample(24 * u.hour) &
       a.Time('2011-11-01', '2011-11-02'))

wave = a.Wavelength(19.5 * u.nm, 19.5 * u.nm)

res = Fido.search(wave, aia | euvi)
files = Fido.fetch(res)

Out:

Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/connector.py:964: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  hosts = await asyncio.shield(self._resolve_host(
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/connector.py:964: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  hosts = await asyncio.shield(self._resolve_host(
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/locks.py:21: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._event = asyncio.Event(loop=loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/locks.py:21: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._event = asyncio.Event(loop=loop)

Files Downloaded:  50%|#####     | 1/2 [00:00<00:00,  7.38file/s]
Files Downloaded: 100%|##########| 2/2 [00:00<00:00,  3.15file/s]
Files Downloaded: 100%|##########| 2/2 [00:00<00:00,  2.28file/s]

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(files)
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)

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
ax2 = fig.add_subplot(1, 2, 2, projection=map_euvi)
map_euvi.plot(axes=ax2)
AIA $193 \; \mathrm{\mathring{A}}$ 2011-11-01 00:00:07, EUVI-A $195 \; \mathrm{\mathring{A}}$ 2011-11-01 00:05:31

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

<matplotlib.image.AxesImage object at 0x7fac7995c9d0>

We now need to construct an output WCS. We build a custom header using sunpy.map.make_fitswcs_header() but we use a lot of the map_aia properties to do it.

out_header = sunpy.map.make_fitswcs_header(
    out_shape,
    map_aia.reference_coordinate,
    scale=u.Quantity(map_aia.scale),
    instrument="EUVI",
    observatory="AIA Observer",
    wavelength=map_euvi.wavelength
)

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. '

Next we construct an WCS object from the header. Currently WCS does not understand the observer position, so we manually set that.

out_wcs = WCS(out_header)
out_wcs.heliographic_observer = map_aia.reference_coordinate.observer

We can now reproject the EUVI map to this output WCS. Here we are using the fastest but least accurate method of reprojection, reproject.reproject_interp(). A more accurate but slower method is reproject.reproject_adaptive().

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

We can now plot the STEREO/EUVI image as seen from the position of SDO, next to the AIA image.

outmap = sunpy.map.Map(output, out_header)
outmap.plot_settings = map_euvi.plot_settings

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
ax2 = fig.add_subplot(1, 2, 2, projection=outmap)
outmap.plot(axes=ax2, title='EUVI image as seen from SDO')
AIA $193 \; \mathrm{\mathring{A}}$ 2011-11-01 00:00:07, EUVI image as seen from SDO

Out:

<matplotlib.image.AxesImage object at 0x7fac7c26ba60>

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)

To generate a target WCS, we first need an appropriate reference coordinate, which is similar to the one for AIA, except now with the observer at Mars.

mars_ref_coord = SkyCoord(map_aia.reference_coordinate.Tx,
                          map_aia.reference_coordinate.Ty,
                          obstime=map_aia.reference_coordinate.obstime,
                          observer=mars,
                          frame="helioprojective")

We then create the WCS header.

mars_header = sunpy.map.make_fitswcs_header(
    out_shape,
    mars_ref_coord,
    scale=u.Quantity(map_aia.scale),
    rotation_matrix=map_aia.rotation_matrix,
    instrument="AIA",
    wavelength=map_aia.wavelength
)

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. '

Once again we need to generate a WCS object and then manually set the observer location.

mars_wcs = WCS(mars_header)
mars_wcs.heliographic_observer = mars

output, footprint = reproject_interp(map_aia, mars_wcs, out_shape)

We generate the output map and plot it next to the original image.

outmap = sunpy.map.Map((output, mars_header))
outmap.plot_settings = map_aia.plot_settings

fig = plt.figure()

ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
outmap.draw_grid(color='w')

ax2 = fig.add_subplot(1, 2, 2, projection=outmap)
outmap.plot(axes=ax2, title='AIA observation as seen from Mars')
outmap.draw_grid(color='w')

plt.show()
AIA $193 \; \mathrm{\mathring{A}}$ 2011-11-01 00:00:07, AIA observation as seen from Mars

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

Gallery generated by Sphinx-Gallery