aiapy: A SunPy affiliated package for analyzing data from the Atmospheric Imaging Assembly¶

The Atmospheric Imaging Assembly (AIA) instrument onboard the NASA Solar Dynamics Observatory (SDO) spacecraft has observed the full- disk of the Sun, nearly continuously, for the last ten years. It is one of three instruments on SDO, along with the Helioseismic and Magnetic Imager (HMI) and the Extreme Ultraviolet Variability Experiment (EVE). AIA is a narrowband imaging instrument comprised of four separate telescopes that collectively observe the full-disk of the Sun at ten different wavelengths: seven extreme ultraviolet (EUV) wavelengths, two far UV wavelengths, and one visible wavelength.

aiapy is a SunPy-affiliated Python package for analyzing calibrated (level 1) imaging data from AIA. It includes capabilities for aligning images between channels, deconvolving images with the instrument point-spread function (PSF), computing channel sensitivity as a function of wavelength, and correcting images for telescope degradation, among others. These capabilities are divided between the three subpackages of aiapy:

• aiapy.calibrate
• aiapy.psf
• aiapy.response

In this blog post, we will demonstrate the functions in each of these subpackages on a set of example level 1 AIA images.

Before we begin, we make a note of what versions of astropy, sunpy, and aiapy were used in the writing of this blog post.

Obtaining AIA Data¶

First, we will use the Fido client to fetch two AIA images from the VSO: one from the 171 Å channel and one from the 211 Å channel.

We can then create sunpy.map.Map objects from the downloaded files.

Note that both the sunpy and aiapy documentation pages include more detailed examples of how to obtain AIA data from the JSOC as well

• sunpy: Querying and Downloading Data from JSOC
• sunpy: Requesting cutouts of AIA images from the JSOC
• aiapy: Requesting specific AIA images from the JSOC

More complex queries (e.g. searching by exposure time or quality keyword) can be constructed by searching the JSOC directly via the drms package. See the last example listed above for an example of this.

PSF Deconvolution¶

AIA images are subject to convolution with the instrument PSF due to effects introduced by the filter mesh of the telescope and the CCD, among others. This has the effect of “blurring” the image. The PSF diffraction pattern may also be particularly noticeable during the impulsive phase of a flare where the intensity enhancement is very localized. To remove these artifacts, the PSF must be deconvolved from the level 1 image. This is because the PSF itself is defined on the pixel grid of the CCD.

We'll calculate the PSF for the 171 Å channel using the psf function in the aiapy.psf subpackage. The PSF model accounts for several different effects, including diffraction from the mesh grating of the filters, charge spreading, and jitter. This implementation of the PSF calculation follows Grigis et al (2012). Further details on how the PSF is calculated can be found in that document. Currently, this only works for $4096\times4096$ full frame images.

Note: The following cell will be significantly slower (many minutes) if you do not have a NVIDIA GPU and the cupy package installed. The GPU-enabled calculation has been benchmarked at ~5 s on an NVIDIA RTX 4090. Use of the GPU can be explicitly disabled by setting use_gpu=False.

We can then visualize the PSF. The diffraction “arms” extending from the center pixel can often be seen in flare observations due to the intense, small-scale brightening.

We can then use the calculated PSF to deconvolve our 171 Å image using the deconvolve function in aiapy.psf. This deconvolution algorithm is an implementation of the well-known Richardson-Lucy algorithm. If the PSF is not explicitly passed in, the PSF for that wavelength will be computed automatically. When deconvolving many images, it is highly recommend to precompute the PSF function for each desired wavelength.

NOTE: The deconvolution of the image can also be optionally accelerated with a NVIDIA GPU and cupy. On a Titan RTX, deconvolution of a full $4096\times4096$ image has been benchmarked at $\le1$ s. On a CPU, the equivalent operation should take ~25 s.

We can now zoom in on some loop structures that extend off the limb to examine the differences between our image before and after deconvolution.

Note that comparing the images on the left and right, deconvolving with the PSF removes the appearance of the diffraction pattern near the bright looptop. This also brings out more detail in the strands of the post-flare arcade.

We can see this more clearly by plotting the intensity of the original and deconvolved images across the top of the looptop.

Respiking Level 1 Images¶

The aiapy.calibrate module contains a number of useful functions, including the ability to "respike" images through the respike function. AIA level 1 images have been corrected for hot-pixels (commonly referred to as “spikes”) using an automated correction algorithm which detects them, removes them, and replaces the “holes” left in the image via interpolation. However, for certain research topics, this automated hot-pixel removal process may result in unwanted removal of bright points which may be physically meaningful. The respike function allows you to revert this removal by putting back all the removed pixel values. This corresponds to the aia_respike.pro IDL procedure as described in the SDO Analysis Guide.

We will create a "respiked" version of our level 1 171 Å image. This function will query the JSOC for the relevant "spikes" file which contains the locations (in pixel coordinates) and associate intensities for each of the hot pixels. The hot pixel values are then replaced in the image at the appropriate locations.

Looking at the two images, it is not obvious by visual inspection what the differences are between the despiked and respiked maps. We can use the fetch_spikes functions (which is called automatically by respike) to find the intensity values and the associated pixel coordinates

We can then use the pixel coordinates to find the "despiked" intensity values and plot the distributions of the "despiked" and "respiked" values.

Additionally, we can plot the coordinates of each of our hot pixels on our "respiked" map. We can convert the above pixel coordinates to world coordinates using the .pixel_to_world() method on our level 1 map. Note that we can also return the world coordinates from the fetch_spikes function automatically by passing in the as_coords=True keyword argument.

For a more detailed example of respiking level 1 images, see the respiking example in the aiapy example gallery.

Transforming Level 1 Images to Level 1.5¶

One of the most common operations performed on level 1 AIA images is "aiaprep," commonly performed with the aia_prep.pro procedure in IDL. This prep process, which promotes level 1 images to level 1.5, is a combination of two steps:

1. update pointing and
2. image registration.

The first step is accomplished through the update_pointing function in aiapy.calibrate. This function queries the 3-hourly master pointing table stored at the JSOC to update relevant FITS keywords.

We can verify this step by examining the reference pixel before and after the pointing update.

The second step is accomplished through register function in aiapy.calibrate. This function performs an affine transform on the level 1 pixel grid to remove the roll angle, scale the image to resolution of 0.6 arcseconds-per-pixel, and shift the center of the image to the center of the Sun. The original intensity values are then reinterpolated onto this transformed grid.

We can verify that the image has been transformed appropriately by comparing at the scale and rotation_matrix attributes before and after registration. Note that the rotation_matrix is now diagonalized, meaning that the world and pixel grids are aligned.

Note that we can also combine these two operations into a "prep" function to imitate aiaprep.pro.

Applying this "prep" procedure means that our 171 Å and 211 Å images are now interpolated to the same pixel grid.

Before performing any operations between two maps from different channels, one should always register each image such that the pixel grids of each image are aligned. For more information on transforming level 1 images to level 1.5, see the following example in the aiapy example gallery:

• Registering and aligning level 1 data

Degradation Correction¶

The performance of the AIA telescope is unfortunately degrading over time, leading to the resulting images becoming increasingly dim. We can correct for this by modeling the degradation over time using the degradation function. This function calculates the degradation factor $d(t)$ as,

where the leading coefficient is the ratio between the effective area at some epoch time $t_e$ and the initial time $t_0$. The term in the parentheses represents an interpolation from $t_e$ to the selected time $t$, where $\delta t=t-t_e$. $p_1,p_2,p_3$ as well as $A_{eff}(t_e)$ and $A_{eff}(t_0)$ are retrieved from the JSOC (by default) or a user-supplied table.

We can use this function to calculate the degradation of the 211 Å channel from the start of the SDO mission to now.

First, we'll query the degradation calibration table from the JSOC. Doing this ahead of time avoids repeated network calls when we want to repeatedly calculate the degradation.

We can then calculate the 211 Å degradation using the above correction table.

Note that we can also use the .date attribute on our 211 Å map to compute the degradation at the time of that map

The degradation calibration is routinely updated (e.g. using data from sounding rocket flights), with new versions added to the table returned from the JSOC. We can look at older versions of the calibration to see how the predicted degradation has changed using the calibration_version keyword argument. Note that the default version in aiapy v0.10.0 is 10.

And compare the two degradation curves as a function of time over the lifetime of the mission.

Additionally, aiapy.calibrate also includes the correct_degradation function which accepts a sunpy.map.Map, calculates the degradation factor, divides the image data by this degradation factor, and returns a new corrected map.

Computing Uncertainties¶

One of the newest additions to aiapy is the function for computing uncertainties on AIA intensities. This is done through the estimate_error function in aiapy.calibrate. estimate_error is an exact port of the SSWIDL function aia_bp_estimate_error.pro. This calculation includes errors from the shot noise, read noise, quantization, and onboard compression. The calculation can also optionally include contributions from the photometric calibration and errors in the atomic data.

We can pass in the (unitful) data array attached to our 171 Å map and compute the intensities.

To easily visualize the spatial distribution of these uncertainties, we can create a new sunpy.map.Map object from these errors.

We can also use the include_chianti flag to understand how the presence of atomic data errors from CHIANTI affect the uncertainty.

Similarly, we can look at how the photometric calibration, using the calibration from EVE, affects the uncertainties.

We can then look at the distribution of these uncertainties with these different combinations of keyword arguments to understand how these different options affect the resulting distribution of errors. Note that the inclusion of the atomic data results in the largest errors because there is an additional 25% (for 171 Å) uncertainty assumed on the input intensities.

Exposure Time Normalization¶

For normalizing images by the exposure time:

This option still properly accounts for the change in units but does not reset the exposure time.

Wavelength Response Functions¶

Lastly, the aiapy.response subpackage includes the Channel object for computing the wavelength response function of each channel as well as providing various pieces of metadata for each channel.

The wavelength response is a combination of the (wavelength-dependent) effective area and the gain of the telescope such that $R(\lambda,t) = A_{eff}(\lambda)G(\lambda)$. According to Section 2 of Boerner et al. (2012), the effective area is given by,

where $A_{geo}$ is the geometrical collecting area, $R_{P,S}$ are the reflectances of the primary and secondary mirrors, respectively, $T_{E,F}$ are the transmission efficiency of the entrance and focal-plane filters, respectively, $D$ is the contaminant transmittance of the optics and $Q$ is the quantum efficiency of the CCD.

We can compute the wavelength response for the 211 Å channel using the following method.

Each of the individual components of the effective area are also available as properties on the Channel object.

The .wavelength_response method also includes several keyword arguments for including the effects of degradation and EVE cross-calibration in the computation of the wavelength response function.

Note that these wavelength responses, when combined with atomic data (e.g. from CHIANTI), can be used to compute the temperature response functions of each AIA channel using Section 2.5 of Boerner et al. (2012). More information about this subpackage can be found in the Computing wavelength response functions example in the aiapy example gallery.

Resources¶

In addition to this blog post, there are also a number of resources for learning more about aiapy:

• Source code on GitLab
• Documentation
• Example Gallery
• Paper in the Journal of Open Source Software (JOSS)

aiapy is under active development and, just as with the core sunpy package and all other affilitated packages, contributions from the community are always welcome. To report a bug, request a feature, or simply ask a question, feel free to open an issue on GitLab.