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.

from __future__ import print_function, division

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)

The results from VSO query:

print(res)

Out:

Results from 2 Providers:

1 Results from the VSOClient:
   Start Time [1]       End Time [1]    Source ...   Type   Wavelength [2]
                                               ...             Angstrom
       str19               str19         str3  ...   str8      float64
------------------- ------------------- ------ ... -------- --------------
2011-01-01 00:00:08 2011-01-01 00:00:09    SDO ... FULLDISK 304.0 .. 304.0

1 Results from the VSOClient:
   Start Time [1]       End Time [1]     Source  ...   Type   Wavelength [2]
                                                 ...             Angstrom
       str19               str19          str8   ...   str8      float64
------------------- ------------------- -------- ... -------- --------------
2011-01-01 00:07:01 2011-01-01 00:07:05 STEREO_B ... FULLDISK 304.0 .. 304.0

Download the files:

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

Out:

['/home/docs/sunpy/data/aia_lev1_304a_2011_01_01t00_00_08_12z_image_lev1.0.fits', '/home/docs/sunpy/data/20110101_000615_n4eub.0.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)}

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.wcs)
    m.plot(axes=ax)
../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_001.png

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

Create a submap of this area

subaia = maps['AIA'].submap(aia_bottom_left, aia_top_right)
subaia.peek()
../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_003.png

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-01 00:00:08.120000, rsun=696000000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2011-01-01 00:00:08.120000): (lon, lat, radius) in (deg, deg, m)
    ( 0., -2.974563,   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:

<SkyCoord (Helioprojective: obstime=2011-01-01 00:07:01.811000, rsun=695508000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2011-01-01 00:07:01.811000): (lon, lat, radius) in (deg, deg, m)
    (-89.58098377,  6.76626853,   1.58345292e+11)>): (Tx, Ty, distance) in (arcsec, arcsec, km)
    [[( 422.41887548, -387.542582  ,   1.57805472e+08)],
     [( 642.7557377 , -376.41356048,   1.57947468e+08)],
     [( 513.36351987, -161.07792184,   1.57784564e+08)],
     [( 709.59775655, -149.17797383,   1.57926514e+08)]]>

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

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

subeuvi = maps['EUVI'].submap(hpc_B[0], hpc_B[3])
subeuvi.peek()
../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_005.png

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)
    m.plot(axes=ax)
../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_006.png

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

Generated by Sphinx-Gallery