plotting the sun's polar field with python + d3

In this notebook, we will be plotting the average magnetic field around the solar poles to determine when, exactly, the Sun's polar field flips. We'll use data from the Helioseismic and Magnetic Imager instrument on NASA's Solar Dynamics Observatory (SDO) satellite, which takes about 1.5 terabytes of data a day (more data than any satellite in NASA history). We'll grab the data from the JSOC database via a JSON API and plot it using mpld3, a python interface for d3.

In [1]:
import json, urllib, requests, numpy as np, matplotlib.pylab as plt, matplotlib.ticker as mtick, mpld3
from mpld3 import plugins
from datetime import datetime as dt_obj
from matplotlib.dates import *
%matplotlib inline
%config InlineBackend.figure_format='retina'

Now we'll execute a sample JSON query using the json and urllib modules to gather the data -- specifically, two keywords that hold the value of the mean radial component of the magnetic field in the polar region (greater than 60 degrees latitude). In theory, it is possible to parse JSON queries using the read_json() function in the pandas library, but, in practice, the object returned by the JSON API for the SDO database, jsoc_info, isn't formatted in a way that pandas.read_json() can understand easily.

In [2]:
# first get the most recent timestamp, or t_rec
url = "http://jsoc.stanford.edu/cgi-bin/ajax/jsoc_info?ds=hmi.meanpf_720s[$]&op=rs_list&key=T_REC"
response = requests.get(url)
data = response.json()
t_last = data['keywords'][0]['values'][0]

# now gather all the data between SDO launch and the most recent t_rec
url = "http://jsoc.stanford.edu/cgi-bin/ajax/jsoc_info?ds=hmi.meanpf_720s[2010.05.01_00_TAI-"+t_last+"@12h][? QUALITY=0 ?]&op=rs_list&key=T_REC,CAPN2,CAPS2"
response = requests.get(url)
data = response.json()

Now we'll create some empty lists to hold the data:

In [3]:
capn2 = [] # this holds the keyword CAPN2 (Mean radial field in N60-N90) 
caps2 = [] # this holds the keyword CAPS2 (Mean radial field in S60-S90)
t_rec  = [] # this holds the keyword T_REC

data is of type dict, so we'll get the number of keyword elements this way:

In [4]:
n_elements = len(data['keywords'][1]['values'])

And then we'll populate our empty lists to hold the keyword values:

In [5]:
for i in range(n_elements):
    capn2.append(float(data['keywords'][1]['values'][i]))

for i in range(n_elements):
    caps2.append(float(data['keywords'][2]['values'][i]))

for i in range(n_elements):
    t_rec.append(data['keywords'][0]['values'][i])

Let's convert our lists to numpy arrays, because they're easier to work with:

In [6]:
capn2 = np.array(capn2)
caps2 = np.array(caps2)
t_rec  = np.array(t_rec,dtype='S16') # this dtype is ok to set as the format of T_REC will never change

The error in the value of CAPN2 and CAPS2 at any point in time is defined by the standard deviation in this quantity over a 30 day period (beginning 15 days before the time of interest and ending 15 days after the time of interest):

In [7]:
err_capn2 = np.zeros(n_elements)
err_caps2 = np.zeros(n_elements)

for i in range(0,n_elements):
    err_capn2[i] = np.std(capn2[i-30:i+30])
    
for i in range(0,n_elements):
    err_caps2[i] = np.std(caps2[i-30:i+30])

T_REC is formatted in a way that matplotlib.pyplot will not understand, so let's convert the numpy array into a datetime object:

In [8]:
def parse_tai_string(tstr,datetime=True):
    year   = int(tstr[:4])
    month  = int(tstr[5:7])
    day    = int(tstr[8:10])
    hour   = int(tstr[11:13])
    minute = int(tstr[14:16])
    if datetime: return dt_obj(year,month,day,hour,minute)
    else: return year,month,day,hour,minute

x = np.array([parse_tai_string(t_rec[i],datetime=True) for i in range(t_rec.size)])

Now for some interactive plotting with mpld3! The home, pan, and zoom buttons on the lower left-hand corner of the plot allow you to explore the data.
To throw the plot below up on a webpage, simply print the html by calling the mpld3 function save_html().

In [9]:
fig, ax = plt.subplots(figsize=(10,8))      # define the size of the figure
orangered = (1.0,0.27,0,1.0)                # create an orange-red color
cornblue  = (0.39,0.58,0.93,1.0)            # create an cornflower blue color

# define some style elements
marker_style = dict(linestyle='-', linewidth=4, fillstyle='full',color=orangered,markeredgecolor=orangered)
marker_style2 = dict(linestyle='-', linewidth=4, fillstyle='full',color=cornblue,markeredgecolor=cornblue)
text_style = {'color':'black','weight':'normal','size': 16}

ax.tick_params(labelsize=14)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))

# ascribe the data to the axes
ax.plot(x, capn2, label='north pole', **marker_style)
ax.plot(x, caps2, label='south pole', **marker_style2)
ax.set_xlim(x[0],x[-1])
ax.set_ylim([-6,6])

# add the error bars
ax.fill_between(x, capn2-err_capn2, capn2+err_capn2, facecolor=orangered, edgecolor='None', alpha=0.4, interpolate=True)
ax.fill_between(x, caps2-err_caps2, caps2+err_caps2, facecolor=cornblue, edgecolor='None', alpha=0.4, interpolate=True)

# format the x-axis with universal time
locator = AutoDateLocator()
locator.intervald[DAILY] = [1] # only show every year
formatter = DateFormatter('%Y')
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)

# label the axes and the plot
ax.set_xlabel('time', fontdict = text_style)
ax.set_ylabel('mean radial field strength in gauss', fontdict=text_style, labelpad = 5)
ax.set_title('mean polar field', fontdict = text_style) # annotate the plot with a start time
fig.set_size_inches(10,4)
legend = plt.legend(loc='upper right', fontsize=12, framealpha=0.0,title='')
legend.get_frame().set_linewidth(0.0)

Now we'll smooth the plot by plotting the running average of CAPN2 and CAPS2 over a 30 day window:

In [10]:
avg_capn2 = np.zeros(n_elements)
avg_caps2 = np.zeros(n_elements)

for i in range(0,n_elements):
    avg_capn2[i] = np.mean(capn2[i-30:i+30])
    
for i in range(0,n_elements):
    avg_caps2[i] = np.mean(caps2[i-30:i+30])
In [11]:
fig, ax = plt.subplots(figsize=(10,8))      # define the size of the figure
orangered = (1.0,0.27,0,1.0)                # create an orange-red color
cornblue  = (0.39,0.58,0.93,1.0)            # create an cornflower blue color

# define some style elements
marker_style = dict(linestyle='-', linewidth=4, fillstyle='full',color=orangered,markeredgecolor=orangered)
marker_style2 = dict(linestyle='-', linewidth=4, fillstyle='full',color=cornblue,markeredgecolor=cornblue)
text_style = {'color':'black','weight':'normal','size': 16}

ax.tick_params(labelsize=14)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))

# ascribe the data to the axes
ax.plot(x, avg_capn2, label='north pole', **marker_style)
ax.plot(x, avg_caps2, label='south pole', **marker_style2)
ax.set_ylim([-6,6])
ax.set_xlim(x[0],x[-1])

# add the error bars
ax.fill_between(x, avg_capn2-err_capn2, avg_capn2+err_capn2, facecolor=orangered, edgecolor='None', alpha=0.4, interpolate=True)
ax.fill_between(x, avg_caps2-err_caps2, avg_caps2+err_caps2, facecolor=cornblue, edgecolor='None', alpha=0.4, interpolate=True)

# format the x-axis with universal time
locator = AutoDateLocator()
locator.intervald[DAILY] = [1] # only show every year
formatter = DateFormatter('%Y')
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)

# label the axes and the plot
ax.set_xlabel('time', fontdict = text_style)
ax.set_ylabel('mean radial field strength in gauss', fontdict=text_style, labelpad = 5)
ax.set_title('mean polar field', fontdict = text_style) # annotate the plot with a start time
fig.set_size_inches(10,4)
legend = plt.legend(loc='upper right', fontsize=12, framealpha=0.0,title='')
legend.get_frame().set_linewidth(0.0)

To see a daily-updated and interactive version of this plot, check out the HMI Polar Field site. To create your own html file, use the following command:

In [12]:
mpld3.save_html(fig,'12hcadence_interp.html')