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,  3.18file/s]SSL error in data received
protocol: <asyncio.sslproto.SSLProtocol object at 0x7f8781ddc160>
transport: <_SelectorSocketTransport fd=20 read=polling write=<idle, bufsize=0>>
Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/asyncio/sslproto.py", line 526, in data_received
    ssldata, appdata = self._sslpipe.feed_ssldata(data)
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/asyncio/sslproto.py", line 207, in feed_ssldata
    self._sslobj.unwrap()
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/ssl.py", line 767, in unwrap
    return self._sslobj.shutdown()
ssl.SSLError: [SSL: KRB5_S_INIT] application data after close notify (_ssl.c:2609)

Files Downloaded: 100%|##########| 2/2 [00:01<00:00,  1.65file/s]
['/home/docs/sunpy/data/aia_lev1_304a_2011_01_01t00_00_08_12z_image_lev1.fits', '/home/docs/sunpy/data/secchi_l0_b_img_euvi_20110101_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)}

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

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)
fig = plt.figure()
subaia.plot()
../../../_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-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:

<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, km)
    [[(422.21047096, -387.55701629, 1.57805382e+08)],
     [(642.59938241, -376.43551334, 1.57947332e+08)],
     [(513.1552792 , -161.09530724, 1.57784455e+08)],
     [(709.44158877, -149.2019928 , 1.57926363e+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)
    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])
fig = plt.figure()
plt.subplot(projection=subeuvi)
subeuvi.plot()
../../../_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)
    m.plot(axes=ax)
plt.show()
../../../_images/sphx_glr_SDO_to_STEREO_Coordinate_Conversion_006.png

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

Gallery generated by Sphinx-Gallery