6. A low-level guide to IRIS with Python

In some cases and for some people it is better to focus on simple and efficient methods that do not make use of higher-level abstractions or IRIS-specific packages such as IRISpy or Sunpy. The motivation could be two-fold: To detail a “barebones” interface when no other packages are available, or when performance is the goal, and secondly to provide a simple guide while IRISpy is not mature enough and its API is changing rapidly.

One can use the earlier chapters in this manual as a starting point for downloading (Looking for IRIS Level 2 data and downloading them) and reading IRIS data (Working with IRIS Level 2 data in Python) into a Python session, but once that is done (or if one wants to read the iris raster and SJI FITS files directly using Astropy) one can work with the data directly.

A practical guide has been written explaining how to read and analyse data from IRIS, using Python at this simple level is available Low level guide to IRIS using Python. The guide has the very nice feature that the tutorial itself is written as a Jupyter notebook, so it can be run in that environment, or if one prefers in e.g., an IPython session.

But we repeat a slightly shortened version of the tutorial here, where we instead of using Astropy FITS reader directly use the methods described in the first chapters of this manual.

6.1. Requirements

To follow the Low level guide to IRIS using Python guide only standard python scientific packages are necessary:

  • Numpy

  • Scipy

  • Matplotlib

  • and in addition as stated above Astropy for its FITS reader (and WCS module).

The guide is written for Python 3. It may not work under Python 2.X.

Further details and instructions are given in the tutorial itself.

6.2. Reading IRIS files

The IRIS level 2 FITS files are the standard science data product, and can be read with any standard FITS reader. Here we make use of the routines discussed in chapters 1, 2, and 3 of this manual, but using astropy.io.fits directly is also possible (though the raster data will then be transposed in x and y direction (which can be changed using a .T operator) compared to the examples show below).

For this part we will make use of a dataset from 2018-01-02. Please download the data files and unpack them in your work directory. The commands used here will assume you are at the directory with the IRIS FITS files.

We’ll start by setting up matplotlib and importing the level2 reader. Here we use the inline backend, but for interactive use you should probably use ipympl.

>>>  #%matplotlib inline
>>> import numpy as np
>>> import iris_lmsalpy.hcr2fits as hcr2fits
>>> import iris_lmsalpy.extract_irisL2data as extract_irisL2data
>>> import matplotlib.pyplot as plt
>>> # Set up some default matplotlib options
>>> plt.rcParams['figure.figsize'] = [10, 6]
>>> plt.rcParams['xtick.direction'] = 'out'
>>> plt.rcParams['image.origin'] = 'lower'
>>> plt.rcParams['image.cmap'] = 'viridis'

6.2.1. Reading level 2 SJI files

Let us open a SJI file from the 1400 filter.

>>> SJI_filename = 'iris_l2_20180102_153155_3610108077_SJI_1400_t000.fits'
>>> sji = extract_irisL2data.load(SJI_filename)
The provided file is a SJI IRIS Level 2 data file containing SJI_1400 data.

Extracting information from file iris_l2_20180102_153155_3610108077_SJI_1400_t000.fits...

Available data with size Y x X x Image are stored in a windows labeled as:

-------------------------------------------------------------
Index --- Window label --- Y x X x Im --- Spectral range [AA]
-------------------------------------------------------------
  0         SJI_1400       548x845x80     1380.00 - 1420.00
-------------------------------------------------------------

Observation description:  Very large dense 320-step raster 105.3x175 320s   Deep x 8 Spatial x


The SJI data are passed to the output variable
(e.g iris_data) as a dictionary with the following entries:

- iris_data.SJI['SJI_1400']

The data itself and some metadata is contained in the sji object as shown earlier, but we can extract the main header (or auxiliary headers using the extent keyword; see ITN 26 on structure of level 2 files).

We can have a look at the main header, which is in the first HDU:

>>> hd = extract_irisL2data.only_header(SJI_filename, verbose = True)
Reading file iris_l2_20180102_153155_3610108077_SJI_1400_t000.fits...
Filename: iris_l2_20180102_153155_3610108077_SJI_1400_t000.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
0  PRIMARY       1 PrimaryHDU     162   (845, 548, 80)   int16 (rescales to float32)
1                1 ImageHDU        38   (31, 80)   float64
2                1 TableHDU        33   80R x 5C   [A10, A10, A3, A66, A55]

Observation description:  Very large dense 320-step raster 105.3x175 320s   Deep x 8 Spatial x

Extension No. 1 stores data and header of SJI_1400: 1380.00 - 1420.00 AA

To get the main header use :
hdr = extract_iris2level.only_header(filename)

To get header corresponding to data of SJI_1400 use :
hdr = extract_irisL2data.only_header(filename, extension = 0)

To get the data of SJI_1400 use :
data = extract_irisL2data.only_data(filename, extension = 0, memmap = True)

Returning header for extension No. 0

>>> hd
SIMPLE  =                    T / Written by IDL:  Tue Mar 20 23:39:27 2018
BITPIX  =                   16 / Number of bits per data pixel
NAXIS   =                    3 / Number of data axes
NAXIS1  =                  845 /
NAXIS2  =                  548 /
NAXIS3  =                   80 /
EXTEND  =                    T / FITS data may contain extensions
DATE    = '2018-03-21'         / Creation UTC (CCCC-MM-DD) date of FITS header
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode 2001A&A...376..359H
TELESCOP= 'IRIS    '           /
INSTRUME= 'SJI     '           /
DATA_LEV=              2.00000 /
LVL_NUM =              2.00000 /
VER_RF2 = 'L12-2017-08-04'     /
DATE_RF2= '2018-03-21T06:00:07.640' /
DATA_SRC=              1.51000 /
ORIGIN  = 'SDO     '           /
BLD_VERS= 'V9R1X   '           /
LUTID   =              0.00000 /
OBSID   = '3610108077'         /
OBS_DESC= 'Very large dense 320-step raster 105.3x175 320s   Deep x 8 Spatial x'
OBSLABEL= '        '           /
OBSTITLE= '        '           /
DATE_OBS= '2018-01-02T15:32:23.340' /
DATE_END= '2018-01-02T16:20:52.740' /
STARTOBS= '2018-01-02T15:31:55.700' /
ENDOBS  = '2018-01-02T16:21:01.960' /
OBSREP  =                    1 /
CAMERA  =                    2 /
STATUS  = 'Quicklook'          /
BTYPE   = 'Intensity'          /
BUNIT   = 'Corrected DN'       /
BSCALE  =                 0.25 /
BZERO   =                 7992 /
HLZ     =                    0 /
SAA     = '           0'       /
SAT_ROT =          0.000227430 /
AECNOBS =                    0 /
AECNRAS =                    0 /
DSUN_OBS=          1.47094E+11 /
IAECEVFL= 'NO      '           /
IAECFLAG= 'NO      '           /
IAECFLFL= 'NO      '           /
TR_MODE = '        '           /
FOVY    =              182.320 /
FOVX    =              281.132 /
XCEN    =             -329.321 /
YCEN    =              286.652 /
SUMSPTRL=                    2 /
SUMSPAT =                    2 /
EXPTIME =              8.00006 /
EXPMIN  =              8.00000 /
EXPMAX  =              8.00011 /
DATAMEAN=              27.0596 /
DATARMS =              21.6273 /
DATAMEDN=              21.8717 /
DATAMIN =             -3219.22 /
DATAMAX =              25706.2 /
DATAVALS=             20441608 /
MISSVALS=             16603192 /
NSATPIX =                    0 /
NSPIKES =                    0 /
TOTVALS =             37044800 /
PERCENTD=              55.1808 /
DATASKEW=              7.27218 /
DATAKURT=              234.259 /
DATAP01 =              8.83011 /
DATAP10 =              12.9206 /
DATAP25 =              16.3524 /
DATAP75 =              31.4368 /
DATAP90 =              46.9377 /
DATAP95 =              61.2195 /
DATAP98 =              81.1593 /
DATAP99 =              97.7437 /
NEXP_PRP=              1.00000 /
NEXP    =                   80 /
NEXPOBS =                  960 /
NRASTERP=                   80 /
RASTYPDX=                    1 /
RASTYPNX=                    1 /
RASRPT  =                    1 /
RASNRPT =                    1 /
CADPL_AV=              36.8290 /
CADPL_DV=              0.00000 /
CADEX_AV=              36.8284 /
CADEX_DV=            0.0124631 /
MISSOBS =                    1 /
MISSRAS =                    0 /
IPRPVER =        1.96000003815 /
IPRPPDBV=        10.0000000000 /
IPRPDVER=             20130925 /
IPRPBVER=             20180310 /
PC1_1   =       0.999970614910 /
PC1_2   =      0.0112785818055 /
PC2_1   =     -0.0112785818055 /
PC2_2   =       0.999970614910 /
PC3_1   =        0.00000000000 /
PC3_2   =        0.00000000000 /
NWIN    =                    1 /
TDET1   = 'SJI     '           /
TDESC1  = 'SJI_1400'           /
TWAVE1  =        1400.00000000 /
TWMIN1  =        1380.00000000 /
TWMAX1  =        1420.00000000 /
TDMEAN1 =              27.0596 /
TDRMS1  =              21.6273 /
TDMEDN1 =              21.8717 /
TDMIN1  =             -3219.22 /
TDMAX1  =              25706.2 /
TDVALS1 =             20441608 /
TMISSV1 =             16603192 /
TSATPX1 =                    0 /
TSPIKE1 =                    0 /
TTOTV1  =             37044800 /
TPCTD1  =              55.1808 /
TDSKEW1 =              7.27218 /
TDKURT1 =              234.259 /
TDP01_1 =              8.83011 /
TDP10_1 =              12.9206 /
TDP25_1 =              16.3524 /
TDP75_1 =              31.4368 /
TDP90_1 =              46.9377 /
TDP95_1 =              61.2195 /
TDP98_1 =              81.1593 /
TDP99_1 =              97.7437 /
TSR1    =                    1 /
TER1    =                  548 /
TSC1    =                    4 /
TEC1    =                  512 /
IPRPFV1 =                  441 /
IPRPGV1 =                    4 /
IPRPPV1 =                  441 /
CDELT1  =             0.332700 /
CDELT2  =             0.332700 /
CDELT3  =              36.8284 /
CRPIX1  =              423.000 /
CRPIX2  =              274.500 /
CRPIX3  =              40.0000 /
CRVAL1  =             -329.321 /
CRVAL2  =              286.652 /
CRVAL3  =              1463.91 /
CTYPE1  = 'HPLN-TAN'           /
CTYPE2  = 'HPLT-TAN'           /
CTYPE3  = 'Time    '           /
CUNIT1  = 'arcsec  '           /
CUNIT2  = 'arcsec  '           /
CUNIT3  = 'seconds '           /
KEYWDDOC= 'http://www.lmsal.com/iris_science/irisfitskeywords.pdf' /
HISTORY iris_prep  Set 1 saturated pixels to Inf
HISTORY iris_prep  Dark v20130925; T=[-59.1,-59.0,-58.3,-59.3,1.4,2.7,2.7,2.7,2.
HISTORY iris_prep  Flat fielded with recnum     440
HISTORY iris_prep  Set permanently bad pixels to 0 prior to warping
HISTORY iris_prep_geowave_roi  ran with rec_num 4
HISTORY iris_prep_geowave_roi  boxwarp set to 1
HISTORY iris_prep_geowave_roi  updated WCS parameters with iris_isp2wcs
HISTORY iris_prep  Used iris_mk_pointdb ver 10
HISTORY iris_prep       Fiducial midpoint shift X,Y [pix]:      0.22    -0.17
HISTORY iris_prep  Used INF_POLY_2D for warping
HISTORY iris_prep VERSION:  1.96
HISTORY iris_prep   ran on 20180321_063639
HISTORY level2  Version L12-2017-08-04

Header entries can be accessed like a dictionary. For example, for the start of the slit-jaw sequence:

>>> hd['STARTOBS']
'2018-01-02T15:31:55.700'

The data are saved under sji.SJI['SJI_1400'].data. You can see that NAXIS1=845 (spatial x), NAXIS2=548, (spatial y), and NAXIS3=80 (time). Because in Python the arrays are in C order, the dimension order will be reversed:

>>> sji.SJI['SJI_1400'].data.shape
(80, 548, 845)

We can plot the first timestep to see what it looks like:

>>> plt.imshow(sji.SJI['SJI_1400'].data[...,0])
>>> plt.show()
_images/reading_iris_11_1.png

The image looks washed out because of the automatic scaling. We can get a better scaling by adjusting the maximum and minimum range:

>>> plt.imshow(sji.SJI['SJI_1400'].data[...,0],vmin=0, vmax=100)
>>> plt.show()
_images/reading_iris_13_1.png

In some cases, loading the whole data into memory may not be practical or too slow. For these, it is possible to use an approach based on numpy.memmap that allows us to load only parts of the data, and only when necessary. Refer to the first chapters of the manual, or the original tutorial on low-level python processing for details.

But let us show another method for scaling the data based on histogram optimization.

>>> import iris_lmsalpy.im_view as im

>>> lev=im.histo_opt(sji.SJI['SJI_1400'].data[...,0])
-2.00e+02/3.24e+02 image min/image max
5.71e+00/1.40e+02 histo_opt min/histo_opt max

>>> plt.imshow(sji.SJI['SJI_1400'].data[...,0],vmin=lev[0], vmax=lev[1])
>>> plt.show()
_images/reading_iris_15_2.png

6.2.1.1. Coordinates

So far, we have used plt.imshow directly, which shows the data in pixel values. Coordinate information is available in the SJI header, with the World Coordinate System (WCS) keywords. The astropy package has an wcs module that allows plotting the image with the proper solar (x, y) coordinates.

First we get the WCS data from the SJI header:

>>> from astropy.wcs import WCS
>>> wcs = WCS(hd)
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'seconds' -> 's'. [astropy.wcs.wcs]

Inspecting this object we see:

>>> wcs
WCS Keywords

Number of WCS axes: 3
CTYPE : 'HPLN-TAN'  'HPLT-TAN'  'Time'
CRVAL : -0.09147805555555556  0.07962555555555555  1463.91
CRPIX : 423.0  274.5  40.0
PC1_1 PC1_2 PC1_3  : 0.99997061491  0.0112785818055  0.0
PC2_1 PC2_2 PC2_3  : -0.0112785818055  0.99997061491  0.0
PC3_1 PC3_2 PC3_3  : 0.0  0.0  1.0
CDELT : 9.241666666666666e-05  9.241666666666666e-05  36.8284
NAXIS : 845  548  80

We can see that the first two axes have type of ‘HPLN’ and ‘HPLT’, meaning helioprojective-cartesian longitude and latitude respectively, or solar X and Y. The last axis is time in seconds.

To plot the coordinate axes in matplotlib, we do the following:

>>> ax = plt.subplot(projection=wcs.dropaxis(-1))
>>> ax.imshow(sji.SJI['SJI_1400'].data[...,0], vmin=lev[0], vmax=lev[1])
>>> plt.show()
_images/reading_iris_21_1.png

We do wcs.dropaxis(-1) because we do not want to plot the time coordinate.

We can also overplot the coordinates grid by calling ax.grid. For the final figure with grid, we do:

>>> ax = plt.subplot(projection=wcs.dropaxis(-1))
>>> ax.grid(color='w', ls=':')
>>> ax.imshow(sji.SJI['SJI_1400'].data[...,0], vmin=lev[0], vmax=lev[1])
>>> plt.show()
_images/reading_iris_24_1.png

You’ll notice that the solar coordinate grid is slightly tilted from the image axes. This is normal. With different IRIS roll angles, this difference will be even more obvious.

6.2.2. Reading level 2 raster files

IRIS level 2 spectrograph files have a different structure from the SJI files. The main header is similar, but the data are all saved in extension HDUs. The reason for this is that spectrograph level 2 files are organised by spectral windows, which may differ for various observations.

Using the raster file from our example above, we can open the files with the methods described above and in previous chapters:

>>> raster_filename = 'iris_l2_20180102_153155_3610108077_raster_t000_r00000.fits'
>>> sp = extract_irisL2data.load(raster_filename)
The provided file is a raster IRIS Level 2 data file.

Extracting information from file /Users/asainz/iris_l2_20180102_153155_3610108077_raster_t000_r00000.fits...
Extracting information from file iris_l2_20180102_153155_3610108077_raster_t000_r00000.fits...

Available data with size Y x X x Wavelength are stored in windows labeled as:

--------------------------------------------------------------------
Index --- Window label --- Y x X x WL --- Spectral range [AA] (band)
--------------------------------------------------------------------
  0      C II 1336         548x320x335     1331.77 - 1340.44  (FUV)
  1      O I 1356          548x320x414     1346.75 - 1357.47  (FUV)
  2      Si IV 1394        548x320x195     1391.20 - 1396.14  (FUV)
  3      Si IV 1403        548x320x338     1398.12 - 1406.70  (FUV)
  4      2832              548x320x60      2831.34 - 2834.34  (NUV)
  5      2814              548x320x76      2812.65 - 2816.47  (NUV)
  6      Mg II k 2796      548x320x380     2790.65 - 2809.95  (NUV)
--------------------------------------------------------------------

Observation description:  Very large dense 320-step raster 105.3x175 320s   Deep x 8 Spatial x


The selected data are passed to the output variable
(e.g iris_data) as a dictionary with the following entries:

 - iris_data.raster['Mg II k 2796']

(The same comment made about memmap applies as for the SJI files, and is more important in the raster files because they tend to be larger.)

The first FITS HDU has the main header. We can use it to find out what kind of raster this is:

>>> hd = extract_irisL2data.getheader(raster_filename)
>>> hd['OBS_DESC']
'Very large dense 320-step raster 105.3x175 320s   Deep x 8 Spatial x'

The number of FITS units in sp is the number of spectral windows plus 3 (with additional metadata). In this case we have 10 FITS units, so 7 spectral windows. We can also look at the NWIN keyword in the main header:

Let us work with the Mg II k 2796 window, which in the fits file is found in extension 7 (= window index + 1). The shape of this data array is as follows:

>>> sp.raster['Mg II k 2796'].data.shape
(548, 320, 380)

The first axis is the y coordinate (space along slit), the second axis the number of raster positions (or time if sit and stare), and the last axis is the wavelength. We can make use of astropy.wcs to use the WCS header information to convert from pixels to coordinates. For example, for obtaining the wavelength scale we do:

>>> from astropy.wcs import WCS
>>> wcs=WCS(extract_irisL2data.getheader(raster_filename,ext=7))
>>> m_to_nm = 1e9  # convert wavelength to nm
>>> nwave = nwave=sp.raster['Mg II k 2796'].data.shape[2]
>>> wavelength = wcs.all_pix2world(np.arange(nwave), [0.], [0.], 0)[0] * m_to_nm

And we can now plot an example spectrum, say pixel [100, 200]:

>>> plt.plot(wavelength,sp.raster['Mg II k 2796'].data[200,100])
>>> plt.xlabel("Wavelength (nm)")
>>> plt.ylabel("Intensity (DN)")
>>> plt.show()
_images/reading_iris_36_1.png

The points at the edges of the plot have a special value of -200 to denote regions not recorded by the detector.

6.2.2.1. Spectroheliograms

Because in this case we have a 320-step raster, we can use it to build a spectroheliogram at a fixed wavelength position. Suppose we want to do this for the core of the Mg II k line. Its vacuum wavelength is 279.6351 nm. Let’s find the wavelength point closest to this:

>>> mg_index = np.argmin(np.abs(wavelength - 279.6351))

Now we can plot the data with an imshow:

>>> plt.imshow(sp.raster['Mg II k 2796'].data[...,mg_index], vmin=100, vmax=2000)
>>> plt.show()
_images/reading_iris_40_1.png

Let’s go a step farther and, as in the SJI example, make use of astropy.wcs to get the correct axes when plotting:

>>> figure = plt.figure(figsize=(6, 10))
>>> ax = plt.subplot(projection=wcs.dropaxis(0), slices=('y', 'x'))
>>> ax.set_xlabel("Solar X")
>>> ax.set_ylabel("Solar Y")
>>> ax.imshow(sp.raster['Mg II k 2796'].data[...,mg_index], vmin=100, vmax=2000)
>>> plt.show()
_images/reading_iris_42_1.png

You can compare this image with the ones from the SJI above.

You’ll notice we had to swap the axes in the WCS projection call using slices=('y', 'x') since the rasters are stored in a transposed manner in the FITS files, this has been corrected in the data during extraction.

6.2.2.2. Time of observations and exposure times

As in the SJI level 2 files, the exposure times are saved both in the main header and in the auxiliary metadata (as a function of time). In most cases, the exposure times are fixed for all scans in a raster. However, when automatic exposure compensation (AEC) is switched on and there is a very energetic event (e.g. a flare), IRIS will automatically use a lower exposure time to prevent saturation in the detectors. You can see the default exposure time, as well as the maximum and minimum exposure times in the header:

>>> print(hd['EXPTIME'], hd['EXPMIN'], hd['EXPMAX'])
7.99914 7.99909 7.99925

In this case these are all the same, meaning that no change in exposure time took place.

If the exposure time varies, you can get the time-dependent exposure times in seconds from the auxiliary metadata, second to last HDU in the file (after the last line window;`hd[‘NWIN’]+1`.

>>> aux_hd = extract_irisL2data.only_header(raster_filename,extension=hd['NWIN']+1)
>>> aux = extract_irisL2data.only_data(raster_filename,extension=hd['NWIN']+1)
>>> fuv_exptime = aux[:,aux_hd['EXPTIMEF']]
>>> nuv_exptime = aux[:,aux_hd['EXPTIMEN']]

The times of each raster scan are also saved in the auxiliary metadata. These are the time in seconds since the start of the observations, given by

>>> hd['DATE_OBS']
'2018-01-02T15:31:55.870'

To get an array with the full date/times for individual timesteps, we can make use of numpy.datetime64:

>>> time_diff = aux[:, aux_hd['TIME']]
>>> times = np.datetime64(hd['DATE_OBS']) + time_diff * np.timedelta64(1, 's')

6.2.3. Reading level 3 files

IRIS level 3 FITS files combine several spectrograph raster files into two files. They combine the temporal domain into raster data, resulting in 4-D arrays of dimensions (nx, ny, nwave, ntime). The main application of level 3 files is the CRISPEX data analysis software. Level 3 files are not available for download, and the users are expected to create them by combining level 2 spectrograph files.

We will include a descriptions of how to read the level 3 data at a later date. A description using the astropy’s standard fits reader can be found in the low-level python processing tutorial.

6.3. Some worked examples

We present some worked examples to show how one can access and use auxiliary data contained in the IRIS level 2 FITS files in data analysis. To run these examples you will need to download the IRIS data, and make sure that you are in the same directory as the files to be able to run this code. We start by importing the relevant packages and setting out some plotting options:

>>> #%matplotlib inline
>>> import numpy as np
>>> import iris_lmsalpy.extract_irisL2data as extract_irisL2data
>>> import matplotlib.pyplot as plt
>>> from astropy.wcs import WCS
>>>  # Set up some default matplotlib options
>>> plt.rcParams['figure.figsize'] = [10, 6]
>>> plt.rcParams['xtick.direction'] = 'out'
>>> plt.rcParams['image.origin'] = 'lower'
>>> plt.rcParams['image.cmap'] = 'viridis'

6.3.1. Mg II Dopplergrams

In this example we are going to produce a Dopplergram for the Mg II k line from an IRIS 400-step raster. The Dopplergram is obtained by subtracting the intensities at symmetrical velocity shifts from the line core (e.g. ±50 km/s). For this kind of analysis we need a consistent wavelength calibration for each step of the raster.

This very large dense raster took more than three hours to complete the 400 scans (30 s exposures), which means that the spacecraft’s orbital velocity changes during the observations. This means that any precise wavelength calibration will need to correct for those shifts.

Start by downloading an IRIS dataset from 2014 July 8: follow this link and download the raster file (726 Mb). Download it to a directory of your choice. Untar it, e.g. on a UNIX system:

$ tar zxvf iris_l2_20140708_114109_3824262996_raster.tar.gz

Let’s open the file (by default the Mg II k line window is downloaded, if another window was desired then the window_info keyword needs to be set):

>>> raster_filename = 'iris_l2_20140708_114109_3824262996_raster_t000_r00000.fits'
>>> sp = extract_irisL2data.load(raster_filename)
>>> hd = extract_irisL2data.getheader(raster_filename)
The provided file is a raster IRIS Level 2 data file.

Extracting information from file iris_l2_20140708_114109_3824262996_raster_t000_r00000.fits...

Available data with size Y x X x Wavelength are stored in windows labeled as:

------------------------------------------------------------
Index in the output variable --- Window label --- Y x X x WL
------------------------------------------------------------
             0               C II 1336       1095x400x189
             1               Fe XII 1349     1095x400x125
             2               O I 1356        1095x400x172
             3               Si IV 1394      1095x400x209
             4               Si IV 1403      1095x400x304
             5               2832            1095x400x120
             6               2814            1095x400x153
             7               Mg II k 2796    1095x400x537
------------------------------------------------------------


The selected data are passed to the output variable
(e.g iris_data) as a dictionary with the following entries:

 - iris_data.raster['Mg II k 2796']

Note that this information is retrievable from the header:

>>> print('Window. Name      : wave start - wave end\n')
>>> for i in range(hd['NWIN']):
>>>     win = str(i + 1)
>>>     print('{0}. {1:15}: {2:.2f} - {3:.2f} Å'
>>>           ''.format(win, hd['TDESC' + win], hd['TWMIN' + win], hd['TWMAX' + win]))
Window. Name      : wave start - wave end

1. C II 1336      : 1332.70 - 1337.58 Å
2. Fe XII 1349    : 1347.68 - 1350.90 Å
3. O I 1356       : 1352.25 - 1356.69 Å
4. Si IV 1394     : 1390.90 - 1396.19 Å
5. Si IV 1403     : 1398.63 - 1406.34 Å
6. 2832           : 2831.23 - 2834.26 Å
7. 2814           : 2812.54 - 2816.41 Å
8. Mg II k 2796   : 2792.98 - 2806.63 Å

We see that Mg II is in window 8, so let’s load those data and get the WCS information to construct the wavelength arrays and get coordinates:

>>> data=sp.raster['Mg II k 2796'].data
>>> wcs=WCS(extract_irisL2data.only_header(raster_filename, extension = 8))

And let’s calculate the wavelength array:

>>>  m_to_nm = 1e9  # convert wavelength to nm
>>>  wave = wcs.all_pix2world(np.arange(wcs._naxis[0]), [0.], [0.], 1)[0] * m_to_nm

We can see how the the spatially averaged spectrum looks like:

>>> plt.plot(wave, data.mean((0, 1)))
>>> plt.show()
_images/tutorials_11_1.png

To better understand the orbital velocity problem let us look at how the line intensity varies for a strong Mn I line at around 280.2 nm, in between the Mg II k and h lines. For this dataset, the line core of this line falls around index 350. To plot a spectroheliogram in the correct orientation we will transpose the data:

>>> fig = plt.figure(figsize=(6, 10))
>>> import iris_lmsalpy.im_view as im
>>> lev = im.histo_opt(data[...,350])
>>> plt.imshow(data[..., 350], vmin=lev[0], vmax=lev[1], aspect=0.5)
>>> plt.show()
 -2.00e+02/5.39e+02 image min/image max
 -1.06e+00/1.00e+02 histo_opt min/histo_opt max
_images/tutorials_13_2.png

You can see a regular bright-dark pattern along the x axis, an indication that its intensities are not taken at the same position in the line because of wavelength shifts. The shifts are caused by the orbital velocity changes, and we can find these in the auxiliary metadata which are to be found in the extension past the “last” window in the data set:

>>>aux=extract_irisL2data.only_data(raster_filename,extension=hd['NWIN']+1)
>>> aux_hd=extract_irisL2data.only_header(raster_filename,extension=hd['NWIN']+1)
>>> v_obs=aux[:,aux_hd['OBS_VRIX']]
>>> v_obs /= 1000.  # convert to km/s
>>> plt.ylabel("Orbital velocity (km/s)")
>>> plt.xlabel("Scan number")
>>> plt.plot(v_obs)
>>> plt.show()
_images/tutorials_15_1.png

To look at intensities at any given scan we only need to subtract this velocity shift from the wavelength scale, but to look at the whole image at a given wavelength we must interpolate the original data to take this shift into account. Here is a way to do it (note that array dimensions apply to this specific set only!):

>>> from scipy.constants import c
>>> from scipy.interpolate import interp1d
>>> c_kms = c / 1000.
>>> wave_shift = - v_obs * wave[350] / (c_kms)
>>> # linear interpolation in wavelength, for each scan
>>> for i in range(hd['NRASTERP']):
>>>     tmp = interp1d(wave - wave_shift[i], data[:,i,:], bounds_error=False)
>>>     data[:,i,:] = tmp(wave)

And now we can plot the shifted data to see that the large scale shifts have disappeared:

>>> fig = plt.figure(figsize=(6, 10))
>>> plt.imshow(data[..., 350], vmin=lev[0], vmax=lev[1], aspect=0.5)
>>> plt.show()
_images/tutorials_19_1.png

Some residual shift remains, but we will not correct for it here. A more elaborate correction can be obtained by the IDL routine iris_prep_wavecorr_l2, but this has not yet been ported to Python see the IDL version of this tutorial for more details.

We can use the calibrated data for example to calculate Dopplergrams. A Dopplergram is here defined as the difference between the intensities at two wavelength positions at the same (and opposite) distance from the line core. For example, at +/- 50 km/s from the Mg II k3 core. To do this, let us first calculate a velocity scale for the k line and find the indices of the -50 and +50 km/s velocity positions (here using the convention of negative velocities for up flows):

>>> mg_k_centre = 279.6351  # in nm
>>> pos = 50  # in km/s around line centre
>>> velocity =  (wave - mg_k_centre) * c_kms / mg_k_centre
>>> index_p = np.argmin(np.abs(velocity - pos))
>>> index_m = np.argmin(np.abs(velocity + pos))
>>> doppl = data[..., index_m] - data[..., index_p]

And now we can plot this as before (intensity units are again arbitrary because of the unscaled DNs):

>>> fig = plt.figure(figsize=(6, 10))
>>> lev = im.histo_opt(doppl)
>>> plt.imshow(doppl, cmap='gist_gray', aspect=0.5,
>>>            vmin=-1.*np.max(np.abs(lev)), vmax=np.max(np.abs(lev)))
>>> plt.show()
-1.79e+03/7.48e+02 image min/image max
-6.32e+01/2.53e+02 histo_opt min/histo_opt max
_images/tutorials_23_2.png

6.3.2. Time series analysis

As a final example we are going to work with spectra and slit-jaw images to study dynamical phenomena. The subject of this example is umbral flashes.

Start by downloading this IRIS dataset from 2013 September 2. Download the 1400 slit-jaw file and the raster file (about 900 Mb in total) to a directory of your choice. Unzip all the files, e.g. on a UNIX system:

$ gunzip *fits.gz
$ tar zxvf *.tar.gz

Now we load the spectrograph FITS file and get the header:

>>> raster_filename='iris_l2_20130902_163935_4000255147_raster_t000_r00000.fits'
>>> lines = extract_irisL2data.show_lines(raster_filename)
>>> hd = extract_irisL2data.only_header(raster_filename)
Extracting information from file iris_l2_20130902_163935_4000255147_raster_t000_r00000.fits...

Available data with size Y x X x Wavelength are stored in windows labeled as:

------------------------------------------------------------
Index in the output variable --- Window label --- Y x X x WL
------------------------------------------------------------
             0               C II 1336       417x1600x229
             1               Fe XII 1349     417x1600x53
             2               Cl I 1352       417x1600x36
             3               O I 1356        417x1600x53
             4               Si IV 1394      417x1600x150
             5               Si IV 1403      417x1600x262
             6               2786            417x1600x171
             7               Mg II k 2796    417x1600x555
             8               2830            417x1600x42
------------------------------------------------------------

We can now get the data and wavelength for the Mg II and C II windows, and the array of times (since the start of the observations):

>>> sp = extract_irisL2data.load(raster_filename,window_info=lines[[0,7]])
>>> data_mg = sp.raster['Mg II k 2796'].data
>>> hd_mg = extract_irisL2data.only_header(raster_filename,extension=8)
>>> data_c = sp.raster['C II 1336'].data
>>> hd_c = extract_irisL2data.only_header(raster_filename,extension=1)

>>> from astropy.wcs import WCS
>>> m_to_nm = 1e9

>>> hd_c['CDELT3'] = 1.e-9  # Fix WCS for sit-and-stare
>>> wcs_c = WCS(hd_c)
>>> wave_c = wcs_c.all_pix2world(np.arange(data_c.shape[-1]),
>>>                              [0.], [0.], 0)[0] * m_to_nm
>>> hd_mg['CDELT3'] = 1.e-9  # Fix WCS for sit-and-stare
>>> wcs_mg = WCS(hd_mg)
>>> wave_mg = wcs_mg.all_pix2world(np.arange(data_mg.shape[-1]),
>>>                               [0.], [0.], 0)[0] * m_to_nm

>>> aux_hd=extract_irisL2data.only_header(raster_filename,extension=hd['NWIN']+1)
>>> aux=extract_irisL2data.only_data(raster_filename,extension=hd['NWIN']+1)
>>> time_diff=aux[:,aux_hd['TIME']]
>>> times = np.datetime64(hd['DATE_OBS']) + time_diff * np.timedelta64(1, 's')
The provided file is a raster IRIS Level 2 data file.

Extracting information from file iris_l2_20130902_163935_4000255147_raster_t000_r00000.fits...

Available data with size Y x X x Wavelength are stored in windows labeled as:

------------------------------------------------------------
Index in the output variable --- Window label --- Y x X x WL
------------------------------------------------------------
             0               C II 1336       417x1600x229
             1               Fe XII 1349     417x1600x53
             2               Cl I 1352       417x1600x36
             3               O I 1356        417x1600x53
             4               Si IV 1394      417x1600x150
             5               Si IV 1403      417x1600x262
             6               2786            417x1600x171
             7               Mg II k 2796    417x1600x555
             8               2830            417x1600x42
------------------------------------------------------------


The selected data are passed to the output variable
(e.g iris_data) as a dictionary with the following entries:

 - iris_data.raster['C II 1336']
 - iris_data.raster['Mg II k 2796']

For this dataset the spectral cadence is about 3 seconds. The Mg II k3 core is located around wavelength pixel 103. We can use this information to make a space-time image of the Mg II k3 wavelength:

>>> lev = im.histo_opt(data_mg[..., 103])
>>> plt.imshow(data_mg[..., 103], vmin=lev[0], vmax=lev[1], aspect='auto')
>>> plt.show()
-2.00e+02/5.75e+02 image min/image max
4.68e+00/2.32e+02 histo_opt min/histo_opt max
_images/tutorials_30_2.png

In the above, the axis correspond to pixel numbers. It is also possible to plot with the time and coordinate axis:

>>> # y scale, convert from degrees to arcsec
>>> y_arcsec = wcs_mg.all_pix2world([0.], np.arange(data_mg.shape[0]), [0.], 0)[1] * 3600.

>>> import matplotlib.dates as dates
>>> from datetime import datetime
>>> times_num = dates.date2num(times.astype(datetime))

>>> fig, ax = plt.subplots()
>>> ax.xaxis_date()
>>> date_format = dates.DateFormatter('%H:%M')
>>> ax.xaxis.set_major_formatter(date_format)
>>> ax.set_xlabel("Time")
>>> ax.set_ylabel("Slit position, solar Y (arcsec)")
>>> ax.imshow(data_mg[..., 103], vmin=0, vmax=200, aspect='auto',
>>>           extent=[times_num[0], times_num[-1], y_arcsec[0], y_arcsec[-1]])
>>> plt.show()
_images/tutorials_32_1.png

The middle section between 60”-75” is on the umbra of a sunspot, even though it is not obvious from this image. One can see very clearly the umbral oscillations, with a clear regular pattern of dark/bright streaks. Let us now load the 1400 slit-jaw and plot it for context:

>>> sji_filename = 'iris_l2_20130902_163935_4000255147_SJI_1400_t000.fits'
>>> sji = extract_irisL2data.load(sji_filename)
>>> sji_hd=extract_irisL2data.only_header(sji_filename)
>>> wcs = WCS(sji_hd)
>>> ax = plt.subplot(projection=wcs.dropaxis(-1))
>>> ax.coords[0].set_major_formatter('s.s')
>>> ax.coords[1].set_major_formatter('s.s')
>>> lev = im.histo_opt(sji.SJI['SJI_1400'].data[:,:,0],cutoff=5.e-2)
>>> ax.imshow(sji.SJI['SJI_1400'].data[:,:,0], vmin=lev[0], vmax=lev[1])
>>> plt.show()
The provided file is a SJI IRIS Level 2 data file containing SJI_1400 data.

Extracting information from file iris_l2_20130902_163935_4000255147_SJI_1400_t000.fits...

Available data with size Y x X x Image are stored in a windows labeled as:

------------------------------------------------------------
Index in the output variable --- Window label --- Y x X x Im
------------------------------------------------------------
             0               SJI_1400         417x388x400
------------------------------------------------------------




The SJI data are passed to the output variable
(e.g iris_data) as a dictionary with the following entries:

 - iris_data.SJI['SJI_1400']


-2.00e+02/1.25e+04 image min/image max
9.27e+00/1.50e+02 histo_opt min/histo_opt max

WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'seconds' -> 's'. [astropy.wcs.wcs]
_images/tutorials_34_3.png

The slit pixel 220 is a location on the sunspot’s umbra. We will use it to get some plots. For example, let’s plot the k3 intensity (spectral pixel 103 of data_mg) and the core of the brightest C II line (spectral pixel 90 of data_c) vs. time in minutes (showing first ~10 minutes only):

>>> plt.plot(times[:200], data_mg[220,:200, 103])
>>> plt.plot(times[:200], data_c[220,:200, 90] * 2)
>>> plt.ylim(0, 85)
>>> axis_lim = plt.axis()
>>> plt.show()
_images/tutorials_36_0.png

In the above we are multiplying the C II data by two to get the two lines closer. Image now you wanted to compare these oscillations with the intensity on the slit-jaw. How to do it? The slit-jaws are typically taken at a different cadence, so you will need to load the corresponding time array for the 1400 slit-jaw:

>>> sji_aux_hd = extract_irisL2data.only_header(sji_filename,extension=1)
>>> sji_aux = extract_irisL2data.only_data(sji_filename,extension=1)
>>> time_diff = sji_aux[:, sji_aux_hd['TIME']]
>>> times_sji = np.datetime64(sji_hd['DATE_OBS']) + time_diff * np.timedelta64(1, 's')

And now we can plot both spectral lines and slit-jaw for a pixel close to the slit at the same y position (index 220):

>>> plt.plot(times[:200], data_mg[220, :200, 103])
>>> plt.plot(times[:200], data_c[220, :200, 90] * 2)
>>> plt.plot(times_sji,sji.SJI['SJI_1400'].data[190, 220, :], 'y-')
>>> plt.axis(axis_lim)
>>> plt.show()
(735113.693821412, 735113.7013329861, 0.0, 85.0)
_images/tutorials_40_1.png

You are now ready to explore all the correlations, anti-correlations, and phase differences.