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.coordinates
import sunpy.map
from sunpy.net import Fido
from sunpy.net import attrs as a
from sunpy.sun import constants

The first step is to download some data. We download an image from STEREO-A and an image from SDO, which are separated in longitude.

stereo = (a.Source('STEREO_A') &
          a.Instrument("EUVI") &
          a.Time('2021-01-01 00:06', '2021-01-01 00:07'))
aia = (a.Instrument.aia &
       a.Sample(24 * u.hour) &
       a.Time('2021-01-01 00:06', '2021-01-02 00:06'))
wave = a.Wavelength(30 * u.nm, 31 * u.nm)
res = Fido.search(wave, aia | stereo)
print(res)

Out:

Results from 2 Providers:

1 Results from the VSOClient:
Source: http://vso.stanford.edu/cgi-bin/search

       Start Time       ...
                        ...
----------------------- ...
2021-01-01 00:06:05.000 ...

1 Results from the VSOClient:
Source: http://vso.stanford.edu/cgi-bin/search

       Start Time               End Time        ...             Info
                                                ...
----------------------- ----------------------- ... ----------------------------
2021-01-01 00:06:15.000 2021-01-01 00:06:19.000 ... EUVI ;  ; NORMAL ; 2048x2048

Download the files.

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

Out:

['/home/docs/sunpy/data/aia_lev1_304a_2021_01_01t00_06_05_13z_image_lev1.fits', '/home/docs/sunpy/data/20210101_000615_n4eua.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)}
maps['AIA'].plot_settings['vmin'] = 0  # set the minimum plotted pixel value

Out:

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]

We will be transforming coordinates where the formation height of 304 A emission makes a difference, so we set the reference solar radius to 4 Mm above the solar surface (see Alissandrakis 2019).

for m in maps.values():
    m.meta['rsun_ref'] = (constants.radius + 4*u.Mm).to_value('m')

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)
AIA $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:05, EUVI-A $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:15

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

aia_bottom_left = SkyCoord(-800 * u.arcsec,
                           -300 * u.arcsec,
                           frame=maps['AIA'].coordinate_frame)
aia_top_right = SkyCoord(-600 * u.arcsec,
                         -50 * u.arcsec,
                         frame=maps['AIA'].coordinate_frame)

Plot a rectangle around the region we want to crop.

fig = plt.figure()
ax = fig.add_subplot(111, projection=maps['AIA'])
maps['AIA'].plot(axes=ax)
maps['AIA'].draw_quadrangle(aia_bottom_left, top_right=aia_top_right)
AIA $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:05

Out:

<astropy.visualization.wcsaxes.patches.Quadrangle object at 0x7f40bae60f90>

Create a submap of this area and draw an X at a specific feature of interest.

subaia = maps['AIA'].submap(aia_bottom_left, top_right=aia_top_right)
fig = plt.figure()
ax = fig.add_subplot(projection=subaia)
subaia.plot()

feature_aia = SkyCoord(-706 * u.arcsec,
                       -181 * u.arcsec,
                       frame=maps['AIA'].coordinate_frame)
ax.plot_coord(feature_aia, 'bx', fillstyle='none', markersize=20)
AIA $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:05

Out:

[<matplotlib.lines.Line2D object at 0x7f40baafe6d0>]

We can transform the coordinate of the feature to see its representation as seen by STEREO EUVI. The original coordinate did not contain a distance from the observer, so it is converted to a 3D coordinate assuming that it has a radius from Sun center equal to the reference solar radius of the coordinate frame, which we set earlier to be the formation height of 304 A emission.

print(feature_aia.transform_to(maps['EUVI'].coordinate_frame))

Out:

<SkyCoord (Helioprojective: obstime=2021-01-01T00:06:15.007, rsun=699700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2021-01-01T00:06:15.007, rsun=699700.0 km): (lon, lat, radius) in (deg, deg, m)
    (-56.45339738, 3.82395723, 1.43975156e+11)>): (Tx, Ty, distance) in (arcsec, arcsec, m)
    (155.86057277, -284.95344914, 1.43313023e+11)>

Now we can plot this box on both the AIA and EUVI images. Note that using draw_quadrangle() means that the plotted rectangle will be automatically warped appropriately to account for the different coordinate frames.

fig = plt.figure(figsize=(10, 4))

ax1 = fig.add_subplot(1, 2, 1, projection=maps['AIA'])
maps['AIA'].plot(axes=ax1)
maps['AIA'].draw_quadrangle(aia_bottom_left, top_right=aia_top_right, axes=ax1)

ax2 = fig.add_subplot(1, 2, 2, projection=maps['EUVI'])
maps['EUVI'].plot(axes=ax2)
maps['AIA'].draw_quadrangle(aia_bottom_left, top_right=aia_top_right, axes=ax2)
AIA $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:05, EUVI-A $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:15

Out:

<astropy.visualization.wcsaxes.patches.Quadrangle object at 0x7f40bad39a10>

We can now zoom in on the region in the EUVI image, and we also draw an X at the feature marked earlier. We do not need to explicitly transform the feature coordinate to the matching coordinate frame; that is performed automatically by plot_coord().

fig = plt.figure(figsize=(15, 5))

ax1 = fig.add_subplot(1, 2, 1, projection=subaia)
subaia.plot(axes=ax1)
ax1.plot_coord(feature_aia, 'bx', fillstyle='none', markersize=20)

subeuvi = maps['EUVI'].submap(aia_bottom_left, top_right=aia_top_right)

ax2 = fig.add_subplot(1, 2, 2, projection=subeuvi)
subeuvi.plot(axes=ax2)
maps['AIA'].draw_quadrangle(aia_bottom_left, top_right=aia_top_right, axes=ax2)
ax2.plot_coord(feature_aia, 'bx', fillstyle='none', markersize=20)

plt.show()
AIA $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:05, EUVI-A $304 \; \mathrm{\mathring{A}}$ 2021-01-01 00:06:15

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

Gallery generated by Sphinx-Gallery