In this notebook, we will be plotting keywords and images, from data taken by the Helioseismic and Magnetic Imager (HMI) instrument on NASA's Solar Dynamics Observatory (SDO) satellite. SDO takes about a terabyte and a half of data a day, which is more data than any other satellite in NASA history.
Data from the HMI and Atmospheric Imaging Assembly (AIA) instruments aboard SDO are stored at Stanford University. The metadata are stored in a pSQL database called the Data Record Management System, or DRMS. The image data are stored separately, in storage units called the Storage Unit Management System, or SUMS. Data are merged together, upon export from both systems, as FITS files. DRMS and SUMS together constitute the Joint Science Operations Center, or JSOC.
The easiest way to access SDO HMI and AIA data is via the python
drms module, available at PyPI. In addition to the numerous tutorials on both the Read the Docs and Github, all the examples below utilize the
drms module. First we'll import the module, and some others:
import drms import json, numpy as np, matplotlib.pylab as plt, matplotlib.ticker as mtick from datetime import datetime as dt_obj import urllib from astropy.io import fits from sunpy.cm import color_tables as ct from matplotlib.dates import * import matplotlib.image as mpimg from IPython.display import Image %matplotlib inline %config InlineBackend.figure_format='retina'
The first step in querying for SDO HMI and AIA data is to establish a connection to JSOC. This can be done with the
Client() class.
import drms c = drms.Client()
The
Client() class allows one to access both metadata and image data simultaneously via a data series. A data series contains all of particular type of data — e.g. there is a series for continuum intensity data, another for magnetic field data, and so forth. Read Section 4 of the SDO Analysis Guide for more information about how to build a data series query. For example, to find all the SHARP data series, execute the following regular expression query:
c.series(r'hmi\.sharp_')
['hmi.sharp_720s', 'hmi.sharp_720s_nrt', 'hmi.sharp_cea_720s', 'hmi.sharp_cea_720s_nrt']
The query below retrieves both metadata and image data for active region 11158, which produced an X2.2-class flare on February 15, 2011 at 1:56 UT, from the SHARP data series. The SHARP data series include patches of vector magnetic field data taken by the HMI instrument. These patches encapsulate automatically-detected active regions that are tracked throughout the entirety of their disk passage. The
c.query() method takes three arguments:
hmi.sharp_cea_720s, appended with two prime keys: the HARP number (377) and the time range (2011.02.14_15:00:00/12h). A prime key, or set of prime keys, is a unique identifier.
hmi.sharp_cea_720s:
si = c.info('hmi.sharp_cea_720s') si.keywords
hmi.sharp_cea_720s:
si = c.info('hmi.sharp_cea_720s') si.segments
keys, segments = c.query('hmi.sharp_cea_720s[377][2011.02.14_15:00:00/12h][? (QUALITY<65536) ?]', key='T_REC, USFLUX, ERRVF', seg='Br')
To convert the keyword
T_REC into a datetime object, we can use the function below.
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
t_rec = np.array([parse_tai_string(keys.T_REC[i],datetime=True) for i in range(keys.T_REC.size)])
Now for some plotting! matplotlib.pyplot generates two objects: a figure and axes. The data are ascribed to the axes. The time axes in particular requires some formatting; in order to free it of clutter, we'll plot tick marks every three hours and label them accordingly.
fig, ax = plt.subplots(figsize=(8,7)) # define the size of the figure orangered = (1.0,0.27,0,1.0) # create an orange-red color # define some style elements marker_style = dict(linestyle='', markersize=8, fillstyle='full',color=orangered,markeredgecolor=orangered) text_style = dict(fontsize=16, fontdict={'family': 'monospace'}) ax.tick_params(labelsize=14) ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2f')) # ascribe the data to the axes ax.plot(t_rec, (keys.USFLUX)/(1e22),'o',**marker_style) ax.errorbar(t_rec, (keys.USFLUX)/(1e22), yerr=(keys.ERRVF)/(1e22), linestyle='',color=orangered) # format the x-axis with universal time locator = AutoDateLocator() locator.intervald[HOURLY] = [3] # only show every 3 hours formatter = DateFormatter('%H') ax.xaxis.set_major_locator(locator) ax.xaxis.set_major_formatter(formatter) # set yrange ax.set_ylim([2.4,2.8]) # label the axes and the plot ax.set_xlabel('time in UT',**text_style) ax.set_ylabel('maxwells x 1e22',**text_style) ax.set_title('total unsigned flux starting at '+str(t_rec[0])+' UT',**text_style) # annotate the plot with a start time
Text(0.5, 1.0, 'total unsigned flux starting at 2011-02-14 15:00:00 UT')
Since we also queried for the image data — specifically, for the radial component of the magnetic field — we can view these data by downloading the unmerged FITS file, which means that the header metadata has not been merged with the image data, via
astropy. This is the quickest and easiest option as the
drms.Client() class can query the header and image data independently, there is no need to store FITS files locally, and there is no need to wait for the exportdata system to generate FITS files.
url = '' + segments.Br[54] # add the jsoc.stanford.edu suffix to the segment name photosphere_image = fits.open(url) # download the data
hmimag = plt.get_cmap('hmimag') plt.imshow(photosphere_image[1].data,cmap=hmimag,origin='lower',vmin=-3000,vmax=3000) print('The dimensions of this image are',photosphere_image[1].data.shape[0],'by',photosphere_image[1].data.shape[1],'.')
The dimensions of this image are 377 by 744 .
The JSOC allows users to export data as FITS files, jpg images and movies in mp4 or mpg format.
n.b. The code below will only work with a valid e-mail address. In order to obtain one, users must register on the JSOC exportdata website. It is worthwhile to have a read around that webpage before using the
drms module to export data.
c = drms.Client(email=email, verbose=True)
We can easily export any image data as a jpg using one of several color tables (see the expordata website under the jpg protocol for a list of color tables). Here is an example with NOAA Active Region 11158 at the time of the X2.2-class flare.
protocol_args = {'ct': 'HMI_mag.lut', 'min': -3000, 'max': 3000}
r = c.export('hmi.sharp_cea_720s[377][2011.02.14_02:00:00]{Br}', protocol='jpg', protocol_args=protocol_args, email=email)
print(r.request_url)
Export request pending. [id="JSOC_20200728_2607", status=2] Waiting for 5 seconds... Export request pending. [id="JSOC_20200728_2607", status=1] Waiting for 5 seconds...
image_url = r.urls['url'][0] print('The image is at this url: ',image_url,'.')
The image is at this url: .
We can also easily create a movie using any image data! The movie below tracks NOAA Active Region 11158 over the same time interval as the plot above. In this case, the segment information, Br, is notated in the curly brackets within the data series query.
r = c.export('hmi.sharp_cea_720s[377][2011.02.14_15:00:00/12h]{Br}', protocol='mpg', protocol_args=protocol_args, email=email)
print(r.request_url)
Export request pending. [id="JSOC_20200728_2608", status=2] Waiting for 5 seconds... Export request pending. [id="JSOC_20200728_2608", status=1] Waiting for 5 seconds... Export request pending. [id="JSOC_20200728_2608", status=1] Waiting for 5 seconds...
print("You can download the movie from here:", r.urls['url'][0])
You can download the movie from here: