SunPy: Python for Solar Physics

A Quick Overview

Loading and plotting images in SunPy
Loading and plotting timeseries in SunPy

First things first again, get this data downloading:



In [ ]:

    
!wget http://sunpy.cadair.com/data/aia_lev1_171a_2014_01_01t00_02_23_34z_image_lev1_fits.fits



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

1. Solar Images with SunPy

SunPy has a Map type that supports 2D images, it makes it simple to read data in from any filetype supported in sunpy.io which is currently FITS, JPEG2000 and ANA files. You can also create maps from any (data, metadata) pair.

SunPy contains a few test data files, under the sunpy namespace:

sunpy.AIA_171_IMAGE
sunpy.CALLISTO_IMAGE
sunpy.EIT_195_IMAGE
sunpy.RHESSI_EVENT_LIST
sunpy.RHESSI_IMAGE

These are all just links to FITS files contained in the SunPy distribution.



In [2]:

    
import sunpy
print sunpy.AIA_171_IMAGE









    



/home/stuart/GitHub/sunpy/sunpy/data/sample/AIA20110319_105400_0171.fits



In [3]:

    
import sunpy.map

aiamap = sunpy.map.Map(sunpy.AIA_171_IMAGE)
aiamap









    Out[3]:





SunPy AIAMap
---------
Observatory:	 SDO
Instrument:	 AIA_3
Detector:	 AIA
Measurement:	 171
Obs Date:	 2011-03-19T10:54:00.34
dt:		 1.999601
Dimension:	 [1024, 1024]
[dx, dy] =	 [2.400000, 2.400000]

array([[ 0.3125, -0.0625, -0.125 , ...,  0.625 , -0.625 ,  0.    ],
       [ 1.    ,  0.1875, -0.8125, ...,  0.625 , -0.625 ,  0.    ],
       [-1.1875,  0.375 , -0.5   , ..., -0.125 , -0.625 , -1.1875],
       ..., 
       [-0.625 ,  0.0625, -0.3125, ...,  0.125 ,  0.125 ,  0.125 ],
       [ 0.5625,  0.0625,  0.5625, ..., -0.0625, -0.0625,  0.    ],
       [ 0.5   , -0.125 ,  0.4375, ...,  0.6875,  0.6875,  0.6875]])

Maps contain both the image data and the metadata associated with the image, this metadata currently does not deviate much from the standard FITS WCS keywords, but presented in a instrument-independant manner.

Plotting Maps

1. The quick way

Most SunPy datatypes have a 'quick view' method called peek() which returns a Figure object and presents an overview of the data.



In [4]:

    
fig = aiamap.peek()









    



/usr/lib/python2.7/site-packages/matplotlib/figure.py:371: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "



In [5]:

    
plt.show()

2. The powerful way

All SunPy datatypes have a plot() method that integrates with matplotlib in a way that allows you to customise your plot.



In [6]:

    
#Create a Figure and an axes
fig, ax = plt.subplots(1)
#Plot the aiamap, aiamap.plot() behaves like plt.imshow()
im = aiamap.plot(axes=ax)
#Add a colorbar
plt.colorbar(im)
#Overlay a grid
aiamap.draw_grid(axes=ax)









    



/home/stuart/GitHub/sunpy/sunpy/wcs/wcs.py:180: RuntimeWarning: invalid value encountered in sqrt
  distance = q - np.sqrt(distance)






    Out[6]:





<matplotlib.axes.AxesSubplot at 0x47d2950>

Using Maps

Commonly used meta data and data manipulation methods are part of the map class, below I highlight a couple. To see the full features of Map look at the GenericMap documentation.

First we read in a Level 1 AIA image:



In [7]:

    
aiamap2 = sunpy.map.Map('aia_lev1_171a_2014_01_01t00_02_23_34z_image_lev1_fits.fits')



In [8]:

    
fig = aiamap2.peek()

The properties of the map provide a standard way of reading meta data. For instance the pixel scale can be accessed with the .scale attribute:



In [9]:

    
aiamap2.scale









    Out[9]:





{'x': 0.599489, 'y': 0.599489}

It is also possible to rotate and rescale an image:



In [10]:

    
aiamapR = aiamap2.rotate(angle=np.pi)



In [11]:

    
fig = aiamapR.peek() # Oh noes the Sun is upside down.

Warning, this method will not change the meta data:



In [12]:

    
aiamapR.rotation_angle









    Out[12]:





{'x': 0.0, 'y': 0.019339}

Exercise 3:

Using map.submap() cut out an active region of your choice, over plot a lon/lat grid at a sutible resolution depending on your box size and rescale the image limits when you plot the output with a colour bar.

Solution:



In [13]:

    
aiasub = aiamap2.submap([800,1200],[-400,100])



In [14]:

    
fig, ax = plt.subplots(1, figsize=(7,7))
im = aiasub.plot(axes=ax, vmin=0, vmax=3300)
ax = aiasub.draw_grid(grid_spacing=10)
cbar = plt.colorbar()

2. Timeseries and Lightcurves

SunPy supports many different sources of time series data via the Lightcurve data type. Currently, the LightCurve class is compatible with the following data sources: the GOES X-ray Sensor (XRS), PROBA2/LYRA, and SDO/EVE

The first Lightcurve we are going to create is for 1st January 2014 from the GOES XRS instrument.



In [15]:

    
import sunpy.lightcurve



In [16]:

    
goeslc = sunpy.lightcurve.GOESLightCurve.create('2014/01/01T00:00:00','2014/01/01T23:59:59')









    



/home/stuart/GitHub/sunpy/sunpy/lightcurve/lightcurve.py:223: RuntimeWarning: Using existing file rather than downloading, use overwrite=True to override.
  warnings.warn("Using existing file rather than downloading, use overwrite=True to override.", RuntimeWarning)



In [17]:

    
fig = goeslc.peek()

Currently there is not a huge amount of functionality actually in the Lightcurve object, the power of Lightcurve comes from the pandas DataFrame stored in the Lightcurve.



In [ ]:

    
goeslc.data

If we want to just plot the $0.5$ - $4.0$ A or the 'A_FLUX' column we can do this:



In [18]:

    
ax = goeslc.data['A_FLUX'].plot()

To customise this plot you can specify lot's of arguments to the pandas DataFrame .plot() method, which is described here.

Exercise 4:

Create a new LightCurve object for the same day in January, but this time using the SDO/EVE instrument. Use EVELightcurve. Pay special attention to the parameters accepted by the EVE .create() method.

Once you have the lightcurve, plot it and explore the data contained in the pandas DataFrame. Then plot the $17.1$ channel on it's own.

Solution:



In [19]:

    
evelc = sunpy.lightcurve.EVELightCurve.create('2014/01/01')



In [20]:

    
fig = evelc.peek()



In [21]:

    
ax = evelc.data['17.1ESP'].plot()

We can see that the missing data in this file hasn't been properly removed. (This is probably a bug in EVELightCurve). We can check the header:



In [22]:

    
evelc.header[8]









    



/home/stuart/GitHub/sunpy/sunpy/lightcurve/lightcurve.py:85: Warning: lightcurve.header has been renamed to lightcurve.meta
for compatability with map, please use meta instead
  for compatability with map, please use meta instead""", Warning)






    Out[22]:





'; Missing data: -1.00e+00\n'

We can now use pandas to convert the missing data to NaN:



In [23]:

    
evelc.data[evelc.data == -1] = np.nan



In [ ]:

    
evelc.data



In [24]:

    
evelc.data['17.1ESP'].plot()









    Out[24]:





<matplotlib.axes.AxesSubplot at 0x8df9750>

Which is much nicer!!

Exercise 5:

Using matplotlib plot both the $17.1$nm EUV flux and the GOES X-ray flux on the same graph.

Hints/Bonuses:

The timeseries data in a Pandas DataFrame is held in the .index attribute, use this as the x axis for your plots.
plt.twinx will give you two independant scales using the same x axis.
Plot a legend.
Use fig.autofmt_xdate() to magically format temporal x axis better.
Use ax.semilogy() to set the GOES axis to a log y scale.

Solution:



In [25]:

    
#Create the start Figure and Axes
fig, ax = plt.subplots(1)

#Plot the EVE 17.1 line
line1 = ax.plot(evelc.data.index, evelc.data['17.1ESP'], color='blue')

#Split the x axis between two y axes, returns a new Axes instance.
ax2 = ax.twinx()
#Plot the GOES flux
line2 = ax2.plot(goeslc.data.index, goeslc.data['A_FLUX'], color='red')
#Set the y scaling to logarithmic
ax2.semilogy()

#Autoformat the date
fig.autofmt_xdate()

#Manually set the legend because it will not auto-detect over two axes.
plt.legend((line1[0],line2[0]), ("EVE $17.1$nm Flux", "GOES $0.5$ - $4.0$ $\AA$ Flux"), loc=2)









    Out[25]:





<matplotlib.legend.Legend at 0x4703550>