Overplotting HEK feature/event polygons on a map

How to overplot HEK outlines on a map.

import numpy as np

import matplotlib.pyplot as plt

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

import sunpy.map
import sunpy.data.sample
from sunpy.net import hek
from sunpy.time import parse_time
from sunpy.coordinates import frames
from sunpy.physics.differential_rotation import solar_rotate_coordinate

We start with the sample data

aia_map = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)

Look for coronal holes detected using the SPoCA feature recognition method:

hek_client = hek.HEKClient()
start_time = aia_map.date - TimeDelta(2*u.hour)
end_time = aia_map.date + TimeDelta(2*u.hour)
responses = hek_client.search(hek.attrs.Time(start_time, end_time),
                              hek.attrs.CH, hek.attrs.FRM.Name == 'SPoCA')

Let’s find the biggest coronal hole within 80 degrees north/south of the equator:

area = 0.0
for i, response in enumerate(responses):
    if response['area_atdiskcenter'] > area and np.abs(response['hgc_y']) < 80.0:
        area = response['area_atdiskcenter']
        response_index = i

Next let’s get the boundary of the coronal hole

ch = responses[response_index]
p1 = ch["hpc_boundcc"][9:-2]
p2 = p1.split(',')
p3 = [v.split(" ") for v in p2]
ch_date = parse_time(ch['event_starttime'])

The coronal hole was detected at different time than the AIA image was taken so we need to rotate it to the map observation time.

ch_boundary = SkyCoord(
    [(float(v[0]), float(v[1])) * u.arcsec for v in p3],
    obstime=ch_date,
    frame=frames.Helioprojective)
rotated_ch_boundary = solar_rotate_coordinate(ch_boundary, time=aia_map.date)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3.dev1+gece6362cd-py3.7-linux-x86_64.egg/sunpy/physics/differential_rotation.py:153: UserWarning: Using 'time' assumes an Earth-based observer.
  warnings.warn("Using 'time' assumes an Earth-based observer.")
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/coordinates/solar_system.py:364: ResourceWarning: unclosed file <_io.BufferedReader name='/home/docs/.astropy/cache/download/py3/d66c2153188958472fec67016b4b4e7a'>
  get_velocity=False)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/coordinates/solar_system.py:364: ResourceWarning: unclosed file <_io.BufferedReader name='/home/docs/.astropy/cache/download/py3/d66c2153188958472fec67016b4b4e7a'>
  get_velocity=False)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/coordinates/solar_system.py:364: ResourceWarning: unclosed file <_io.BufferedReader name='/home/docs/.astropy/cache/download/py3/d66c2153188958472fec67016b4b4e7a'>
  get_velocity=False)

Now let’s plot the rotated coronal hole boundary on the AIA map, and fill it with hatching.

fig = plt.figure()
ax = plt.subplot(projection=aia_map)
aia_map.plot(axes=ax)
ax.plot_coord(rotated_ch_boundary, color='c')
ax.set_title('{:s}\n{:s}'.format(aia_map.name, ch['frm_specificid']))
plt.colorbar()
plt.show()
../../../_images/sphx_glr_overplot_hek_polygon_001.png

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

Gallery generated by Sphinx-Gallery