Drawing the AIA limb on a STEREO EUVI image

In this example we use a STEREO-B and an SDO image to demonstrate how to overplot the limb as seen by AIA on an EUVI-B image. Then we overplot the AIA coordinate grid on the STEREO image.

import matplotlib.pyplot as plt

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

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

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)
result = Fido.search(wave, aia | stereo)

Let’s inspect the result and download the files.

print(result)
downloaded_files = Fido.fetch(result)
print(downloaded_files)

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


['/home/docs/sunpy/data/20210101_000615_n4eua.fts']
Errors:
(error(filepath_partial=functools.partial(<function VSOClient.mk_filename at 0x7f40d2b99cb0>, '/home/docs/sunpy/data/{file}', <QueryResponseRow index=0>
       Start Time               End Time        Source Instrument Wavelength [2] Provider  Physobs  Wavetype Extent Width Extent Length Extent Type   Size                              Info                                    fileid
                                                                     Angstrom                                                                        Mibyte
         object                  object          str3     str3       float64       str4      str9     str6       str4          str4         str8    float64                            str57                                    str24
----------------------- ----------------------- ------ ---------- -------------- -------- --------- -------- ------------ ------------- ----------- -------- --------------------------------------------------------- ------------------------
2021-01-01 00:06:05.000 2021-01-01 00:06:06.000    SDO        AIA 304.0 .. 304.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844 AIA level 1, 4096x4096 [2.901 exposure] [100.00 percentd] aia__lev1:304:1388534804), url='https://sdo7.nascom.nasa.gov/cgi-bin/drms_export.cgi?series=aia__lev1;compress=rice;record=304_1388534804-1388534804', exception=FailedDownload()
 https://sdo7.nascom.nasa.gov/cgi-bin/drms_export.cgi?series=aia__lev1;compress=rice;record=304_1388534804-1388534804 <ClientResponse(https://sdo7.nascom.nasa.gov/cgi-bin/drms_export.cgi?series=aia__lev1;compress=rice;record=304_1388534804-1388534804) [500 Internal Server Error]>
<CIMultiDictProxy('Date': 'Thu, 02 Dec 2021 16:16:03 GMT', 'Server': 'Apache', 'Strict-Transport-Security': 'max-age=31536000; includeSubDomains', 'Content-Length': '541', 'Connection': 'close', 'Content-Type': 'text/html; charset=iso-8859-1')>
))

Let’s create a dictionary with the two maps, which we crop to full disk.

maps = {m.detector: m.submap(SkyCoord([-1100, 1100], [-1100, 1100],
                                      unit=u.arcsec, frame=m.coordinate_frame))
        for m in sunpy.map.Map(downloaded_files)}
maps['AIA'].plot_settings['vmin'] = 0  # set the minimum plotted pixel value
Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/checkouts/stable/examples/units_and_coordinates/AIA_limb_STEREO.py", line 47, in <module>
    for m in sunpy.map.Map(downloaded_files)}
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/checkouts/stable/examples/units_and_coordinates/AIA_limb_STEREO.py", line 45, in <dictcomp>
    maps = {m.detector: m.submap(SkyCoord([-1100, 1100], [-1100, 1100],
  File "/home/docs/checkouts/readthedocs.org/user_builds/sunpy/envs/stable/lib/python3.7/site-packages/sunpy/map/mapbase.py", line 271, in __getitem__
    "The ability to index Map by physical"
NotImplementedError: The ability to index Map by physical coordinate is not yet implemented.

Now, let’s plot both maps, and we draw the limb as seen by AIA onto the EUVI image. We remove the part of the limb that is hidden because it is on the far side of the Sun from STEREO’s point of view.

fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1, 2, 1, projection=maps['AIA'])
maps['AIA'].plot(axes=ax1)
maps['AIA'].draw_limb()

ax2 = fig.add_subplot(1, 2, 2, projection=maps['EUVI'])
maps['EUVI'].plot(axes=ax2)
visible, hidden = maps['AIA'].draw_limb()
hidden.remove()

Let’s plot the helioprojective coordinate grid as seen by SDO on the STEREO image in a cropped view. Note that only those grid lines that intersect the edge of the plot will have corresponding ticks and tick labels.

fig = plt.figure()
ax = plt.subplot(projection=maps['EUVI'])

maps['EUVI'].plot()

# Crop the view using pixel coordinates
ax.set_xlim(500, 1300)
ax.set_ylim(100, 900)

# Shrink the plot slightly and move the title up to make room for new labels.
ax.set_position([0.1, 0.1, 0.8, 0.7])
ax.set_title(ax.get_title(), pad=45)

# Change the default grid labels and line properties.
stereo_x, stereo_y = ax.coords
stereo_x.set_axislabel("Helioprojective Longitude (STEREO B) [arcsec]")
stereo_y.set_axislabel("Helioprojective Latitude (STEREO B) [arcsec]")
ax.coords.grid(color='white', linewidth=1)

# Add a new coordinate overlay in the SDO frame.
overlay = ax.get_coords_overlay(maps['AIA'].coordinate_frame)
overlay.grid()

# Configure the grid:
x, y = overlay

# Wrap the longitude at 180 deg rather than the default 360.
x.set_coord_type('longitude', 180.)

# Set the tick spacing
x.set_ticks(spacing=250*u.arcsec)
y.set_ticks(spacing=250*u.arcsec)

# Set the ticks to be on the top and left axes.
x.set_ticks_position('tr')
y.set_ticks_position('tr')

# Change the defaults to arcseconds
x.set_major_formatter('s.s')
y.set_major_formatter('s.s')

# Add axes labels
x.set_axislabel("Helioprojective Longitude (SDO) [arcsec]")
y.set_axislabel("Helioprojective Latitude (SDO) [arcsec]")

plt.show()

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

Gallery generated by Sphinx-Gallery