#!/usr/bin/env python
# coding: utf-8

# # TESS 2015 - Introduction to Solar Data Analysis in Python

# <img src="https://raw.github.com/sunpy/sunpy-logo/master/generated/sunpy_logo_compact_192x239.png">
# # SunPy!
# ## Author: Steven Christe
# Email: steven.christe@nasa.gov

# Through tutorials and presentations we will demonstrate how free and open-source Python and SunPy library can be used to analyze solar data. Depending on interest, dinner may be ordered after the main presentation (roughly an hour) and a hands-on help session will take place for the remainder of the evening.  Installing the following software is recommended but not required: Anaconda Python and SunPy.

# Schedule: 19:00 to 22:00
# + 18:00-18:15 pm: Organization (chatroom, organize dinner)?
# + 18:15-19:00 pm: Intro to SunPy (presentation)
# + 19:00-19:15 pm: Break
# + 19:15-20:15: Dinner & Installation Workshop
# + 20:15-22:00: SunPy Workshop

# ## What is SunPy?
# ### A community-developed, free and open-source solar data analysis environment for Python.

# + website: [http://www.sunpy.org](http://www.sunpy.org)<br>
# + documentation: [http://docs.sunpy.org](http://docs.sunpy.org)<br>
# + code (version control!): [https://github.com/sunpy/sunpy](https://github.com/sunpy/sunpy)
# + Mailing list: https://groups.google.com/forum/#!forum/sunpy**
# + IRC: #sunpy on freenode.net [web client](https://kiwiirc.com/client/irc.freenode.net/#SunPy)

# SunPy is built upon foundational libraries which enable scientific computing in Python which includes
# 
# + [NumPy](http://numpy.org)
# + [SciPy](http://scipy.org)
# + [matplotlib](http://matplotlib.org)
# + [AstroPy](http://astropy.org)

#   <p>Supported observations
#       <ul>Images
#           <li>SDO AIA and HMI</li>
#           <li>SOHO EIT, LASCO, MDI</li>
#           <li>STEREO EUVI and COR</li>
#           <li>TRACE</li>
#           <li>Yohkoh SXT</li>
#           <li>RHESSI mapcubes (beta)</li>
#           <li>PROBA2 SWAP</li>
#           <li>IRIS Slit-Jaw (beta)</li>
#         </ul>
#     <ul>Time Series
#             <li>GOES XRS</li>
#           <li>PROBA2 LYRA</li>
#           <li>Fermi GBM</li>
#           <li>SDO EVE</li>
#         <li>RHESSI Summary Lightcurves</li>
#           <li>Nobeyama Radioheliograph LightCurve</li>
#         <li>NOAA Solar Cycle monthly indices</li>
#         <li>NOAA Solar Cycle Prediction</li>
#     </ul>
#     <ul>Spectra
#         <li>Callisto</li>
#         <li>STEREO SWAVES</li>
#     </ul>
#       
# <p>Supported data retrieval services 
#     <ul>
#         <li>Virtual Solar Observatory (VSO)</li>
#         <li>JSOC</li>
#         <li>Heliophysics Events Knowledgebase (HEK)</li>
#         <li>Helio</li>
#         <li>Helioviewer</li>
#     </ul>
# </p>
#  
# <p>Supported file formats include
# <ul>
#     <li>FITS (read/write)</li>
#     <li>Comma-separated files, text files (read/write)</li>
#     <li>ANA (read/write)</li>
#     <li>JPG2 (read/write)</li>
# </ul>
# </p>

# # What is this webby notebooky magic?!
# ## ipython notebook (now known as jupyter)
# similar to matlab, mathematica, maple

# Some setup...

# In[ ]:


import numpy as np
import matplotlib.pyplot as plt
import matplotlib
get_ipython().run_line_magic('matplotlib', 'inline')
matplotlib.rc('savefig', dpi=120)
import warnings
warnings.simplefilter("ignore", Warning)
from matplotlib import dates


# In[ ]:


import sunpy


# In[ ]:


sunpy.system_info()


# SunPy version (stable) 0.5

# # Let's study a flare!

# ## Searching for events in the HEK

# In[ ]:


from sunpy.net import hek
client = hek.HEKClient()
tstart, tend = '2014/01/01 00:00:00', '2014/01/02 00:00:00'
result = client.query(hek.attrs.Time(tstart, tend), 
                      hek.attrs.EventType('FL'), 
                      hek.attrs.FRM.Name=='SSW Latest Events')
len(result)


# In[ ]:


result[0]


# In[ ]:


for res in result:
    print(res.get('fl_goescls'))


# In[ ]:


result = client.query(hek.attrs.Time(tstart, tend), 
                      hek.attrs.EventType('FL'), 
                      hek.attrs.FRM.Name=='SSW Latest Events', 
                      hek.attrs.FL.GOESCls>'M')


# In[ ]:


len(result)


# In[ ]:


result


# We can find out when this event occured

# In[ ]:


result[0].get('event_peaktime')


# and where it occurred

# In[ ]:


result[0].get('hpc_coord')


# # Lightcurves!
# ## Let's look at the GOES curve for this event.
# First some time manipulation!

# In[ ]:


from sunpy.time import TimeRange, parse_time
from datetime import timedelta


# In[ ]:


tmax = parse_time(result[0].get('event_peaktime'))


# In[ ]:


tmax


# In[ ]:


tr = TimeRange(tmax - timedelta(minutes=30), tmax + timedelta(minutes=30))


# In[ ]:


tr


# In[ ]:


from sunpy.lightcurve import GOESLightCurve
goes = GOESLightCurve.create(tr)


# In[ ]:


goes.peek()


# The data is stored in a standard place

# In[ ]:


goes.data


# This is a pandas dataframe! Provides lots of additional functionality. For example

# In[ ]:


print('The max flux is {flux:2.5f} at {time}'.format(flux=goes.data['xrsb'].max(), time=goes.data['xrsb'].idxmax()))


# Compares well to the official max from the HEK

# In[ ]:


str(tmax)


# In[ ]:


goes.peek()
plt.axhline(goes.data['xrsb'].max())
plt.axvline(goes.data['xrsb'].idxmax())


# Meta data is also stored in a standard place

# In[ ]:


goes.meta


# This is a dictionary like the hek results so...

# In[ ]:


goes.meta.get('COMMENT')


# In[ ]:


goes.data.resample('10s', how='mean')


# # Solar Images in 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.

# Let's download an AIA image of this flare from the vso

# In[ ]:


from sunpy.net import vso
client=vso.VSOClient()


# Then do a search.

# In[ ]:


recs = client.query(vso.attrs.Time(tr), vso.attrs.Instrument('AIA'))


# It works, now lets see how many results we have!

# In[ ]:


recs.num_records()


# That's way too many!

# So that is every image that SDO/AIA had on that day. Let us reduce that amount.
# 
# To do this, we will limit the time search for a specify a wavelength.

# In[ ]:


recs = client.query(vso.attrs.Time('2014/01/01 18:52:08', '2014/01/01 18:52:15'), 
                    vso.attrs.Instrument('AIA'),
                    vso.attrs.Wave(171,171))


# In[ ]:


recs.num_records()


# In[ ]:


recs.show()


# Let's also grab another wavelength for later.

# In[ ]:


recs = client.query(vso.attrs.Time('2014/01/01 18:52:08', '2014/01/01 18:52:15'), 
                    vso.attrs.Instrument('AIA'),
                    vso.attrs.Wave(94,171))


# In[ ]:


recs.num_records()


# Let's download this data!

# In[ ]:


f = client.get(recs, methods = ('URL-FILE_Rice')).wait()


# In[ ]:


f


# For SunPy the top level name-space is kept clean. Importing SunPy does not give you access to much. You need to import specific names. SciPy is the same. 
# 
# So, the place to start here is with the core SunPy data object. It is called Map. 

# In[ ]:


from sunpy.map import Map


# In[ ]:


aia = Map(f[1])


# In[ ]:


aia


# 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-independent manner.

# In[ ]:


aia.peek()


# # SunPy Maps!
# Maps are the same for each file, **regardless of source**. It does not matter if the source is SDO or SOHO, for example.
# 
# The most used attributes are as follows (some of them will look similar to NumPy's Array):

# In[ ]:


aia.data


# The data (stored in a numpy array)

# In[ ]:


type(aia.data)


# In[ ]:


aia.mean(),aia.max(),aia.min()


# Because it is just a numpy array you have access to all of those function

# The standard deviation

# In[ ]:


aia.data.std()


# In[ ]:


aia.data.shape


# The original metadata (stored in a dictionary)

# In[ ]:


aia.meta


# In[ ]:


aia.meta.keys()


# In[ ]:


aia.meta.get('rsun_obs')


# We also provide quick access to some key metadata values as object variables (these are shortcuts)

# In[ ]:


print(aia.date, aia.coordinate_system, aia.detector, aia.dsun)


# Maps also provide some nice map specific functions such as submaps. Let's zoom in on the flare location which was given to us by the HEK.

# In[ ]:


result[0].get('hpc_coord')


# In[ ]:


point = [665.04, -233.4096]
dx = 50
dy = 50
xrange = [point[0] - dx, point[0] + dx]
yrange = [point[1] - dy, point[1] + dy]
aia.submap(xrange,yrange).peek()
plt.plot(point[0], point[1], '+')
plt.xlim(xrange)
plt.ylim(yrange)


# The default image scale is definitely not right. Let's fix that so we can see the flare region better.

# In[ ]:


smap = aia.submap(xrange,yrange)

import matplotlib.colors as colors
norm = colors.Normalize(0, 3000)
smap.plot(norm=norm)
plt.plot(point[0], point[1], '+')
smap.draw_grid(grid_spacing=1)
plt.colorbar()


# ## Composite Maps
# Let's plot the two channels we downloaded from AIA together! Composite map is the way to do this. Let's check out the other map we got.

# In[ ]:


aia131 = Map(f[0])
aia131


# In[ ]:


smap131 = aia131.submap(xrange, yrange)
smap131.peek()


# In[ ]:


norm = colors.Normalize(0, 4000)
smap131.plot(norm=norm)
plt.colorbar()


# In[ ]:


smap171 = smap


# In[ ]:


compmap = Map(smap171, smap131, composite=True)


# In[ ]:


levels = np.arange(0,100,5)
print(levels)


# In[ ]:


compmap.set_levels(1, levels, percent=True)
compmap.set_mpl_color_normalizer(0, norm)
compmap.set_colors(1, plt.cm.Reds)


# In[ ]:


compmap.plot(norm=norm)
plt.show()


# Some other topics...

# ## Solar Constants

# In[ ]:


from sunpy.sun import constants as solar_constants


# In[ ]:


solar_constants.mass


# In[ ]:


print(solar_constants.mass)


# In[ ]:


(solar_constants.mass/solar_constants.volume).cgs


# In[ ]:


solar_constants.volume + solar_constants.density


# Not all constants have a short-cut assigned to them (as above). The rest of the constants are stored in a dictionary. The following code grabs the dictionary and gets all of the keys.:

# In[ ]:


solar_constants.physical_constants.keys()


# In[ ]:


type(solar_constants.mass)


# These are astropy constants which are a subclass of Quantities (numbers with units) which are a great idea.

# In[ ]:


from astropy import units as u


# In[ ]:


u.keV


# In[ ]:


u.keV.decompose()


# This has been a simple overview of AstroPy units. IF you want to read more, see http://astropy.readthedocs.org/en/latest/units/.

# ## More References
# 
# + [Python for Data Analysis](http://shop.oreilly.com/product/0636920023784.do) by Wes McKinney is great
# + [Lectures on scientific computing with Python](http://jrjohansson.github.io) Great lectures on scientific computing in Python

# ## Consider [contributing](http://sunpy.org/contribute/) to SunPy!
# 
# + Provide feedback on the [mailing list](https://groups.google.com/forum/#!forum/sunpy)
# + Report [bugs](https://github.com/sunpy/sunpy/issues)
# + Provide Code (see our [developer guide](http://sunpy.readthedocs.org/en/stable/dev.html))