AIA to STEREO coordinate conversion

How to convert a point of a source on an AIA image to a position on a STEREO image.

import matplotlib.pyplot as plt

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

import sunpy.map
import sunpy.coordinates
import sunpy.coordinates.wcs_utils
from sunpy.net import Fido, attrs as a

The first step is to download some data, we are going to get an image from early 2011 when the STEREO spacecraft were roughly 90 deg seperated from the Earth.

stereo = (a.vso.Source('STEREO_B') &
          a.Instrument('EUVI') &
          a.Time('2011-01-01', '2011-01-01T00:10:00'))

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

wave = a.Wavelength(30 * u.nm, 31 * u.nm)


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

Download the files:

files = Fido.fetch(res)
print(files)

Out:

Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]
Files Downloaded:  50%|#####     | 1/2 [00:00<00:00,  8.24file/s]
Files Downloaded: 100%|##########| 2/2 [00:03<00:00,  1.01s/file]
Files Downloaded: 100%|##########| 2/2 [00:03<00:00,  1.60s/file]
['/home/docs/sunpy/data/aia_lev1_304a_2011_01_01t00_00_08_12z_image_lev1.fits', '/home/docs/sunpy/data/20110101_000615_n4eub.fts']

Create a dictionary with the two maps, cropped down to full disk.

maps = {m.detector: m.submap(SkyCoord([-1100, 1100]*u.arcsec,
                                      [-1100, 1100]*u.arcsec,
                                      frame=m.coordinate_frame))
        for m in sunpy.map.Map(files)}

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 = -4.13136390199
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 = -4.13136390199
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 = -4.13136390199
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 = -4.13136390199
keyword looks very much like CROTAn but isn't.
  colsel=colsel)

Plot both maps

fig = plt.figure(figsize=(10, 4))
for i, m in enumerate(maps.values()):
    ax = fig.add_subplot(1, 2, i+1, projection=m)
    m.plot(axes=ax)
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_001.png

Out:

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

We are now going to pick out a region around the south west corner:

aia_width = 200 * u.arcsec
aia_height = 250 * u.arcsec
aia_bottom_left = SkyCoord([[-800, -300]] * u.arcsec,
                           frame=maps['AIA'].coordinate_frame)
aia_top_right = SkyCoord(aia_bottom_left.Tx + aia_width,
                         aia_bottom_left.Ty + aia_height,
                         frame=maps['AIA'].coordinate_frame)

Plot a rectangle around the region we want to crop

m = maps['AIA']
fig = plt.figure()
ax = fig.add_subplot(111, projection=m)
m.plot(axes=ax)
m.draw_rectangle(aia_bottom_left, aia_width, aia_height)
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_002.png

Out:

[<matplotlib.patches.Rectangle object at 0x7f142a02ea10>]

Create a submap of this area

subaia = maps['AIA'].submap(aia_bottom_left, aia_top_right)
fig = plt.figure()
subaia.plot()
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_003.png

Out:

<matplotlib.image.AxesImage object at 0x7f142a038190>

We now want to crop out this same area on the STEREO EUVI image. First, we create a SkyCoord object with the four corners of the box. When we create this object, we use Map.coordinate_frame so that the location parameters of SDO are correctly set.

corners = ([aia_bottom_left.Tx, aia_bottom_left.Ty],
           [aia_bottom_left.Tx + aia_width, aia_bottom_left.Ty],
           [aia_bottom_left.Tx, aia_bottom_left.Ty + aia_height],
           [aia_top_right.Tx, aia_top_right.Ty])

hpc_aia = SkyCoord(corners, frame=maps['AIA'].coordinate_frame)

print(hpc_aia)

Out:

<SkyCoord (Helioprojective: obstime=2011-01-01T00:00:08.120, rsun=696000000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2011-01-01T00:00:08.120): (lon, lat, radius) in (deg, deg, m)
    (-0.01598098, -2.97460829, 1.47100674e+11)>): (Tx, Ty) in arcsec
    [[(-800., -300.)],
     [(-600., -300.)],
     [(-800.,  -50.)],
     [(-600.,  -50.)]]>

We can now transform to from the AIA frame to the EUVI frame. This transformation first transforms to Heliographic Stonyhurst coordinates and then into the EUVI frame.

hpc_B = hpc_aia.transform_to(maps['EUVI'].coordinate_frame)
print(hpc_B)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/wcs/wcs.py:466: FITSFixedWarning: CROTA = -4.13136390199
keyword looks very much like CROTAn but isn't.
  colsel=colsel)
<SkyCoord (Helioprojective: obstime=2011-01-01T00:07:01.811, rsun=695700000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2011-01-01T00:07:01.811): (lon, lat, radius) in (deg, deg, m)
    (-89.58098377, 6.76626853, 1.58345292e+11)>): (Tx, Ty, distance) in (arcsec, arcsec, m)
    [[(422.14140721, -387.56099958, 1.57805355e+11)],
     [(642.54613611, -376.44172931, 1.57947290e+11)],
     [(513.08620817, -161.10016561, 1.57784422e+11)],
     [(709.38834261, -149.2088207 , 1.57926317e+11)]]>

Now we can plot this box on both the AIA and EUVI images:

fig = plt.figure(figsize=(10, 4))
for i, (m, coord) in enumerate(zip([maps['EUVI'], maps['AIA']],
                                   [hpc_B, hpc_aia])):
    ax = fig.add_subplot(1, 2, i+1, projection=m)
    m.plot(axes=ax)

    # coord[3] is the top-right corner coord[0] is the bottom-left corner.
    w = (coord[3].Tx - coord[0].Tx)
    h = (coord[3].Ty - coord[0].Ty)
    m.draw_rectangle(coord[0], w, h,
                     transform=ax.get_transform('world'))
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_004.png

Out:

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

We can now zoom in on the region in the EUVI image:

subeuvi = maps['EUVI'].submap(hpc_B[0], hpc_B[3])
fig = plt.figure()
plt.subplot(projection=subeuvi)
subeuvi.plot()
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_005.png

Out:

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

<matplotlib.image.AxesImage object at 0x7f1403ac0610>

Putting them together:

fig = plt.figure(figsize=(15, 5))
for i, m in enumerate((subeuvi, subaia)):
    ax = fig.add_subplot(1, 2, i+1, projection=m)
    m.plot(axes=ax)
plt.show()
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_006.png

Out:

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

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

Gallery generated by Sphinx-Gallery