Tracking an Active Region Across the Solar Disk#

This example demonstrates how to track an active region as it rotates across the solar disk and make cutouts around that active region at each time step to build a tracked datacube.

import matplotlib.pyplot as plt

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.visualization import ImageNormalize, SqrtStretch

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

First, let’s download a series of images in time using Fido. In this example, we will download a series of AIA 171 Å images observed over the course of half of a day at a cadence of 1 image every 1 hour.

query = Fido.search(a.Time('2018-05-30 00:00:00', '2018-05-30 12:00:00'),
                    a.Instrument.aia,
                    a.Wavelength(171*u.angstrom),
                    a.Sample(1*u.h))
print(query)
files = Fido.fetch(query)

Results from 1 Provider:

12 Results from the VSOClient:
Source: https://sdac.virtualsolar.org/cgi/search
Data retrieval status: https://docs.virtualsolar.org/wiki/VSOHealthReport
Total estimated size: 813.466 Mbyte

       Start Time               End Time        Source Instrument   Wavelength   Provider  Physobs  Wavetype Extent Width Extent Length Extent Type   Size
                                                                     Angstrom                                                                        Mibyte
----------------------- ----------------------- ------ ---------- -------------- -------- --------- -------- ------------ ------------- ----------- --------
2018-05-30 00:00:09.000 2018-05-30 00:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 01:00:09.000 2018-05-30 01:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 02:00:09.000 2018-05-30 02:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 03:00:09.000 2018-05-30 03:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 04:00:09.000 2018-05-30 04:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 05:00:09.000 2018-05-30 05:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 06:00:09.000 2018-05-30 06:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 07:00:09.000 2018-05-30 07:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 08:00:09.000 2018-05-30 08:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 09:00:09.000 2018-05-30 09:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 10:00:09.000 2018-05-30 10:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844
2018-05-30 11:00:09.000 2018-05-30 11:00:10.000    SDO        AIA 171.0 .. 171.0     JSOC intensity   NARROW         4096          4096    FULLDISK 64.64844

Now that we have a set of images in time, we can create a MapSequence to hold all of them and animate that sequence in time.

aia_sequence = sunpy.map.Map(files, sequence=True)

fig = plt.figure()
ax = fig.add_subplot(projection=aia_sequence[0])
norm = norm=ImageNormalize(vmin=0, vmax=3e3, stretch=SqrtStretch())
ani = aia_sequence.plot(axes=ax, norm=norm)

Next, let’s crop one of the maps in our sequence to the active region of interest.

corner = SkyCoord(Tx=-250*u.arcsec, Ty=0*u.arcsec, frame=aia_sequence[6].coordinate_frame)
cutout_map = aia_sequence[6].submap(corner, width=500*u.arcsec, height=500*u.arcsec)

fig = plt.figure()
ax = fig.add_subplot(projection=cutout_map)
cutout_map.plot(axes=ax)

$AIA $171 \; \mathrm{\mathring{A}}$ 2018-05-30 07:00:09$

<matplotlib.image.AxesImage object at 0x7fb194ba5ac0>

Now that we have our cutout around our active region, we can reproject each map in our sequence to the WCS of that cutout. Additionally, we will use the propagate_with_solar_surface context manager to adjust the field of view of the cutout with the rotation of the solar surface.

with propagate_with_solar_surface():
    aia_sequence_aligned = sunpy.map.Map([m.reproject_to(cutout_map.wcs) for m in aia_sequence], sequence=True)

Finally, we can animate our sequence of reprojected cutouts to confirm that we’ve tracked our active region of interest as a function of time across the solar disk.

fig = plt.figure()
ax = fig.add_subplot(projection=aia_sequence_aligned[0])
ani = aia_sequence_aligned.plot(axes=ax, cmap='sdoaia171', norm=norm)

plt.show()