Reprojecting Images to Different Observers

This example demonstrates how you can reproject images to the view from different observers, we use both AIA and STEREO A data to demonstrate this.

You will need reproject v0.6 or higher installed.

# sphinx_gallery_thumbnail_number = 2

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 STEREO A, when their separation was around 90 degrees.

stereo = (a.vso.Source('STEREO_A') &
          a.Instrument('EUVI') &
          a.Time('2010-08-19', '2010-08-19T00:10:00'))
aia = (a.Instrument('AIA') &
       a.vso.Sample(24 * u.hour) &
       a.Time('2010-08-19', '2010-08-19T00:10:00'))
wave = a.Wavelength(17 * u.nm, 18 * u.nm)

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

Out:

Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]


20100819_000715_n4eua.fts:   0%|          | 0.00/8.41M [00:00<?, ?B/s]


20100819_000715_n4eua.fts:   7%|6         | 558k/8.41M [00:00<00:01, 5.56MB/s]


20100819_000715_n4eua.fts:  40%|####      | 3.38M/8.41M [00:00<00:00, 7.33MB/s]


20100819_000715_n4eua.fts:  81%|########1 | 6.85M/8.41M [00:00<00:00, 9.60MB/s]


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

Create a map for each image.

map_aia, map_stereo = sunpy.map.Map(sorted(files))

# We downsample these maps to reduce memory consumption, you can comment this
# out.
map_aia = map_aia.resample((512, 512) * u.pix)
map_stereo = map_stereo.resample((512, 512) * 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_stereo)
map_stereo.plot(axes=ax2)
../../../_images/sphx_glr_reprojection_different_observers_001.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/io/fits/card.py:1003: VerifyWarning: Card is too long, comment will be truncated.
  VerifyWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/io/fits/card.py:1003: VerifyWarning: Card is too long, comment will be truncated.
  VerifyWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

<matplotlib.image.AxesImage object at 0x7f1403c5c990>

We now need to construct an output WCS. Because we want to downsample the image (for performance) we build a custom header using sunpy.map.make_fitswcs_header but we use a lot of the map_aia properties to do it. We use the reference coordinate from the AIA image, and make the scale 4x larger to compensate for the fact we are making the resolution 4x lower.

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

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

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

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

We can now reproject the STEREO 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.

output, footprint = reproject_interp(map_stereo, out_wcs, out_shape)

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

outmap = sunpy.map.Map(output, out_header)
outmap.plot_settings = map_stereo.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)
../../../_images/sphx_glr_reprojection_different_observers_002.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

<matplotlib.image.AxesImage object at 0x7f14086588d0>

AIA as Seen from Mars

We can also change the observer of the AIA image to any observer coordinate. SunPy provides a function which can get the observer coordinate for any known body. In this example we are going to use Mars.

mars = get_body_heliographic_stonyhurst('mars', map_aia.date)

We now generate a target WCS, to do this we need to generate a reference coordinate, which is similar to the aia frame, but scaled down by 4x and with the observer at Mars. To do this we generate a new reference coordinate.

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

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

then a header

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

Once again we need to generate a WCS 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)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

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)
outmap.draw_grid(color='w')

plt.show()
../../../_images/sphx_glr_reprojection_different_observers_003.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -6.17129447065
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

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

Gallery generated by Sphinx-Gallery