AIA to STEREO Coordinate Conversion

In this example we demonstrate how you can identify a point or region on the surface of the Sun in an AIA image and then convert that point to a point in a STEREO image.


This example requires WCSAxes which is an optional SunPy dependency.

import matplotlib.pyplot as plt

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

import sunpy.coordinates
import sunpy.coordinates.wcs_utils
from import vso

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 = (vso.attrs.Source('STEREO_B') &
          vso.attrs.Instrument('EUVI') &
          vso.attrs.Time('2011-01-01', '2011-01-01T00:10:00'))

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

wave = vso.attrs.Wave(30 * u.nm, 31 * u.nm)

vc = vso.VSOClient()
res = vc.query(wave, aia | stereo)

The results from VSO query:



Start Time [1]       End Time [1]     Source  Instrument   Type
------------------- ------------------- -------- ---------- --------
2011-01-01 00:07:01 2011-01-01 00:07:05 STEREO_B     SECCHI FULLDISK
2011-01-01 00:00:08 2011-01-01 00:00:09      SDO        AIA FULLDISK

Download the files:

files = vc.get(res).wait()

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

maps = {m.detector: m.submap((-1100, 1100) * u.arcsec,
                             (-1100, 1100) * u.arcsec) for m in}

Plot both maps

fig = plt.figure(figsize=(15, 5))
for i, m in enumerate(maps.values()):
    ax = fig.add_subplot(1, 2, i+1, projection=m.wcs)

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 = (-800, -300) * u.arcsec

Plot a rectangle around the region we want to crop

m = maps['AIA']
fig = plt.figure()
ax = fig.add_subplot(111, projection=m.wcs)
m.draw_rectangle(aia_bottom_left, aia_width, aia_height)

Create a submap of this area

subaia = maps['AIA'].submap(u.Quantity((aia_bottom_left[0],
                                        aia_bottom_left[0] + aia_width)),
                                        aia_bottom_left[1] + aia_height)))

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.

hpc_aia = SkyCoord((aia_bottom_left,
                    aia_bottom_left + u.Quantity((aia_width, 0 * u.arcsec)),
                    aia_bottom_left + u.Quantity((0 * u.arcsec, aia_height)),
                    aia_bottom_left + u.Quantity((aia_width, aia_height))),



<SkyCoord (Helioprojective: D0=1.47100673949e+11 m, dateobs=2011-01-01 00:00:08.120000, L0=0.0 deg, B0=-2.974563 deg, rsun=695508.0 km): (Tx, Ty) in arcsec
    [(-800., -300.), (-600., -300.), (-800.,  -50.), (-600.,  -50.)]>

Now we convert these coordinates into Heliographic Stonyhurst coordinates, which are on the Sun, with the zero meridian facing the Earth.

hgs = hpc_aia.transform_to('heliographic_stonyhurst')


<SkyCoord (HeliographicStonyhurst: dateobs=2011-01-01 00:00:08.120000): (lon, lat, radius) in (deg, deg, km)
    [(-60.17483354, -19.3723875 ,  695507.99999902),
     (-40.76422128, -20.11635134,  695508.00000253),
     (-55.16235377,  -4.63378036,  695507.99999544),
     (-37.99160921,  -5.27440002,  695508.00000112)]>

Now we need to provide the position information from the STEREO Imager:

hgs.D0 = maps['EUVI'].dsun
hgs.L0 = maps['EUVI'].heliographic_longitude
hgs.B0 = maps['EUVI'].heliographic_latitude

We do this on the Heliographic frame because when in a Heliographic frame these parameters have no effect on the frame, but they are used when the frame is converted back to Helioprojective. And now we can convert back to Helioprojective, but this time from the view-point of STEREO B:

hpc_B = hgs.transform_to('helioprojective')


<SkyCoord (Helioprojective: D0=1.58345292023e+11 m, dateobs=2011-01-01 00:00:08.120000, L0=-89.5809837657 deg, B0=6.76626852815 deg, rsun=695508.0 km): (Tx, Ty, distance) in (arcsec, arcsec, km)
    [( 421.08835996, -387.47532689,   1.57805467e+08),
     ( 641.87198769, -376.36878696,   1.57947465e+08),
     ( 512.23675413, -161.01960073,   1.57784560e+08),
     ( 708.78219966, -149.13566688,   1.57926511e+08)]>

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

fig = plt.figure(figsize=(15, 5))
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.wcs)

    # 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(u.Quantity((coord[0].Tx, coord[0].Ty)), w, h,

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

subeuvi = maps['EUVI'].submap(u.Quantity((hpc_B[0].Tx, hpc_B[3].Tx)),
                              u.Quantity((hpc_B[0].Ty, hpc_B[3].Ty)))

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.wcs)

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

Generated by Sphinx-Gallery