Monday, August 22, 2016

AIA Response Summary

This is my final blog post for GSoC. It contains a summary of the work that has been completed over the summer, and some points of how to extend this module further.

Project Goal: A SunPy AIA response module that includes temperature and wavelength response functions that will utilize ChiantiPy and the CHIANTI database.

Main PR Links:

Project Abstract:

Solar physics uses the CHIANTI atomic physics database to obtain properties about various elements and ionisation states. By using observed elemental abundances and ionization states, one can use CHIANTI to obtain synthetic spectra of solar plasma of various features which informs a response function used by the observational instruments themselves. This response function is vital to understanding observations.

The Atmospheric Imaging Assembly (AIA) is a multi-wavelength imager on the Solar Dynamics Observatory, specifically looking at the solar corona to understand magnetic processes.

This program uses SunPy and ChiantiPy  to calculate two response functions that infer plasma properties like temperature and density for the AIA:
  • Wavelength response functions: calculate the amount of flux per wavelength measured by AIA
  • Temperature response functions: calculate the sensitivity of light from the plasma per temperature measured by AIA

    Major Resources:

    Project based on Research Paper - Boerner et al. (2012)
    SunPy -  Sunpy User Guide

    To Run the Reponse Module:

    >>>  from sunpy.instr.aia import Response
    >>> path = 'path/to/ssw/aia/response/data.genx'
    >>> channel_list = [94,131,171,193,211,304,335] 
    >>> response = Response(path_to_genx_dir = path, version = 6, channel_list = channel_list)

    This starts the Response module, and then can be used to calculate and look at the wavelength response or the temperature response of the instrument.  

    To look at the calibration of one part of the instrument, one can use a simple loop for multiple channels, or like this for one channel:
    >>> primary_mirror = response.get_channel_data(94)['primary_mirror_reflectance']
    >>> secondary_mirror = response.get_channel_data(94)['secondary_mirror_reflectance']
    >>> total_reflectance = primary_mirror + secondary_mirror

    The transmittance , reflectance, quantum efficiency, and ccd contamination plotted versus the wavelengths (response.get_channel_data(94)['wavelength_range']) will return calibration plots. 

    To view an instrument response, the Effective area vs wavelength is one of the more important features as this defines the wavelength response for each channel. For that, just calculate the response, and then use the attribute for each channel. Using the values from a .genx file is recommended for best results, but a calculated effective area is available.

    >>> response.calculate_effective_area(94, use_genx_effarea = True)
    Similarly, the temperature response can be calculated and then called as an attribute. The function takes into account the contribution from each ion transition in small wavelength ranges. Due to using ChiantiPy, the population of a datatable is necessary and is expensive in terms of time. 

    >>> response. calculate_temperature_response(94)

    This function returns a dictionary giving the wavelength and response values that can be used to plot the function on a log-log scale. Alternate temperatures can be provided using the temperature_array keyword. 

    Final Thoughts:


    This project has incorporated many programs (SunPy, SSW, ChiantiPy) and many new ideas. I'm very proud to have constructed the response functions. With these implemented, it provides a way to do calibration using the instrument properties. However, many of the struggles in the project have been from understanding everything that was implemented in the SSW version of these functions, and as a result, many more features could still be added as explained in the next section. I'm very happy with the overall product and feel it is a good base (as I had hoped it would be starting out). A special thanks to my advisor’s Will Barnes and Drew Leonard. They have always been very supportive during this process and provided invaluable advice and direction.

    Next Steps:


    1. The program sunpy/instr/tests/ has housed all of my initial testing parameters that were not implemented in this project. This will need to be flushed out as additional features are implemented.
    2. SunPy could benefit from a general .genx reader. This was discussed in the developer meetings for this project. Currently, I have implemented a program with this module called This can be removed with the addition of a general .genx reader.
    3. SunPy could also benefit from having a defined file storage application. Currently, the temperature response in this code requires the user to calculate a datatable to store for future use. With a file storage application, this could be stored ahead of time, allowing for quick calculations as needed. 
    4. One additional feature that could be implemented is the ability to alter density in the temperature response. Currently, it is the constant 1e+9, but the ability to have this as an array of density = various pressures * temperature would allow for more generalization.
    5. The temperature response creates an ion datatable if one is not already created. This feature would greatly benefit from a way in ChiantiPy to search the transitions that happen in small wavelength ranges without having to create an ion object for each ion. The only way that I found was to use the ion.Intensity to see which transitions occur, but having to make an ion  object for each ion takes a long time.

    Thanks GSoC for this opportunity!


Tuesday, August 2, 2016

Temperature Response Functions


The AIA instrument on the Solar Dynamics Observatory doesn't have any spectroscopic abilities so it has no way to calibrate the data directly. My projects wavelength response and temperature response function are necessary to infer plasma properties like temperature and density using this instruement. Essentially, the temperature response function calculates the sensitivity of light from the plasma at specific wavelengths and temperatures as measured by the AIA. 

These plots are reproductions of those found in Boerner et al. (2012) . Figure 11 in this paper utilized initial calibrations of AIA.  I use last posts wavelength response function (that utilize instrument information as found in Version 6 of SSW data files) to obtain this new function:

Temperature Response Functions

This displays the temperature response for all six extreme ultra violet channels (EUV) observed with the AIA instrument vs Temperature on a log 10 scale. It is calculated by multiplying the wavelength response function times an ion contribution function, and then iterating through all of the wavelengths in the desired region.

The wavelength response function that I previously implemented tell the efficiency of the optics around each instruments channels wavelength center. The new and difficult part of this function was implementing the contribution function of each ion. This is where I've utilyzed ChiantiPy. This python package accesses CHIANTI, an 'Atomic Database for Spectroscopic Diagnostics of Astrophyiscal Plasmas'. As well as having a large range of atomic elements, this code analyses each ion to get attributes such as ionization, intensity of various lines, and emission measurements. The contribution function in ChiantiPy uses these  parameters and is strongly peaked in temperature. With it, we can gauge the sensitivity of light from the plasma seen by the AIA instrument as a function of wavelength and temperature.

Temperature Response Per Channel VS Temperature  

 Channels: 94, 131, 171, 193, 211, 304, 335
Temperature is on a log base 10 scale

The values for these response functions are still a bit higher than anticipated. At first, I tried incorporating the plate scale of the instrument to see if this helped, but unless I also multiply by an assumed density, the units for this do not match to the expected result. 

Finishing up Goals:

  1.  I edited ChiantiPy so that the contribution function would not plot. This allowed me to obtain the values for multiple ions at a time. Unfortunately, this means that I am needing this to be added to ChiantiPy on github if this code is going to be implemented into Sunpy on github. The issue I ran into with this is that I had been editing a lot of my version of ChiantiPy while I was testing out different methods. I plan reinstalling the github version, inserting my changes in just the contribution function, and submitting a pull request asap!
  2. One change that was recommended in the SunPy developer meeting last week is to make a dataframe containing the ion contributions. This will speed up the process if it can be read each time instead of calculating them each time. I think this is a great idea and am in the midst of implementing it.
  3. I'd like to make this code as general as possible. This means adding more keywords so the user can specify their own temperature, density, ions, etc. 
  4. Documentation! This is something I've taken time to read up on more this week because I recognize how wonderful convention can be when making code more readable. It is a process that I am still working through. 

Sunday, July 10, 2016

 AIA Instrument Wavelength Response Working!


The past few weeks started with success in getting code that worked up to a gain calculation, and then was followed with working through bugs while I implemented changes to restructure the code and calculating the gain.



  Newly Implemented Code Features:

  • .genx properties now saved to Astropy QTable structure to utilize storage of data type and units.
  •  Effective Area function now can load for the whole wavelength range or for just the channel specified
  • Instrument Wavelength response outputs values close to what we expect!

 Table Structure:


Previously, I was storing the information gathered from .genx files in a dataframe. This worked, but we wanted a way to store units and data types with the actual values as well. This is possible with a new and exciting thing called Astropy Table structures. (Okay, new and exciting to me.) But, figuring out the best method for filling these was difficult. Here is what worked:

First, I made a table:

> table = QTable(names = array_of_property_names , dtype = array_of_data_types)

Initially, I was trying to fill a whole column of channel data (the inverse axis of what is now implemented). This worked fine, but there isn't a way to save different data type values in one column. So, I filled the table row by row with an array for each channel where each element was a value corresponding to the property in the column name.

> table.add_row(array_of_channel_properties)

I also added in a Column after the table was filled. This effectively inserted the channel names for each row to be used as an index.

> C = Table.Column(data=channel_list, name='channel', dtype=int)
> table.add_column(C, index = 0)
> table.add_index(['channel'])

One issue that came up was that any np.ndarray that was stored in the array would not be loaded into the row as I was requesting.  I found a  way around this was to list the data type for these as 'object'. Then, when they are read in to the program as one object, I simply iterated through them and restored the elements as a np.array. 

 Gain Calculations:

The last post displayed instrument calibration plots. Of those, the effective area function is the most important because it gives us an idea as to how the instrument is responding at wavelengths where we see strong emission spectra from solar features. The values obtained at these wavelengths are multiplied by the wavelength dependent gain of the ccd system to get a wavelength response function, as defined by Boerner et al. (2012). Sounds so simple, but it has taken a while to nail down.

Initially, I used the entire wavelength range for each channel to calculate the Effective Area and then plotted them around specific wavelengths. When I went to use these values for calculating gain, they turned out to be tiny (eg 10^ -12)! After some research, I divided this by the plate scale to get a reasonable number, but it still was off. 

As it turns out, I had found a hidden bug. Any arrays loaded in from ssw packages contain properties of the AIA instrument, and cover some wavelength range. I found that the index of these arrays at the value that corresponded to a desired wavelength was not the same as the wavelength number (which was a subtle but quick assumption on my part). After changing the way the arrays were indexed (to correspond to the desired wavelength),  I obtained nice values from the Effective Area function and worked on calculating the gain of the system. 

Here is a chart listing the wavelength response values for each channel:

Wavelength Response

SSW imported Eff-area Eff-area function Eff-area function

Boerner Table 2 SSW imported gain Table 2 Gain Table 6 ccd gain
channel Reference calculated response calculated response calculated response
94 0.664 0.572 0.617 0.524
131 1.785 1.87 1.933 1.714
171 3.365 3.257 3.389 2.987
193 1.217 1.79 1.906 1.642
211 1.14 1.101 1.184 1.01
304 0.041 0.038 0.0413 0.035
335 0.027 0.0255 0.026 0.023

Column 2 lists the values we expect from this function.
Column 3 is the calculation using only SSW imported values. 
Column 4 is how I tested the effective area function - by using the gain provided by the reference table 2 in Boerner 2012. 
Column 5 are the values obtained using my own gain calculation using the ccd instrument gain.

 Voila! With these values, the wavelength response function is complete! 


Now, my mentors and myself will need to decide which is the best to use moving forward, and as a preview, the next blog will be a detailed report on the next function..... Emissivity!


Friday, June 24, 2016

And Then There Were Plots!


These plots are reproduced from Boerner et al. (2012), which utilized initial calibrations of AIAI have now developed enough of the response module to analyse updated instrumental properties to get my own calibrations using AIA instrument outfiles from SSW.

Filter Calibration

The transmittance is a ratio of the signal received through the instrument filters.  The values match the lowest lines modelled in Figure 2 of the paper.

 Mirror Calibration


This shows the reflectance of the mirrors per wavelength, and it is exactly like Figure 5 in the paper. Huzzah!

CCD Calibration  

Updated: Here the quantum efficiency of the CCD is shown per wavelength, and now it looks like Figure 6 after loading in the UV channels.

Effective Area Functions

Each effective area function displays the efficiency of the optics around the center wavelength of each instrument channel. It is calculated using the instrument calibration properties shown in the other plots.
Updated: Since re-factoring my code, I now get the shapes expected from the paper. These show the calculated effective area in blue and the effective area loaded in from the .genx file dashed in red.


Monday, June 20, 2016

GSoC AIA Response Functions


Project Details:


First:  A bit of a delay...

The past few weeks have been quite busy, yet rewarding. The culmination of finals week, graduation, family visiting the West coast for the first time, and finishing my senior research by doing some final remote observations have been taking up more time than anticipated. I'm sad to say that all of this stunted my progress in the GSoC project.

Then:  Back to Work!

So, in hopes to catch up, this last week I've worked with a technician on UW campus to get SSW working. Several days later, we finally just root installed ftp to get it to download. Now, I have access to an idl version of what I'm trying to redo! This has helped already to see what each piece of the code outputs.


First, I read in the idl .genx files from ssw with python.  My mentor Will sent me documentation on how to do this with I need some properties in these files to make the instrument response function. I was warned that this would make numpy.recarray, and have found this to be a challenging object that works with fields...? It took me a few days to figure out just how nested these objects can be!

My first function I completed was utilizing the information from the .genx files and putting them into pandas dataframes. I used each channel as a key in my next function, which is to load each channel individually.

The past few days I have been working on implementing the response function with these properties. So far, I'm having issues loading in the parameters to calculate an effective area. Mostly, it's been getting past multiplying floats and arrays. My first attempt as this has been to use numpy.ones to make everything the same size, but then the division throws errors. I am still working on it!


  1. I made a pull request in Sunpy to make it easier for the developing team to see my work. My goal for this next week especially is to make the habit of committing code to github as I go, rather than in a big bundle like I did for this past commit. It will help to show my gradual progress, and it will help my mentors to guide the work as well.  
  2. I plan to do another post this week to make up for missing one last week. It will include the pull requests that I have been working through for the past week. 
  3.  Get some plots of the instrument response functions asap for midterm evaluation!!
  4. Be ready and on time for the developers meeting on Wednesday, 9 am my time!!!


Sunday, May 29, 2016

Guidance provided by Sunpy developers guide

  • Found while reading documentation

  • Starting to be utilized while organizing the response modules

I'm really starting to know and love the SunPy developers guide. One guideline, also recommended by Astropy:docs, is setting up a python virtual environment. This makes things easier by allowing for specific versions of packages to be used per environment.
Here's how one would start an environment using conda, which is included in the installation of Anaconda:
>>> conda create name_env python=2.7
Then, any imports for the different requirements match that python, so switching between using python 2.7 (and all of the dependencies that depend on it) to python 3.4 (with new dependencies) is as easy as pie. These commands below are how you quickly go between one environment and another.
>>> source activate name_env
>>> source deactivate
>>> source activate name2_env
Lastly, this first command makes it easy to see what's being loaded in each environment, and the second one is how to install more packages:
>>> conda list
>>> conda install name_of_package

Pretty sweet. This has already helped me in testing ChiantiPy and Sunpy, and by helping me solve a problem I ran into in the first weeks with their respective python versions. I looked into popular/new pythons and found 2.7 and 3.4 are widely used. When I researched more, I found:

- The sunpy documentation states it works best with python 2.7.
- The ChiantiPy documentation has been tested with python 3.3 specifically.

So, I chose to set up an environment with python 2.7 because I've found this loads both modules. Then, I like to run sunpy.self_test() afterwards. 

       I also found similar git guidelines while reading the documentation pages provided by Astropy. I've utilized the section on 'How to make code contribution' quite a bit. Last post I mentioned setting up a branch on github, and the guides I've mentioned helped with that process. Mostly, I've used git to save my own repositories, but working with the development version of a code is much more extensive.

       Sunpy has asked that each week the participants of GSoC look into one pull request and check for validity. This is a new process for me, but I welcome the practice of working in the community-developing environment. While scanning the pull requests on Sunpy, I found one that reminded me of the astropy documentation, and sparked this blog post.

Sunpy Pull Request testing: Updates to #1656

  • First, I recognized .md from my work in preparation for this project where I resolved an issue that led to the development of in Sunpy
  • For this pull request, I read the full changes to the file with and without the notes provided.
  • Few things I liked:
    • The instant connection to people via recommending the IRC channel. That's awesome. 
    • That an example of an issue is provided. Great recommendation.
  • Inputs:
    • Wafels comments covered quite well any thoughts that I had.

Overall, I feel this is a well developed file that is ready to be implemented. I have not tested this yet, as the build process is unfamiliar.  Looking at the details in other developers builds has allowed me to follow their process. It seems to mirror what I've described in this blog of setting up an environment to test the file, and making sure no errors come out. I'm working through this part, and feel this is something to ask about in the SunPy IRC.





Wednesday, May 25, 2016

Community Bonding Period Post

Understanding what goes into AIA Response Functions




Many of the references sent by my advisor’s for this project show interesting applications of temperature response functions that I've been scanning for details, and the science is great. One reference is the Initial Calibrations of AIA on SDO. This paper explains that the instruments on AIA do not give spectroscopic information, so there is no way to calibrate the data directly. They give the instrument-response function below, that I know I'll be utilizing:
"instrument-response function [R(λ,t)=Aeff(λ, t)G(λ)]"

Learning the Coding Environments:

During the community bonding period, I have taken advantage of a few resources to help me get used to the SunPy and ChiantiPy coding environments. One is the ChiantiPy tutorial. I was advised that the contribution functions are vital to this project, so I'll offer a summary of that here. This is just a small piece of what will go into the modules.

The goal is to analyse emission spectra. Chianti.core has the ability to analyse the emissivities of the three strongest lines around 200 Angstrom using a contribution function. After importing this module, I defined the temperature and density of iron 14 to use for testing and applied the contribution function.

>fe14 = ch.ion('fe_14', temperature=t, eDensity=1.e+9, em=1.e+27)
>fe14.gofnt(wvlRange=[210., 220.],top=3)

I really like how this is automatically stored in at dictionary, with keys = [‘gofnt’, ‘temperature’, ‘density’]. This will be quite beneficial later when I am running this on a number of lines

I also tested this for a different line of Iron to show the peaks in temperature (I chose 177.240) at a specific density :

> fe10 = ch.ion('fe_10', temperature=t, eDensity=1.e+9, em=1.e+27)
> fe10.gofnt(wvlRange=[170., 180.], top=3)


On a side note, I successfully made a new branch in github!

REVISION:  Thanks to the recommendations of my advisor Will Barnes, I double checked the versions of Chianti that were being loaded. As it turns out, the version I had downloaded from github was not the one being loaded. This was causing discontinuous spikes in the contribution functions that were not suppose to be there. The plots have been updated to have the right continuous shape after loading the correct ChiantiPy, and I switched to a temp of 174 to 177 for Iron 10.



For the first weeks of coding, here are my goals: 
  • My first task is to implement a new response function module into sunpy/sunpy/instr. This will be on my version of sunpy in github.
  • In this module, I am filling in the different classes and definitions that will need to be developed, and am implementing some documentation for each one. 
  • As part of GSoC, I found out recently that I need to review a pull request in Sunpy each week as well. This will be good practice for me, and I'll try to incorporate this into the blog. 


Sunday, April 24, 2016

GSoC Proposal Accepted!

Approved Proposal Shown here!

I'm so very excited that my GSoC project working with Sunpy has been accepted! 

In the next month, I'll be becoming more acquainted with both Sunpy and CHIANTI. Then, coding will ensue on this project. I'm also excited to be working with mentors Drew Leonard and Will Barnes.

Thanks for this opportunity Google!


Monday, April 4, 2016

Getting into SunPy

Part 2: Pull Request: Spectrogram log y axis #291  


I have found data that can be used with Sunpy's spectrogram class! A solar radio spectrometer called CALLISTO is used internationally for networking solar radio data (with direct access!). Technically, I am using the CallistoSpectrogram class in sunpy.spectra to open data that I'll use to test the Spectrogram class. Using this data, I am confident that the format loaded into the class I'd like to test is correct.

>>> import spectrogram
>>> import matplotlib.pyplot as plt
>>> from sunpy.spectra.sources.callisto import CallistoSpectrogram
>>> spectra = CallistoSpectrogram.from_file('')
>>> gram = spectrogram.Spectrogram(, spectra.time_axis, spectra.freq_axis, spectra.start, spectra.end)
>>> gram.plot()
Woot, data loaded successfully!

Next, the pull request asks if: instead of the frequency being on a linear scale, is it possible to make it a log scale? Looking at the spectra.spectrogram class in sunpy, I see that the linear frequency scale took some work to implement. Simply implementing a logy keyword would be confusing with the linear keyword already in place. So, I am trying to implement the axes.set_yscale('log') after the linear scale is established. Perhaps this could be changed just before plotting for an easy addition to the class. Otherwise, I plan to implement a log frequency axis in a similar way to the linear frequency axis.

To be continued.


reference to helpful sunpy help page

Saturday, March 26, 2016

Getting into SunPy 

Attempting Pull Request: Spectrogram log y axis #291 


I'm love spectra. So, when the GSoC process required that I dive into issues in the open source Python library SunPy, I looked for something regarding spectra. I'm still getting used to GitHub, though, so attempting to change the code of such a large and utilized library is kind of intimidating. Then, I found the pull request linked above using a spectrogram... Great! I have low resolution data that I can use to test this!

As it turns out, this class is still being developed. I also realized that a spectrogram is slightly different than the spectra I've been working with so far.

A spectrogram show a spectrum of frequencies as they vary with time. It is used to visualize sounds plotted with time on the x axis and frequency on the y axis. This is different, but closely related to a spectra or spectrum of light that I am accustomed to, where frequencies (or wavelengths) are plotted on the x axis and the amount of of light (or flux) is plotted on the y axis. This means that I probably shouldn't try to test it with  my visual range spectra. Instead, this is most likely working with radio frequencies.

This makes sense for solar data. Radio bursts from the corona are observed, and need to be analysed.

Getting Started:

I found the code needing to be tested. On first glance, I noticed that the spectrogram class has a plot definition, and thought that a simple 'logy' keyword and an if statement added in the plot definition would suffice, like this:
  if logy:

I need to test this though, and started with loading in the module and inserting a basic 2D array of random data:

>>> import spectrogram
>>> import numpy as np
>>> data = np.ndarray((10,10), buffer = np.array(np.random.random(100)))
>>> freq_axis = np.linespace(30, 300, 10)
>>> time_axis = np.array(range(0, 10, 1))

Then, I set up some arbitrary start and end times (10 seconds) to suffice the keywords:
>>> import datetime
>>> start = datetime.datetime(2016, 2, 2, 12, 30, 30)
>>> end = datetime.datetime(2016, 2, 2, 12, 30, 40)

Put it all together:
>>> gram = spectrogram.Spectrogram(data, time_axis, freq_axis, start, end)

And, just like that I have a spectrogram object. Next, to test the plot function:
>>> gram.plot()

Colorful! Now, let's try this with my addition if statement:
>>> gram.plot(logy=True)

This flopped with many, many errors. and it is not showing a a frequency range. So, I think my data may not be in the right format.

I tried it again, this time with making a sine wave with 10 being the sampling rate.
>>> sine = np.sin(2 * np.pi * (3*10**4) * (np.arange(100)) / 10)
>>> data2 = np.ndarray((10,10), buffer=sine)
>>> gram2 = spectrogram.Spectrogram(data2, time_axis, freq_axis, start, end)

It produces:
Looking less like noise. But, still not making a frequency axis. I'm going to look for data to use so that I can reproduce the use of this module, and then add my changes to test the log y axis change.


Thursday, March 24, 2016

Implementing AIA response functions in SunPy 


Project Abstract:

Solar physics uses the CHIANTI atomic physics database to obtain properties about various elements and ionisation states. By using observed elemental abundances and ionization states, one can use CHIANTI to obtain synthetic spectra of solar plasma of various features which informs a response function used by the observational instruments themselves. This response function is vital to understanding observations.

The Atmospheric Imaging Assembly (AIA) is a multi-wavelength imager on the Solar Dynamics Observatory, specifically looking at the solar corona to understand magnetic processes.

This project aims to use SunPy to infer plasma properties like temperature and density by
developing the routines necessary to calculate two response functions for the AIA using python and ChiantiPy:
  • Wavelength response functions: calculate the amount of flux per wavelength measured by AIA
  • Temperature response functions: calculate the sensitivity of light from the plasma per temperature measured by AIA

Why am I interested?

My goal is to explore the connection between what we observe on the sun (heliophysics) and how it relates to what we see in other stars (astrophysics). My current research has been looking at visual wavelength range spectra for stellar objects such as M Dwarfs, K Giants, and eclipsing binaries. I would like to use this summer to expand my knowledge on solar observations, and the modules and vernacular that come with solar research.


I'm excited to learn about ChiantiPy,  the python interface for astrophysical spectroscopy using the atomic database  CHIANTI. I've observed spectra before, but this project is all about understanding observations by the AIA through python code. This project also wants to get away from using SolarSoftWare (SSW), an idl libray, to calculate these response data structures.  


 In reading Boerner et al 2012, I'm learning that the response functions need to use instrument calibration to better the results. The main idea for this project is to analyse the AIA  instrument responses to output a wavelength response, and then use a spectral model to obtain a temperature response.  Now back to work on my proposal! 



Source of SunPy project idea