PyBats Tutorial: RIM

A colleague of mine asked for a Spacepy/pybats plotting script for looking at results from the Ridley Ionosphere Model, or RIM*. While it’s easy enough to hand someone a quick script to whip up a set of plots, it kind of misses the point of PyBats within Spacepy. Rather than just act as a way to read SWMF output, including RIM files, PyBats is designed to enable interactive analysis and plotting of output and make generation of scripts a trivial task. PyBats is an API, not just a group of scripts. Let’s do an interactive tutorial to demonstrate this.

*Note that “RIM” is shorthand for two separate models: RIM proper, an advanced and flexible ionosphere model, and “Ridley_serial”, the legacy ionosphere model within the Space Weather Modeling Framework and most widely used of the two. Usually, when someone says “RIM”, they mean the latter, even if they don’t know it.

Start with an example data file, which can be found here. This is a typical ionospheric electrodynamics output file from a simple simulation of the magnetosphere-ionosphere-ring current system under purely southward interplanetary magnetic field. Next, crack open IPython, which can/should be installed either via a package manager (including PIP) or via your Anaconda install. IPython makes interactive Python work easier. If you’re not using it, you should start right now.

The basics: Opening and Exploring

How do we open this file? How do we play with what’s inside? Here, we turn to the first two golden rules of PyBats:

1) Always use the submodule that corresponds to the model in question.
2) Always use the object that matches your file type if possible.

In this case, the model in question is RIM- we want to look at RIM output. Next, we’re looking at basic ionosphere output, so we want an “Iono” object. In Python, this will look like so:

# Import the correct submodule: 
from spacepy.pybats import rim 

#Open our file as an Iono object: 
data = rim.Iono('./it150321_032544_000.idl') 

The critical way to learn what’s is available is via the IPython terminal- use tab complete to truly explore what submodules, functions, and classes are available. Once you have an object instantiated (i.e., the file is opened as an object), use tab-complete to explore object attributes and methods.

Our data object is a Spacepy datamodel object, which means that it is easy to explore the values inside. Fundamentally, it works like a dictionary with key-value pairs to store the data inside.

 
data.keys()

 
dict_keys(['n_theta', 's_theta', 'n_psi', 's_psi', 'n_sigmah', 's_sigmah', 'n_sigmap', 's_sigmap', 'n_jr', 's_jr', 'n_phi', 's_phi', 'n_e-flux', 's_e-flux', 'n_ave-e', 's_ave-e', 'n_rt 1/b', 's_rt 1/b', 'n_rt rho', 's_rt rho', 'n_rt p', 's_rt p', 'n_jouleheat', 's_jouleheat', 'n_ionnumflux', 's_ionnumflux', 'n_conjugate dlat', 's_conjugate dlat', 'n_conjugate dlon', 's_conjugate dlon'])

 
data['n_jr']

dmarray([[0.0048284, 0.0048297, 0.004831 , ..., 0.0048258, 0.0048271,
          0.0048284],
         [0.004186 , 0.0044131, 0.0046408, ..., 0.0037345, 0.0039597,
          0.004186 ],
         [0.0037367, 0.0041913, 0.0046464, ..., 0.0028315, 0.0032833,
          0.0037367],
         ...,
         [0.       , 0.       , 0.       , ..., 0.       , 0.       ,
          0.       ],
         [0.       , 0.       , 0.       , ..., 0.       , 0.       ,
          0.       ],
         [0.       , 0.       , 0.       , ..., 0.       , 0.       ,
          0.       ]])

In the first command, we see a list of the variables stored within the data object. Note how, for Iono objects, the prefixes n_ and s_ refer to the northern and southern hemispheres, respectively. Once we see what’s available, we can access each value by name as we would a dictionary. The data inside acts like a Numpy array, so we can slice and manipulate it as such. The usual methods, such as mean and max, all apply here. Use tab-complete in IPython to see what methods are available!

Unlike Numpy arrays, however, Spacepy’s dmarray objects (short for data model arrays) have attributes associated with them. Access them via the attrs object attributes:

 
#Start by looking at the object's file-level attributes:
data.attrs

 
{'file': './it150321_032544_000.idl',<br>'title': 'BATSRUS: Ionospheric Potential Solution, 2015-03-21-03-25-44-000, BTiltDeg= 0.00 0.00',<br>'nvars': 15,<br>'ntheta': 91,<br>'nphi': 181,<br>'time': datetime.datetime(2015, 3, 21, 3, 25, 44),<br>'iter': 5145,<br>'simtime': 5144.0}

 
# Now look at attributes of a single variable:
data['n_jr'].attrs

 
{'units': 'mA/m^2'}

Creating New Values

Because Spacepy and Pybats objects inherit from dictionaries, adding more data works as you would expect– just add more key-data pairs.

 
# This doesn't make physics sense, but that's fine.
data['n_sigmaAll'] = data['n_sigmah']+data['n_sigmap']

# We can do an easy unit change...
data['n_phi_V'] = 1000*data['n_phi']

#But it's good to note the unit change:
data['n_phi_V'].attrs['units'] = 'V'

Spacepy datamodel objects work like dictionaries. Add new key-value pairs when creating new variables.

Some common calculations come built-in. Look for calc_ methods. The easiest way in IPython is to type data.calc then mash the tab button and let tab complete show you what’s available. The docstrings, accessible via the question mark syntax (e.g., data.calc_I?), will fill you in on the details. For rim.Iono objects, there’s currently only one calculation function*: calc_I, which calculates integrated field-aligned current values across each hemisphere.

*This is changing in some of the development branches, with new functions to calculation azimuthal currents and more.

 
# Call the calculation function:
data.calc_I()

# Look for the new values within the object:
data.keys()

 
dict_keys([... 'n_I', 'n_Iup', 'n_Idown', 's_I', 's_Iup', 's_Idown'])

 
# Look at the new values:
data['s_Iup']

 
0.92346768990336

Plotting and Visualization

Users can just pass values to Matplotlib functions and classes to create plots. However, this is typically a very time intensive process, so many common plotting tasks are customized within Pybats objects. Let’s start with a quick example and then talk about what’s happening.

 
# Import matplotlib and turn on interactive plotting:
from matplotlib import pyplot as plt
plt.ion()

# Create a contour plot:
data.add_cont('n_jr', n=51, add_cbar=True)

 
(<Figure size 640x480 with 2 Axes>,
 <matplotlib.axes._subplots.PolarAxesSubplot at 0x7f9eeef96c70>,
 <matplotlib.contour.QuadContourSet at 0x7f9eef0e3940>,
 <matplotlib.colorbar.Colorbar at 0x7f9eef0cf9a0>)

A very quick and easy plot of the field aligned currents! In Pybats, most built-in plotting routines start with add_. Use tab-complete to quickly find the available built-in plotting routines. Note how the plotting method returned many objects: the figure, the subplot, the contour, and the color bar. That allows for further customization after the initial plot is created.

Use object.add_ and tab-complete to quickly find the available built-in plotting routines.

Let’s get slightly more complicated and compare the northern and southern hemisphere on a single plot with a shared color bar:

 
# Create a figure:
fig = plt.figure()

# Add the contours to the figure.
# Note how we control the plot range and destination location:
out1 = data.add_cont('n_jr', zmax=.2, target=fig, loc=121, label='North')
out2 = data.add_cont('s_jr', zmax=.2, target=fig, loc=122, label='South')

# Add a colorbar to the figure:
fig.colorbar(out1[2], cax=fig.add_axes([.25, .08, .5, .05]), orientation='horizontal', label='FAC')

This plot could use some polishing for sure, including changing the awful default tick spacing on color bar using the various Matplotlib ticker classes. But it does get the job done in relatively short order. Let’s break down the critical pieces:

• The target and loc syntax are omnipresent in Spacepy. target sets the location where the plot will be placed and can be a figure or an axes object. If it’s a figure object, loc is a matplotlib-style location identifier (for example, 121 means one row, two columns, and the plot goes in the first location). Use target and loc to make customized multi-panel plots easily.
• Each add_* type method has a ton of keyword arguments for customization. Here, we used zmax to set the plot range and label to set the label at the top of the plot. Explore the docstrings to learn more.
• Note how the contour objects resulting from the plotting methods is saved within out1 and out2 and later used to create the color bar at the bottom of the plot. Pybats’ add_* methods will return lists of objects corresponding to the components of the drawn plot. Knowing how to work with Matplotlib after the initial plot is created is a critical component to good looking figures.

This type of behavior is found in all of the add_* methods inside of Pybats. All of this leads us to our 3rd golden rule of PyBats, which applies everywhere across Python:

3) Exploration yields expertise

Keep digging to find new classes, object methods, and keyword arguments. What you want to do might already be in the library! If not, throw a feature request over on GitHub and we’ll try to get it in the next distribution.

From Prototype to Script

At this point, it’s just a matter of fine-tuning your plots until you get it right where you want it. Then, using the history command, copy-pasta the relevant commands into a script, add a for loop and some glob action and you have a bona-fide plotting script to analyze your simulation. A good example of this final step can be found over on this post here.

There are certainly more details, but these are more module- and output-specific. I’ll cover those in future posts.