<jupyter_start><jupyter_text>Solar Data Processing with Python Part IINow we have a grasp of the basics of python, but the whole reason for downloading python in the first place was to analyze solar data. Let's take a closer look at examples of solar data analysis. We will be using SunPy to access solar data. SunPy is a python package designed to interface between the powerful tools that exist in other Python Libraries with current repositories of solar data. With SunPy we will show how to: download solar data sets from the VSO, calibrate to industry standards, plot and overlay a time series.   Fitting A Gaussian to a Spectral LineOne of the most common data types in solar data processing is a time series. A time series is a measurement of how one physical parameter changes as a function of time. This example shows how to fit a gaussian to a spectral line. In this example, it will be as "real world" as possible.First, let's import some useful libraries.<jupyter_code>from datetime import datetime, timedelta #we saw these in the last tutorial

import numpy as np 
from astropy.io import fits #how to read .fits files
from astropy.modeling import models, fitting #some handy fitting tools from astropy
import matplotlib.pyplot as plt 
from scipy.integrate import trapz #numerical itegration tool
import astropy.units as u #units!!<jupyter_output><empty_output><jupyter_text>Next we need to load in the data set we want to work with:<jupyter_code>filename = '/Users/mskirk/Documents/Conferences/SDO 2016/EVS_L2_2011045_01_005_01.fits'
hdulist = fits.open(filename)<jupyter_output><empty_output><jupyter_text>So what did we get when we opened the file? Let's take a look:<jupyter_code>len(hdulist)<jupyter_output><empty_output><jupyter_text>We got 4 items in the list. Lets take a look at the first one:<jupyter_code>hdulist[0].header<jupyter_output><empty_output><jupyter_text>This doesn't contain any useful information. Let's take a look at the second item:<jupyter_code>hdulist[1].header<jupyter_output><empty_output><jupyter_text>Alright, now we are getting somewhere. This has wavelength information in units of 'nm' and accuracy information without units. Let's take a look at the other elements of the list we got:<jupyter_code>hdulist[2].header
hdulist[3].header<jupyter_output><empty_output><jupyter_text>So it looks like we are working with some wavelength data, spectral information, irradiance data, etc.   Plotting Spectral DataLet's take a look at some of the data we've got.<jupyter_code>len(hdulist[1].data)
hdulist[1].data["WAVELENGTH"]
len(hdulist[3].data)
hdulist[3].data[0]["IRRADIANCE"]
plt.plot(hdulist[1].data["WAVELENGTH"], hdulist[3].data[0]["IRRADIANCE"])
plt.ylim(0, 0.007)
plt.xlim(5, 40)<jupyter_output><empty_output><jupyter_text>So now we have a plot of wavelength vs. irradiance. We can see there is one major spike in the data. Let's filter the data so that we just have that one spike.<jupyter_code>w = np.logical_and(hdulist[1].data["WAVELENGTH"] > 30.2, hdulist[1].data["WAVELENGTH"] < 30.6)
w<jupyter_output><empty_output><jupyter_text>This function, "np.logical_and", is similar to a "where" statement in IDL. We can see that "w" is now an array of true and false values. To take a subsection of our data where our filter is true:<jupyter_code>wavelength = hdulist[1].data["WAVELENGTH"][w]
wavelength<jupyter_output><empty_output><jupyter_text>Now, we need to add some units to this data. The header of the file tells us what the units are:<jupyter_code>wavelength = u.Quantity(wavelength, unit="nm")<jupyter_output><empty_output><jupyter_text>Let's do the same thing to the irradiance data.<jupyter_code>irr = hdulist[3].data[0]["IRRADIANCE"][w]
irr = u.Quantity(irr, unit="W/m**2")<jupyter_output><empty_output><jupyter_text>What have we got? Let's plot it and take a look:<jupyter_code>plt.plot(wavelength, irr, 'o')<jupyter_output><empty_output><jupyter_text>Fit He II 304 line with a GaussianNow that we have extracted the He II line from our total spectrum, we want to fit it with a gaussian. Do do this we will make use of a couple of packages in in astropy. We will initialize the gaussian fit with some approximations (max, center, FWHM):<jupyter_code>g_init = models.Gaussian1D(0.007, 30.4, .2)<jupyter_output><empty_output><jupyter_text>Now let's define a fitting method and produce a fit:<jupyter_code>fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, wavelength, irr)<jupyter_output><empty_output><jupyter_text>Let's take a look at some of the qualities of our fitted gaussian:<jupyter_code>g.mean
g.stddev
g.amplitude
g<jupyter_output><empty_output><jupyter_text>Our guesses wern't too bad, but we over estimated the Standard Deviation by about a factor of 10. The variable 'g' has the fitted parameters of our gaussian but it doesn't actually contain an array. To plot it over the data, we need to create an array of values. We will make an array from 30.2 to 30.6 with 1000 points in it.<jupyter_code>x = np.linspace(30.2, 30.6, 1000)<jupyter_output><empty_output><jupyter_text>To find the values of our fit at each location, it is easy:<jupyter_code>y = g(x)<jupyter_output><empty_output><jupyter_text>Now we can plot it:<jupyter_code>plt.plot(wavelength, irr, 'o')
plt.plot(x, y)<jupyter_output><empty_output><jupyter_text>Ingegrating under the curve. Let's find the area under the curve we just created. We can numerically integrate it easily:<jupyter_code>intensity = trapz(y,x)
intensity<jupyter_output><empty_output><jupyter_text>Produce a Light Curve of He II 304 From All Spectra in FileThe file we downloaded is a hyper-spectrum. This means that the spectrum changes over time. Can we find how the intensity of the line changes over time using the same fitting tools that we just showcased? Sure, we just need to put everything into a loop.<jupyter_code>x = np.linspace(30.2, 30.6, 1000) #define our wavelength array
intensity = [] # initialize the intensity array
w = np.logical_and(hdulist[1].data["WAVELENGTH"] > 30.2, hdulist[1].data["WAVELENGTH"] < 30.6) # isolate the 304 line
wavelength = hdulist[1].data["WAVELENGTH"][w] # filter all of the spectra
fit_g = fitting.LevMarLSQFitter() #fitting algorithm

for spectrum in hdulist[3].data: #looping over all the spectra (360)
    irr = spectrum["IRRADIANCE"][w] #pulling out the irradiance 
    g_init = models.Gaussian1D(0.007, 30.4, .2) #initialize the gaussian fit
    g = fit_g(g_init, wavelength, irr) #do the fitting 
    y = g(x) # create the curve
    intensity.append(trapz(y,x)) # aggregate the results
intensity = u.Quantity(intensity, unit="nm") # adding units. 
plt.plot(intensity)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0889556724267468, 'fraction_numerical': 0.018632607062359127, 'mean_word_length': 5.883143743536712, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 83, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '486394', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2303, 'n_tokens_neox': 2153, 'n_words': 792, 'path': 'SDO2016_Py3.ipynb', 'repo_name': 'MSKirk/tutorials', 'size': '6655'}