Registering and aligning level 1 data

This example demonstrates how to convert AIA images to a common pointing, rescale them to a common plate scale, and remove the roll angle. This process is often referred to as “aia_prep” and the resulting data are typically referred to as level 1.5 data. In this example, we will demonstrate how to do this with aiapy. This corresponds to the aia_prep.pro procedure as described in the SDO Analysis Guide.

import sunpy.map

import aiapy.data.sample as sample_data
from aiapy.calibrate import normalize_exposure, register, update_pointing

Performing multi-wavelength analysis on level 1 data can be problematic as each of the AIA channels have slightly different spatial scales and roll angles. Furthermore, the estimates of the pointing keywords (CDELT1, CDELT2, CRPIX1, CRPIX2, CROTA2) may have been improved due to limb fitting procedures. The Joint Science Operations Center (JSOC) stores AIA image data and metadata separately; when users download AIA data, these two data types are combined to produce a FITS file. While metadata are continuously updated at JSOC, previously downloaded FITS files will not contain the most recent information.

Thus, before performing any multi-wavelength analyses, level 1 data should be updated to the most recent and accurate pointing and interpolated to a common grid in which the y-axis of the image is aligned with solar North.

First, let’s read a level 1 94 Å AIA image from the aiapy sample data into a Map object.

DEBUG: Reading /home/docs/.local/share/aiapy/aia_lev1_94a_2019_01_01t00_00_11_12z_image_lev1.fits [sunpy.map.map_factory]

The first step in this process is to update the metadata of the map to the most recent pointing using the aiapy.calibrate.update_pointing function. This function queries the JSOC for the most recent pointing information, updates the metadata, and returns a sunpy.map.Map with updated metadata.

m_updated_pointing = update_pointing(m)
DEBUG: HCRS->HGS [sunpy.coordinates.transformations]
DEBUG: ├─From: <HCRS Coordinate (obstime=2019-01-01T00:00:11.120): (ra, dec, distance) in (deg, deg, m)
       │           (100.86403557, 23.05574578, 1.47085026e+11)> [sunpy.coordinates.transformations]
DEBUG: ├─To  : <HeliographicStonyhurst Frame (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km)> [sunpy.coordinates.transformations]
DEBUG: └─Out : <HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
                   (-0.0148041, -2.97612324, 1.47085026e+11)> [sunpy.coordinates.transformations]
DEBUG: HPC->HPC [sunpy.coordinates.transformations]
DEBUG: ├─From: <Helioprojective Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
       │           (-0.0148041, -2.97612324, 1.47085026e+11)>): (Tx, Ty) in arcsec
       │           (0., 0.)> [sunpy.coordinates.transformations]
DEBUG: ├─To  : <Helioprojective Frame (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
       │           (-0.0148041, -2.97612324, 1.47085026e+11)>)> [sunpy.coordinates.transformations]
DEBUG: └─Out : <Helioprojective Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
                   (-0.0148041, -2.97612324, 1.47085026e+11)>): (Tx, Ty) in arcsec
                   (0., 0.)> [sunpy.coordinates.transformations]
DEBUG: Running following query: aia.master_pointing3h[2018.12.31_12:00:48_TAI-2019.01.01_12:00:48_TAI] [sunpy.net.jsoc.jsoc]
DEBUG: Requesting following keywords: **ALL** [sunpy.net.jsoc.jsoc]

If we take a look at the plate scale and rotation matrix of the map, we find that the scale is slightly off from the expected value of \(0.6''\) per pixel and that the rotation matrix has off-diagonal entries.

SpatialPair(axis1=<Quantity 0.600109 arcsec / pix>, axis2=<Quantity 0.600109 arcsec / pix>)
[[ 0.99999712  0.0024018 ]
 [-0.0024018   0.99999712]]

We can use the aiapy.calibrate.register function to scale the image to the \(0.6''\) per pixel and derotate the image such that the y-axis is aligned with solar North.

DEBUG: HCRS->HGS [sunpy.coordinates.transformations]
DEBUG: ├─From: <HCRS Coordinate (obstime=2019-01-01T00:00:11.120): (ra, dec, distance) in (deg, deg, m)
       │           (100.86403557, 23.05574578, 1.47085026e+11)> [sunpy.coordinates.transformations]
DEBUG: ├─To  : <HeliographicStonyhurst Frame (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km)> [sunpy.coordinates.transformations]
DEBUG: └─Out : <HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
                   (-0.0148041, -2.97612324, 1.47085026e+11)> [sunpy.coordinates.transformations]
DEBUG: HPC->HPC [sunpy.coordinates.transformations]
DEBUG: ├─From: <Helioprojective Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
       │           (-0.0148041, -2.97612324, 1.47085026e+11)>): (Tx, Ty) in arcsec
       │           (0., 0.)> [sunpy.coordinates.transformations]
DEBUG: ├─To  : <Helioprojective Frame (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
       │           (-0.0148041, -2.97612324, 1.47085026e+11)>)> [sunpy.coordinates.transformations]
DEBUG: └─Out : <Helioprojective Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
                   (-0.0148041, -2.97612324, 1.47085026e+11)>): (Tx, Ty) in arcsec
                   (0., 0.)> [sunpy.coordinates.transformations]
DEBUG: scipy rotating image: 2.353 s [sunpy.image.transform]
DEBUG: scipy clipping image: 0.036 s [sunpy.image.transform]

If we look again at the plate scale and rotation matrix, we should find that the plate scale in each direction is \(0.6''\) per pixel and that the rotation matrix is diagonalized. The image in m_registered is now a level 1.5 data product.

SpatialPair(axis1=<Quantity 0.6 arcsec / pix>, axis2=<Quantity 0.6 arcsec / pix>)
[[ 1.00000000e+00 -2.84766725e-20]
 [-4.62157541e-19  1.00000000e+00]]

Though it is not typically part of the level 1.5 “prep” data pipeline, it is also common to normalize the image to the exposure time such that the units of the image are DN / pixel / s.

m_normalized = normalize_exposure(m_registered)

Finally, we can plot the exposure-normalized map. Note that small negative pixel values are possible because CCD images were taken with a pedestal set at ~100 DN. This pedestal is then subtracted when the JSOC pipeline performs dark (+pedestal) subtraction and flatfielding to generate level 1 files.

AIA $94 \; \mathrm{\mathring{A}}$ 2019-01-01 00:00:11
DEBUG: HCRS->HGS [sunpy.coordinates.transformations]
DEBUG: ├─From: <HCRS Coordinate (obstime=2019-01-01T00:00:11.120): (ra, dec, distance) in (deg, deg, m)
       │           (100.86403557, 23.05574578, 1.47085026e+11)> [sunpy.coordinates.transformations]
DEBUG: ├─To  : <HeliographicStonyhurst Frame (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km)> [sunpy.coordinates.transformations]
DEBUG: └─Out : <HeliographicStonyhurst Coordinate (obstime=2019-01-01T00:00:11.120, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
                   (-0.0148041, -2.97612324, 1.47085026e+11)> [sunpy.coordinates.transformations]

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

Gallery generated by Sphinx-Gallery