'''
A PyBats module for handling input, output, and visualization of
binary SWMF output files taylored to BATS-R-US-type data.
.. currentmodule:: spacepy.pybats.bats
.. rubric:: Classes
.. autosummary::
:template: clean_class.rst
:toctree:
BatsLog
Stream
Bats2d
Mag
MagFile
GeoIndexFile
VirtSat
'''
import sys
import numpy as np
from spacepy.pybats import PbData, IdlFile, LogFile, calc_wrapper
import spacepy.plot.apionly
from spacepy.plot import set_target, applySmartTimeTicks
from spacepy.datamodel import dmarray
#### Module-level variables:
# recognized species:
mass = {'hp':1.0, 'op':16.0, 'he':4.0,
'sw':1.0, 'o':16.0, 'h':1.0, 'iono':1.0, '':1.0}
RE = 6371000 # Earth radius in meters.
#### Module-level functions:
def _calc_ndens(obj):
'''
Given an object of subclass :class:`~spacepy.pybats.PbData` that uses
standard BATS-R-US variable names, calculate the number density from
all mass density variables (i.e., all variables that end in *rho*).
New variables are added to the object. Variable names are constructed
by taking the mass-density variable, *Speciesrho*, and replacing *rho*
with *N*. Total number density is also saved as *N*.
Composition information is also saved by taking each species and
calculating the percent of the total number density via
*fracspecies* = 100% x *speciesN*/*N*.
Many species that are implemented in multi-species and multi-fluid
BATS-R-US equation files have known mass. For example, *oprho* is known
a-priori to be singly ionized oxygen mass density, so *opN* is created
by dividing *oprho* by 16. If the fluid/species is not recognized, it is
assumed to be hyrogen. The single atom/molecule mass is saved in the
attributes of the new variable.
This function should be called by object methods of the same name.
It is placed at the module level because it is used by many different
classes.
Parameters
==========
obj : :class:`~spacepy.pybats.PbData` object
The object on which to act.
Other Parameters
================
Returns
=======
True
Examples
========
>>> import spacepy.pybats.bats as pbs
>>> mhd = pbs.Bats2d('spacepy-code/spacepy/pybats/slice2d_species.out')
>>> pbs._calc_ndens(mhd)
'''
species = []
names = []
# Get name of Rho (case sensitive check):
rho = 'Rho'*('Rho' in obj) + 'rho'*('rho' in obj)
# Find all species: the variable names end or begin with "rho".
# Take care not to double-count.
for k in obj:
# Ends with rho?
if (k[-3:] == rho) and (k != rho) and (k[:-3] + 'N' not in obj):
species.append(k)
names.append(k[:-3])
# Begins with rho?
if (k[:3] == rho) and (k != rho) and (k[3:] + 'N' not in obj):
species.append(k)
names.append(k[3:])
# Individual number density
for s, n in zip(species, names):
if n.lower() in mass:
m = mass[n.lower()]
else:
m = 1.0
obj[n+'N'] = dmarray(obj[s]/m, attrs={'units':'$cm^{-3}$',
'amu mass':m})
# Total N is sum of individual number densities.
obj['N'] = dmarray(np.zeros(obj[rho].shape),
attrs={'units':'$cm^{-3}$'})
if species:
# Total number density:
for n in names: obj['N'] += obj[n+'N']
# Composition as fraction of total per species:
for n in names:
obj[n+'Frac'] = dmarray(100.*obj[n+'N']/obj['N'],
{'units':'Percent'})
else:
# No individual species => no composition, simple ndens.
obj['N'] += dmarray(obj[rho], attrs={'units':'$cm^{-3}$'})
#### Classes:
[docs]class BatsLog(LogFile):
'''
A specialized version of :class:`~spacepy.pybats.LogFile` that includes
special methods for plotting common BATS-R-US log file values, such as
D$_{ST}$.
.. autosummary::
~BatsLog.add_dst_quicklook
.. automethod:: add_dst_quicklook
'''
def fetch_obs_dst(self):
'''
Fetch the observed Dst index for the time period covered in the
logfile. Return *True* on success.
Observed Dst is automatically fetched from the Kyoto World Data Center
via the :mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoDst` object, which holds the observed
Dst, is stored as *self.obs_dst* for future use.
'''
import spacepy.pybats.kyoto as kt
# Return if already obtained:
if hasattr(self, 'obs_dst'): return True
# Start and end time to collect observations:
stime = self['time'][0]
etime = self['time'][-1]
# Attempt to fetch from Kyoto website:
try:
self.obs_dst = kt.fetch('dst', stime, etime)
# Warn on failure:
except BaseException as args:
raise Warning('Failed to fetch Kyoto Dst: ', args)
return False
return True
def fetch_obs_sym(self):
'''
Fetch the observed SYM-H index for the time period covered in the
logfile. Return *True* on success.
Observed SYM-H is automatically fetched from the Kyoto World Data Center
via the :mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoSym` object, which holds the observed
Dst, is stored as *self.obs_sym* for future use.
'''
import spacepy.pybats.kyoto as kt
# Return if already obtained:
if hasattr(self, 'obs_sym'): return True
# Start and end time to collect observations:
stime = self['time'][0]
etime = self['time'][-1]
# Attempt to fetch from Kyoto website:
try:
self.obs_sym = kt.fetch('sym', stime, etime)
# Warn on failure:
except BaseException as args:
raise Warning('Failed to fetch Kyoto Dst: ', args)
return False
return True
[docs] def add_dst_quicklook(self, target=None, loc=111, plot_obs=False,
epoch=None, add_legend=True, plot_sym=False,
dstvar=None, obs_kwargs={'c':'k', 'ls':'--'},
**kwargs):
'''
Create a quick-look plot of Dst (if variable present in file)
and compare against observations.
Like all *add_\* * methods in Pybats, the *target* kwarg determines
where to place the plot.
If kwarg *target* is **None** (default), a new figure is
generated from scratch. If *target* is a matplotlib Figure
object, a new axis is created to fill that figure at subplot location
*loc* (defaults to 111). If target is a matplotlib Axes object,
the plot is placed into that axis at subplot location *loc*.
With newer versions of BATS-R-US, new dst-like variables are included,
named 'dst', 'dst-sm', 'dstflx', etc. This subroutine will attempt
to first use 'dst-sm' as it is calculated consistently with
observations. If not found, 'dst' is used. Users may choose which
value to use via the *dstvar* kwarg.
Observed Dst and SYM-H is automatically fetched from the Kyoto World
Data Center via the :mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoDst` or
:class:`spacepy.pybats.kyoto.KyotoSym` object, which holds the observed
Dst/SYM-H, is stored as *self.obs_dst* for future use.
The observed line can be customized via the *obs_kwargs* kwarg, which
is a dictionary of plotting keyword arguments.
If kwarg *epoch* is set to a datetime object, a vertical dashed line
will be placed at that time.
The figure and axes objects are returned to the user.
'''
from datetime import datetime
import matplotlib.pyplot as plt
# Set the correct dst value to be used.
if not dstvar:
if 'dst_sm' in self:
dstvar = 'dst_sm'
else: dstvar='dst'
if dstvar not in self:
return None, None
fig, ax = set_target(target, figsize=(10, 4), loc=loc)
if 'label' not in kwargs:
kwargs['label'] = 'BATS-R-US $D_{ST}$ (Biot-Savart)'
if 'label' not in obs_kwargs:
obs_kwargs['label'] = 'Obs. Dst'
if plot_sym: obs_kwargs['label'] = 'Obs. SYM-H'
ax.plot(self['time'], self[dstvar], **kwargs)
ax.hlines(0.0, self['time'][0], self['time'][-1],
'k', ':', label='_nolegend_')
applySmartTimeTicks(ax, self['time'])
ax.set_ylabel('D$_{ST}$ ($nT$)')
ax.set_xlabel('Time from ' + self['time'][0].isoformat()+' UTC')
# Add observations (Dst and/or SYM-H):
if(plot_obs):
# Attempt to fetch observations, plot if success.
if self.fetch_obs_dst():
ax.plot(self.obs_dst['time'], self.obs_dst['dst'], **obs_kwargs)
applySmartTimeTicks(ax, self['time'])
if(plot_sym):
# Attempt to fetch SYM-h observations, plot if success.
if self.fetch_obs_sym():
ax.plot(self.obs_sym['time'], self.obs_sym['sym-h'], **obs_kwargs)
applySmartTimeTicks(ax, self['time'])
# Place vertical line at epoch:
if type(epoch) == datetime:
yrange = ax.get_ylim()
ax.vlines(epoch, yrange[0], yrange[1], linestyles='dashed',
colors='k', linewidths=1.5)
ax.set_ylim(yrange)
# Apply legend
if add_legend: ax.legend(loc='best')
if target is None: fig.tight_layout()
return fig, ax
class Extraction(PbData):
'''
A class for creating and visualizing extractions from other
:class:`~spacepy.pybats.PbData` 2D data sets. Bilinear interpolation is
used to obtain data between points. At present, this class only
works with :class:`~spacepy.pybats.bats.Bats2d` objects, but will be
generalized in the future. Though it can be instantiated in typical
fashion, It is best to create these objects via other object methods (e.g.,
:func:`~spacepy.pybats.bats.Bats2d.extract`)
Parameters
----------
x : float or sequence
X value(s) for points at which to extract data values.
y : float or sequence
Y value(s) for points at which to extract data values.
dataset : Bats
:class:`~spacepy.pybats.bats.Bats2d` object from which to extract values
Other Parameters
----------------
var_list : string or sequence of strings
List of values to extract from *dataset*. Defaults to 'all', for all
values within *dataset*.
'''
def __init__(self, x, y, dataset, var_list='all', *args, **kwargs):
# Initialize as a PbData object.
super(Extraction, self).__init__(*args, **kwargs)
# Store x, y internally:
# If our x, y locations are not numpy arrays, fix that.
# Additionally, convert scalars to 1-element vectors.
x, y = np.array(x), np.array(y)
if not x.shape: x = x.reshape(1)
if not y.shape: y = y.reshape(1)
# Default: all variables are extracted except coordinates.
if var_list == 'all':
var_list = list(dataset.keys())
for v in ('x','y','z','grid'):
if v in var_list: var_list.remove(v)
elif type(var_list)==str:
var_list = var_list.split()
elif type(var_list) != list:
var_list=[]
else:
raise TypeError('Kwarg var_list must be string, list, or None.')
# Stash var list, dataset internally:
self._var_list = var_list
self._dataset = dataset
# Extract along trace:
self.extract(x, y)
def extract(self, xpts, ypts):
# Perform actual extraction:
from spacepy.pybats.batsmath import interp_2d_reg
# Some convenience variables:
x, y = xpts, ypts
data = self._dataset
for v, value in zip(data['grid'].attrs['dims'], (x, y)):
self[v] = dmarray(value, attrs=data[v].attrs)
# Create data object for holding extracted values.
# One vector for each value w/ same units as parent object.
for v in self._var_list:
self[v] = dmarray(np.zeros(x.shape), attrs=data[v].attrs)
# Some helpers:
xAll = data[data['grid'].attrs['dims'][0]]
yAll = data[data['grid'].attrs['dims'][1]]
# Navigate the MHD quad tree, interpolating as we go.
for k in data.qtree:
# Only leafs, of course!
if not data.qtree[k].isLeaf: continue
# Find all points that lie in this block.
# If there are none, just keep going.
pts = \
(x >= data.qtree[k].lim[0]) & (x <= data.qtree[k].lim[1]) & \
(y >= data.qtree[k].lim[2]) & (y <= data.qtree[k].lim[3])
if not pts.any(): continue
locs = data.qtree[k].locs
for v in self._var_list:
self[v][pts] = interp_2d_reg(x[pts], y[pts], xAll[locs],
yAll[locs], data[v][locs])
[docs]class Stream(Extraction):
'''
A class for streamlines. Contains all of the information about
the streamline, including extracted variables.
Upon instantiation, the object will trace through the vector
field determined by the "[x/y]field" values and the Bats object
"bats".
Parameters
----------
bats : Bats
:class:`~spacepy.pybats.bats.Bats2d` object through which to trace.
xstart : float
X value of location to start the trace.
ystart : float
Y value of location to start the trace.
xfield : str
Name of variable in ``bats`` which contains X values of the field
yfield : str
Name of variable in ``bats`` which contains Y values of the field
Other Parameters
----------------
style : str
Sets line style, including colors. See :meth:`set_style` for details.
(Default 'mag')
type : str
(Default 'streamline')
method : str
Integration method. The default is Runge-Kutta 4 ('rk4') which gives
a good blend of speed and accuracy. See the test functions in
:mod:`~spacepy.pybats.trace2d` for more info. The other option is
a simple Euler's method approach ('eul'). (Default 'rk4')
extract : bool
(Default: False) Extract variables along stream trace and save within
object.
maxPoints : int
(Default : 20000) Maximum number of integration steps to take.
var_list : string or sequence of strings
(Default : 'all') List of values to extract from *dataset*.
Defaults to 'all', for all values within *bats*.
Notes
-----
.. Not really "notes" but need to keep this section from being parsed
as parameters
.. rubric:: Methods
.. autosummary::
~Stream.set_style
~Stream.treetrace
~Stream.trace
~Stream.plot
.. automethod:: set_style
.. automethod:: treetrace
.. automethod:: trace
.. automethod:: plot
'''
def __init__(self, bats, xstart, ystart, xfield, yfield, style = 'mag',
type='streamline', method='rk4', var_list='all',
extract=False, maxPoints=20000, *args, **kwargs):
# Key values:
self.xstart = xstart #X and Y starting
self.ystart = ystart #points in the field.
self.xvar = xfield #Strings that list the variables
self.yvar = yfield #that will be used for tracing.
# Descriptors:
self.open = True
self.status = 'open'
self.method = method
# Do tracing if a tracing method has been set.
self.treetrace(bats, maxPoints=maxPoints)
# Now, initialize as Extraction using self.x, self.y.
# This handles the extraction as well as converting self.x, y into
# dmarrays with the correct names.
super(Stream, self).__init__(self.x, self.y, bats,
var_list=var_list*extract)
# Place parameters into attributes:
self.attrs['start'] = [xstart, ystart]
self.attrs['trace_var'] = [xfield, yfield]
self.attrs['method'] = method
# Set style
self.set_style(style)
#def __repr__(self):
# pass
#
#def __str__(self):
# pass
[docs] def set_style(self, style):
'''
Set the line style either using a simple matplotlib-type style
string or using a preset style type. Current types include:
'mag' : treat line as a magnetic field line. Closed lines are
white, other lines are black.
'''
import re
# Here, we can set line color and style based on
# line characteristics. Right now, only one preset is
# available.
if style == 'mag':
if self.open:
self.style = 'k-'
self.color = 'k'
else:
self.style = 'w-'
self.color = 'w'
else:
self.style = style
col = re.match('([bgrcmykw])', style)
if col:
self.color=col.groups()[0]
else:
self.color='k'
[docs] def treetrace(self, bats, maxPoints=20000):
'''
Trace through the vector field using the quad tree.
'''
from numpy import array, sqrt, append
if self.method == 'euler' or self.method == 'eul':
from spacepy.pybats.trace2d import trace2d_eul as trc
elif self.method == 'rk4':
from spacepy.pybats.trace2d import trace2d_rk4 as trc
else:
raise ValueError('Tracing method {} not recognized!'.format(
self.method))
# Get name of dimensions in order.
grid = bats['grid'].attrs['dims']
# Find starting block and set starting locations.
block = bats.find_block(self.xstart, self.ystart)
xfwd = [self.xstart]
yfwd = [self.ystart]
xnow, ynow = self.xstart, self.ystart
# Trace forwards.
while(block):
# Grab indices of all values inside current block.
# Ghost cell check is for experimental testing:
if hasattr(bats.qtree[block], 'ghost'):
loc = bats.qtree[block].ghost
# Typical approach is no ghost cell info:
else:
loc = bats.qtree[block].locs
# Trace through this block.
x, y = trc(bats[self.xvar][loc], bats[self.yvar][loc],
xnow, ynow, bats[grid[0]][loc][0,:],
bats[grid[1]][loc][:,0], ds=0.01)
# Update location and block:
xnow, ynow = x[-1], y[-1]
newblock = bats.find_block(xnow,ynow)
# If we didn't leave the block, stop tracing.
# Additionally, if inside rBody, stop.
if(block == newblock) or (xnow**2+ynow**2) < bats.attrs['rbody'] * .8 :
block = False
elif newblock:
block = newblock
else:
block = False
# Append to full trace vectors.
xfwd = np.append(xfwd, x[1:])
yfwd = np.append(yfwd, y[1:])
del(x); del(y)
# It's possible to get stuck swirling around across
# a few blocks. If we spend a lot of time tracing,
# call it quits.
if xfwd.size > maxPoints: block = False
# Trace backwards. Same Procedure as above.
block = bats.find_block(self.xstart, self.ystart)
xbwd = [self.xstart]
ybwd = [self.ystart]
xnow, ynow = self.xstart, self.ystart
while(block):
if hasattr(bats.qtree[block], 'ghost'):
loc = bats.qtree[block].ghost
else:
loc = bats.qtree[block].locs
x, y = trc(bats[self.xvar][loc], bats[self.yvar][loc],
xnow, ynow, bats[grid[0]][loc][0, :],
bats[grid[1]][loc][:, 0], ds=-0.01)
xnow, ynow = x[-1], y[-1]
newblock = bats.find_block(xnow, ynow)
if(block == newblock) or (xnow**2+ynow**2) < bats.attrs['rbody'] * .8:
block = False
elif newblock:
block = newblock
else:
block = False
# Append to full trace vectors.
xbwd = np.append(x[::-1], xbwd)
ybwd = np.append(y[::-1], ybwd)
if xbwd.size > maxPoints:
block = False
# Trim duplicate points created when backwards tracing.
# There's always at least 1 duplicate.
for i in range(1, xbwd.size):
if xbwd[-(i+1)]-xfwd[0] != 0: break
# Combine foward and backward traces.
self.x = np.append(xbwd[:-i], xfwd)
self.y = np.append(ybwd[:-i], yfwd)
# If planetary run w/ body:
# 1) Check if line is closed to body.
# 2) Trim points within body.
if 'rbody' in bats.attrs:
# Radial distance:
r = sqrt(self.x**2.0 + self.y**2.0)
# Closed field line? Lobe line? Set status:
if (r[0] < bats.attrs['rbody']) and (r[-1] < bats.attrs['rbody']):
self.open = False
self.status = 'closed'
elif (r[0] > bats.attrs['rbody']) and (r[-1] < bats.attrs['rbody']):
self.open = True
self.status = 'north lobe'
elif (r[0] < bats.attrs['rbody']) and (r[-1] > bats.attrs['rbody']):
self.open = True
self.status = 'south lobe'
# Trim the fat!
limit = bats.attrs['rbody']*.8
self.x, self.y = self.x[r > limit], self.y[r > limit]
[docs] def trace(self, bats):
'''
Trace through the vector field.
'''
from numpy import array, sqrt
if self.method == 'euler':
from spacepy.pybats.trace2d import trace2d_eul as trc
elif self.method == 'rk4':
from spacepy.pybats.trace2d import trace2d_rk4 as trc
# Get name of dimensions in order.
grid = bats['grid'].attrs['dims']
# Trace forward
x1, y1 = trc(bats[self.xvar], bats[self.yvar],
self.xstart, self.ystart,
bats[grid[0]], bats[grid[1]])
# Trace backward
x2, y2 = trc(bats[self.xvar], bats[self.yvar],
self.xstart, self.ystart,
bats[grid[0]], bats[grid[1]], ds=-0.1)
# Join lines together such that point 0 is beginning of line
# and point -1 is the end (when moving parallel to line.)
self.x = array(x2[::-1].tolist() + x1[1:].tolist())
self.y = array(y2[::-1].tolist() + y1[1:].tolist())
# Check if line is closed to body.
if 'rbody' in bats.attrs:
r1 = sqrt(self.x[ 0]**2.0 + self.y[0]**2.0)
r2 = sqrt(self.x[-1]**2.0 + self.y[-1]**2.0)
if (r1 < bats.attrs['rbody']) and (r2 < bats.attrs['rbody']):
self.open = False
[docs] def plot(self, ax, *args, **kwargs):
'''
Add streamline to axes object "ax".
'''
ax.plot(self.x, self.y, self.style, *args, **kwargs)
[docs]class Bats2d(IdlFile):
'''
A child class of :class:`~pybats.IdlFile` tailored to 2D BATS-R-US output.
Calculations
------------
New values can be added via the addition of new keys. For example,
a user could add radial distance to an equatorial Bats2d object as follows:
>>> import numpy as np
>>> from spacepy.pybats import bats
>>> mhd = bats.Bats2d('z=0_example.out')
>>> mhd['rad'] = np.sqrt( mhd['x']**2 + mhd['y']**2 )
Note, however, that if the user switches the data frame in a *.outs file
to access data from a different epoch, these values will need to be
updated.
An exception to this is built-in `calc_*` methods, which perform common
MHD/fluid dynamic calculations (i.e., Alfven wave speed, vorticity, and
more.) These values are updated when the data frame is switched (see the
`switch_frame` method).
Plotting
--------
While users can employ Matplotlib to plot values, a set of built-in
methods are available to expedite plotting. These are the
`add_<plot type>` methods. These methods always have the following
keyword arguments that allow users to optionally build more complicated
plots: *target* and *loc*. The *target* kwarg tells the plotting method
where to place the plot and can either be a Matplotlib figure or axes
object. If it's an axes object, *loc* sets the subplot location using
the typical matplotlib syntax (e.g., `loc=121`). The default behavior is
to create a new figure and axes object.
This approach allows a user to over-plot contours, field lines, and
other plot artists as well as combine different subplots onto a single
figure. Continuing with our example above, let's plot the grid layout
for our file as well as equatorial pressure and flow streamlines:
>>> import matplotlib.pyplot as plt
>>> fig = plt.Figure(figsize=(8,6))
>>> mhd.add_grid_plot(target=fig, loc=121)
>>> mhd.add_contour('x','y','p', target=fig, loc=122)
>>> mhd.add_stream_scatter('ux', 'uy', target=fig, loc=122)
Useful plotting methods include the following:
| Plot Method | Description |
| ------------------ | ---------------------------------------------- |
| add_grid_plot | Create a quick-look diagram of the grid layout |
| add_contour | Create a contour plot of a given variable |
| add_pcolor | Add a p-color (no-interpolation contour) plot |
| add_stream_scatter | Scatter stream traces (any vector field) |
| add_b_magsphere | Add magnetic field lines for X-Z plane cuts |
| add_planet | Add a simple black/white planet at the origin |
| add_body | Add an inner boundary at the origin |
Extracting and Stream Tracing
-----------------------------
Extracting values via interpolation to arbitrary points and creating
stream traces through any vector field (e.g., velocity or magnetic field)
are aided via the use of the following object methods:
| Method | Description |
| ---------- | -------------------------------------------------- |
| extract | Interpolate to arbitrary points and extract values |
| get_stream | Integrate stream lines through vector fields |
Be sure to read the docstring information of :class:`~pybats.IdlFile` to
see how to handle multi-frame files (*.outs) and for a list of critical
attributes.
'''
# Init by calling IdlFile init and then building qotree, etc.
def __init__(self, filename, *args, **kwargs):
# Create quad tree object attribute:
self._qtree = None
# Read file.
IdlFile.__init__(self, filename, keep_case=False, *args, **kwargs)
# Behavior of output files changed Jan. 2017:
# Check for 'r' instead of 'rbody' in attrs.
if 'r' in self.attrs and 'rbody' not in self.attrs:
self.attrs['rbody'] = self.attrs['r']
def switch_frame(self, *args, **kwargs):
'''
For files that have more than one data frame (i.e., `*.outs` files),
load data from the *iframe*-th frame into the object replacing what is
currently loaded.
'''
# Reset our quad tree before reading our file:
self._qtree = None
# Switch frames using parent method:
super(Bats2d, self).switch_frame(*args, **kwargs)
@property
def qtree(self):
if self._qtree is None:
from numpy import array
import spacepy.pybats.qotree as qo
# Parse grid into quad tree.
if self['grid'].attrs['gtype'] != 'Regular':
xdim, ydim = self['grid'].attrs['dims'][0:2]
try:
self._qtree = qo.QTree(array([self[xdim], self[ydim]]))
except:
from traceback import print_exc
print_exc()
self._qtree = False
self.find_block = lambda: False
else:
self._qtree = False
self.find_block = lambda: False
return self._qtree
@property
def find_block(self):
if self.qtree:
return self._qtree.find_leaf
else:
return lambda: False
####################
# CALCULATIONS
####################
@calc_wrapper
def calc_temp(self, units='eV'):
'''
Calculate plasma temperature for each fluid. Number density is
calculated using *calc_ndens* if it hasn't been done so already.
Temperature is obtained via density and pressure through the simple
relationship P=nkT.
Use the units kwarg to set output units. Current choices are
KeV, eV, and K. Default is eV.
'''
from spacepy.datamodel import dmarray
units = units.lower()
# Create dictionary of unit conversions.
conv = {'ev' : 6241.50935, # nPa/cm^3 --> eV.
'kev': 6.24150935, # nPa/cm^3 --> KeV.
'k' : 72429626.47} # nPa/cm^3 --> K.
# Calculate number density if not done already.
if 'N' not in self:
self.calc_ndens()
# Find all number density variables.
for key in list(self.keys()):
# Next variable if not number density:
if key[-1] != 'N':
continue
# Next variable if no matching pressure:
if not key[:-1]+'p' in self:
continue
self[key[:-1]+'t'] = dmarray(
conv[units] * self[key[:-1]+'p']/self[key],
attrs={'units':units})
@calc_wrapper
def calc_b(self):
'''
Calculates total B-field strength using all three B components.
Retains units of components. Additionally, the unit vector
b-hat is calculated as well.
'''
from numpy import sqrt
self['b'] = sqrt(self['bx']**2.0 + self['by']**2.0 + self['bz']**2.0)
self['b'].attrs = {'units':self['bx'].attrs['units']}
self['bx_hat'] = self['bx'] / self['b']
self['by_hat'] = self['by'] / self['b']
self['bz_hat'] = self['bz'] / self['b']
self['bx_hat'].attrs = {'units':'unitless'}
self['by_hat'].attrs = {'units':'unitless'}
self['bz_hat'].attrs = {'units':'unitless'}
@calc_wrapper
def calc_j(self):
'''
Calculates total current density strength using all three J components.
Retains units of components, stores in self['j']
'''
from numpy import sqrt
self['j'] = sqrt(self['jx']**2.0 + self['jy']**2.0 + self['jz']**2.0)
self['j'].attrs = {'units':self['jx'].attrs['units']}
@calc_wrapper
def calc_E(self):
'''
Calculates the MHD electric field, -UxB. Works for default
MHD units of nT and km/s; if these units are not correct, an
exception will be raised. Stores E in mV/m.
'''
# Some quick declarations for more readable code.
ux = self['ux']; uy=self['uy']; uz=self['uz']
bx = self['bx']; by=self['by']; bz=self['bz']
# Check units. Should be nT(=Volt*s/m^2) and km/s.
if (bx.attrs['units'] != 'nT') or (ux.attrs['units'] != 'km/s'):
raise Exception('Incorrect units! Should be km/s and nT.')
# Calculate; return in millivolts per meter
self['Ex'] = -1.0*(uy*bz - uz*by) / 1000.0
self['Ey'] = -1.0*(uz*bx - ux*bz) / 1000.0
self['Ez'] = -1.0*(ux*by - uy*bx) / 1000.0
self['Ex'].attrs = {'units':'mV/m'}
self['Ey'].attrs = {'units':'mV/m'}
self['Ez'].attrs = {'units':'mV/m'}
# Total magnitude.
self['E'] = np.sqrt(self['Ex']**2+self['Ey']**2+self['Ez']**2)
@calc_wrapper
def calc_ndens(self):
'''
Calculate number densities for each fluid. Species mass is ascertained
via recognition of fluid name (e.g. OpRho is clearly oxygen). A full
list of recognized fluids/species can be found by exploring the
dictionary *mass* found in :mod:`~spacepy.pybats.bats`. Composition is
also calculated as percent of total number density.
New values are saved using the keys *speciesN* (e.g. *opN*) and
*speciesFrac* (e.g. *opFrac*).
'''
# Use shared function.
_calc_ndens(self)
@calc_wrapper
def calc_beta(self):
'''
Calculates plasma beta (ratio of plasma to magnetic pressure,
indicative of who - plasma or B-field - is "in charge" with regards
to flow.
Assumes:
-pressure in units of nPa
-B in units of nT.
Resulting value is unitless.
Values where b_total = zero are set to -1.0 in the final array.
'''
from numpy import pi
if 'b' not in self:
self.calc_b()
mu_naught = 4.0E2 * pi # Mu_0 x unit conversion (nPa->Pa, nT->T)
temp_b = self['b']**2.0
temp_b[temp_b < 1E-8] = -1.0*mu_naught*self['p'][temp_b == 0.0]
temp_beta = mu_naught*self['p']/temp_b
self['beta'] = temp_beta
self['beta'].attrs = {'units':'unitless'}
@calc_wrapper
def calc_jxb(self):
'''
Calculates the JxB force assuming:
-Units of J are uA/m2, units of B are nT.
Under these assumptions, the value returned is force density (nN/m^3).
'''
from numpy import sqrt
from spacepy.datamodel import dmarray
# Unit conversion (nT, uA, cm^-3 -> nT, A, m^-3) to nN/m^3.
conv = 1E-6
# Calculate cross product, convert units.
self['jbx'] = dmarray((self['jy']*self['bz']-self['jz']*self['by'])*conv,
{'units':'nN/m^3'})
self['jby'] = dmarray((self['jz']*self['bx']-self['jx']*self['bz'])*conv,
{'units':'nN/m^3'})
self['jbz'] = dmarray((self['jx']*self['by']-self['jy']*self['bx'])*conv,
{'units':'nN/m^3'})
self['jb'] = dmarray(sqrt(self['jbx']**2 +
self['jby']**2 +
self['jbz']**2), {'units':'nN/m^3'})
@calc_wrapper
def calc_alfven(self):
'''
Calculate the Alfven speed, B/(mu*rho)^(1/2) in km/s. This is
performed for each species and the total fluid.
The variable is saved under key "alfven" in self.data.
'''
from numpy import sqrt, pi
from spacepy.datamodel import dmarray
if 'b' not in self:
self.calc_b()
#M_naught * conversion from #/cm^3 to kg/m^3
mu_naught = 4.0E-7 * pi * 1.6726E-27 * 1.0E6
# Get all rho-like variables. Save in new list.
rho_names = []
for k in self:
if (k[-3:]) == 'rho':
rho_names.append(k)
# Calculate Alfven speed in km/s. Separate step to avoid
# changing dictionary while looping over keys.
for k in rho_names:
self[k[:-3]+'alfven'] = dmarray(self['b']*1E-12 /
sqrt(mu_naught*self[k]),
attrs={'units':'km/s'})
@calc_wrapper
def _calc_divmomen(self):
'''
Calculate the divergence of momentum, i.e.
$\rho(u \dot \nabla)u$.
This is currently exploratory.
'''
from spacepy.datamodel import dmarray
from spacepy.pybats.batsmath import d_dx, d_dy
if self.qtree == False:
raise ValueError('calc_divmomen requires a valid qtree')
# Create empty arrays to hold new values.
size = self['ux'].shape
self['divmomx'] = dmarray(np.zeros(size), {'units':'nN/m3'})
self['divmomz'] = dmarray(np.zeros(size), {'units':'nN/m3'})
# Units!
c1 = 1000. / 6371.0 # km2/Re/s2 -> m/s2
c2 = 1.6726E-21 # AMU/cm3 -> kg/m3
c3 = 1E9 # N/m3 -> nN/m3.
for k in self.qtree:
# Calculate only on leafs of quadtree.
if not self.qtree[k].isLeaf: continue
# Extract values from current leaf.
leaf = self.qtree[k]
ux = self['ux'][leaf.locs]
uz = self['uz'][leaf.locs]
self['divmomx'][leaf.locs] = ux*d_dx(ux, leaf.dx)+uz*d_dy(ux, leaf.dx)
self['divmomz'][leaf.locs] = ux*d_dx(uz, leaf.dx)+uz*d_dy(uz, leaf.dx)
# Unit conversion.
self['divmomx'] *= self['rho']*c1*c2*c3
self['divmomz'] *= self['rho']*c1*c2*c3
@calc_wrapper
def calc_vort(self, conv=1000./RE):
'''
Calculate the vorticity (curl of bulk velocity) for the direction
orthogonal to the cut plane. For example, if output file is
a cut in the equatorial plane (GSM X-Y plane), only the z-component
of the curl is calculated.
Output is saved as self['wD'], where D is the resulting dimension
(e.g., 'wz' for the z-component of vorticity in the X-Y plane).
Parameters
==========
None
Other Parameters
================
conv : float
Required unit conversion such that output units are 1/s.
Defaults to 1/RE (in km), which assumes grid is in RE and
velocity is in km/s.
'''
from spacepy.pybats.batsmath import d_dx, d_dy
if self.qtree is False:
raise ValueError('calc_vort requires a valid qtree')
# Determine which direction to calculate based on what direction
# is not present. Save appropriate derivative operators, order
# useful directions, create new variable name.
dims = self['grid'].attrs['dims']
if 'x' not in dims:
w = 'wx'
dim1, dim2 = 'z', 'y'
dx1, dx2 = d_dy, d_dx
elif 'y' not in dims:
w = 'wy'
dim1, dim2 = 'x', 'z'
dx1, dx2 = d_dy, d_dx
else:
w = 'wz'
dim1, dim2 = 'x', 'y'
dx1, dx2 = d_dy, d_dx
# Create new arrays to hold curl.
size = self['ux'].shape
self[w] = dmarray(np.zeros(size), {'units':'1/s'})
# Navigate quad tree, calculate curl at every leaf.
for k in self.qtree:
# Plot only leafs of the tree.
if not self.qtree[k].isLeaf: continue
# Get location of points and extract velocity:
leaf = self.qtree[k]
u1 = self['u'+dim1][leaf.locs]
u2 = self['u'+dim2][leaf.locs]
# Calculate curl
self[w][leaf.locs] = conv * (dx1(u1, leaf.dx) - dx2(u2, leaf.dx))
@calc_wrapper
def calc_gradP(self):
'''
Calculate the pressure gradient force.
'''
from spacepy.datamodel import dmarray
from spacepy.pybats.batsmath import d_dx, d_dy
if self.qtree is False:
raise ValueError('calc_gradP requires a valid qtree')
if 'p' not in self:
raise KeyError('Pressure not found in object')
# Create new arrays to hold pressure.
dims = self['grid'].attrs['dims']
size = self['p'].shape
self['gradP'] = dmarray(np.zeros(size), {'units':'nN/cm^3'})
for d in dims:
self['gradP_'+d] = dmarray(np.zeros(size), {'units':'nN/m^3'})
for k in self.qtree:
# Plot only leafs of the tree.
if not self.qtree[k].isLeaf: continue
# Extract leaf; place pressure into 2D array.
leaf = self.qtree[k]
z = self['p'][leaf.locs]
# Calculate derivatives; place into new dmarrays.
# Unit conversion: P in nPa => gradP in nN/m3, dx=Re
# Convert to nN/m3 by multiplying by 6378000m**-1.
# Clever indexing maps values back to "unordered" array so that
# values can be plotted like all others.
conv = 1.0/-6378000.0
self['gradP_'+dims[0]][leaf.locs] = d_dx(z, leaf.dx)*conv
self['gradP_'+dims[1]][leaf.locs] = d_dy(z, leaf.dx)*conv
# Scalar magnitude:
for d in dims:
self['gradP'] += self['gradP_'+d]**2
self['gradP'] = np.sqrt(self['gradP'])
@calc_wrapper
def calc_utotal(self):
'''
Calculate bulk velocity magnitude: $u^2 = u_X^2 + u_Y^2 + u_Z^2$.
This is done on a per-fluid basis.
'''
from numpy import sqrt
from spacepy.datamodel import dmarray
species = []
# Find all species, the variable names end in "ux".
for k in self:
if (k[-2:] == 'ux') and (k[:-2]+'u' not in self):
species.append(k[:-2])
units = self['ux'].attrs['units']
for s in species:
self[s+'u'] = dmarray(sqrt(self[s+'ux']**2 +
self[s+'uy']**2 +
self[s+'uz']**2),
attrs={'units':units})
@calc_wrapper
def _calc_Ekin(self, units='eV'):
'''
Calculate average kinetic energy per particle using
$E=\frac{1}{2}mv^2$. Note that this is not the same as energy
density. Units are $eV$.
'''
from numpy import sqrt
from spacepy.datamodel import dmarray
raise Warning("This calculation is unverified.")
conv = 0.5 * 0.0103783625 # km^2-->m^2, amu-->kg, J-->eV.
if units.lower == 'kev':
conv = conv/1000.0
species = []
# Find all species, the variable names end in "rho".
for k in self:
if (k[-3:] == 'rho') and (k[:-3]+'Ekin' not in self):
species.append(k[:-3])
for s in species:
#THIS IS WRONG HERE: 1/2mV**2? Notsomuch.
self[s+'Ekin'] = dmarray(sqrt(self[s+'ux']**2 +
self[s+'uy']**2 +
self[s+'uz']**2)
* conv * mass[s.lower()],
attrs={'units':units})
def calc_all(self, exclude=[]):
'''
Perform all variable calculations (e.g. calculations that
begin with 'calc_'). Any exceptions raised by functions that
could not be peformed (typicaly from missing variables) are
discarded.
'''
for command in dir(self):
if (command[0:5] == 'calc_') and (command != 'calc_all') and (command not in exclude):
try:
eval('self.'+command+'()')
except AttributeError:
raise Warning('Did not perform {0}: {1}'.format(
command, sys.exc_info()[0]))
#####################
# Other calculations
#####################
def gradP_regular(self, cellsize=None, dim1range=-1, dim2range=-1):
'''
Calculate pressure gradient on a regular grid.
Note that if the Bats2d object is not on a regular grid, one of
two things will happen. If kwarg cellsize is set, the value of
cellsize will be used to call self.regrid and the object will
be regridded using a cellsize of cellsize. Kwargs dim1range and
dim2range can be used in the same way they are used in self.regrid
to restrict the regridding to a smaller domain. If cellsize is
not set and the object is on an irregular grid, an exception is
raised.
The gradient is calculated using numpy.gradient. The output units
are force density (N/m^3). Three variables are added to self.data:
gradp(dim1), gradp(dim2), gradp. For example, if the object is an
equatorial cut, the variables gradpx, gradpy, and gradp would be
added representing the gradient in each direction and then the
magnitude of the vector.
'''
from numpy import gradient, sqrt
if self.gridtype != 'Regular':
if not cellsize:
raise ValueError('Grid must be regular or ' +
'cellsize must be given.')
self.regrid(cellsize, dim1range=dim1range, dim2range=dim2range)
# Order our dimensions alphabetically.
newvars = []
for key in sorted(self.grid.keys()):
newvars.append('gradp'+key)
dx=self.resolution*6378000.0 # RE to meters
p =self['p']*10E-9 # nPa to Pa
self[newvars[0]], self[newvars[1]]=gradient(p, dx, dx)
self['gradp']=sqrt(self[newvars[0]]**2.0 + self[newvars[1]]**2.0)
def cfl(self,dt,xcoord='x',ycoord='z'):
"""
Calculate the CFL number in each cell given time step dt.
Result is stored in self['cfl'].
"""
try:
U=self['u']
except KeyError:
self.calc_utotal()
cfl=np.zeros(self['u'].shape)
for key in list(self.qtree.keys()):
child=self.qtree[key]
if not child.isLeaf: continue
pts=np.where(np.logical_and.reduce((self[xcoord]>=child.lim[0],
self[xcoord]<=child.lim[1],
self[ycoord]>=child.lim[2],
self[ycoord]<=child.lim[3]))
)
cfl[pts]=self['u'][pts]*dt/child.dx
self['cfl']=dmarray(cfl,attrs={'units':''})
def vth(self,m_avg=3.1):
"""
Calculate the thermal velocity. m_avg denotes the average ion
mass in AMU.
Result is stored in self['vth'].
"""
m_avg_kg=m_avg*1.6276e-27
ndensity=self['rho']/m_avg*1e6
self['vth']=dmarray(np.sqrt(self['p']*1e-9/ndensity/(m_avg_kg))
/1000,attrs={'units':'km/s'})
def gyroradius(self,velocities=('u','vth'),m_avg=3.1):
"""
Calculate the ion gyroradius in each cell.
velocities to use in calculating the gyroradius are listed by name in
the sequence argument velocities. If more than one variable is given,
they are summed in quadrature.
m_avg denotes the average ion mass in AMU.
Result is stored in self['gyroradius']
"""
if 'vth' in velocities:
try:
self['vth']
except KeyError:
self.vth()
if 'u' in velocities:
try:
U=self['u']
except KeyError:
self.calc_utotal()
velocities_squared_sum=self[velocities[0]]**2
for vname in velocities[1:]:
velocities_squared_sum+=self[vname]**2
v=np.sqrt(velocities_squared_sum)*1000
if 'b' not in self: self.calc_b()
B=self['b']*1e-9
m_avg_kg=m_avg*1.6276e-27
q=1.6022e-19
self['gyroradius']=dmarray(m_avg_kg*v/(q*B)/6378000,
attrs={'units':'Re'})
def plasma_freq(self,m_avg=3.1):
"""
Calculate the ion plasma frequency.
m_avg denotes the average ion mass in AMU.
Result is stored in self['plasm_freq'].
"""
from numpy import sqrt, pi
m_avg_kg=m_avg*1.6276e-27
ndensity=self['rho']/m_avg*1e6
q=1.6022e-19
self['plasma_freq']=dmarray(sqrt(4*pi*ndensity*q**2/m_avg_kg),
attrs={'units':'rad/s'})
def inertial_length(self,m_avg=3.1):
"""
Calculate the ion inertial length.
m_avg denotes the average ion mass in AMU.
Result is stored in self['inertial_length'].
"""
try:
self['plasma_freq']
except KeyError:
self.plasma_freq(m_avg=m_avg)
if 'alfven' not in self: self.calc_alfven()
self['inertial_length']=dmarray(self['alfven']/self['plasma_freq']
/6378000, attrs={'units':'Re'})
def regrid(self, cellsize=1.0, dim1range=-1, dim2range=-1, debug=False):
'''
Re-bin data to regular grid of spacing cellsize. Action is
performed on all data entries in the bats2d object.
'''
from matplotlib.mlab import griddata
if self['grid'].attrs['gtype'] == 'Regular': return
# Order our dimensions alphabetically.
dims = self['grid'].attrs['dims']
if debug: print("Ordered dimensions: ", dims)
# Check to see if dimranges are 2-element lists.
# If not, either set defaults or raise exceptions.
if dim1range == -1:
dim1range = [self[dims[0]].min(),self[dims[0]].max()]
else:
if isinstance(dim1range, ( type(()), type([]) ) ):
if len(dim1range) != 2:
raise ValueError('dim1range must have two elements!')
else: raise TypeError('dim1range must be a tuple or list!')
if dim2range == -1:
dim2range = [self[dims[1]].min(),self[dims[1]].max()]
else:
if isinstance(dim2range, ( type(()), type([]) ) ):
if len(dim2range) != 2:
raise ValueError('dim2range must have two elements!')
else: raise TypeError('dim2range must be a tuple or list!')
if debug:
print('%s range = %f, %f' % (dims[0], dim1range[0], dim1range[1]))
print('%s range = %f, %f' % (dims[1], dim2range[0], dim2range[1]))
# Now, Regrid.
grid1 = np.arange(dim1range[0], dim1range[1]+cellsize, cellsize)
grid2 = np.arange(dim2range[0], dim2range[1]+cellsize, cellsize)
for key in self:
# Skip grid-type entries.
if key in (self['grid'].attrs['dims']+['grid']): continue
self[key] = griddata(self[dims[0]], self[dims[1]],
self[key], grid1, grid2)
# Change grid, gridtype, gridsize, and npoints to match new layout.
self['grid'].attrs['gtype'] = 'Regular'
self['grid'].attrs['npoints'] = len(grid1) * len(grid2)
self['grid'].attrs['resolution'] = cellsize
self[dims[0]] = grid1; self[dims[1]] = grid2
#############################
# EXTRACTION/INTERPOLATION
#############################
def extract(self, x, y, **kwargs):
'''
For x, y of a 1D curve, extract values along that curve
and return slice as a new :class:`~spacepy.pybats.bats.Extraction`
object. Valid keyword arguments are the same as for
:class:`~spacepy.pybats.bats.Extraction`.
'''
return Extraction(x, y, self, **kwargs)
######################
# TRACING TOOLS
######################
def get_stream(self, x, y, xvar, yvar, method='rk4', style='mag',
maxPoints=40000, extract=False):
'''
Trace a 2D streamline through the domain, returning a Stream
object to the caller.
x and y set the starting point for the tracing.
xvar and yvar are string keys to self.data that define the
vector field through which this function traces.
The method kwarg sets the numerical method to use for the
tracing. Default is Runge-Kutta 4 (rk4).
'''
startvals = [x, y]
dims = self['grid'].attrs['dims']
for v, d in zip(startvals, dims):
if v < self[d].min() or v > self[d].max():
raise ValueError('Start value {} out of range for variable {}.'
.format(v, d))
stream = Stream(self, x, y, xvar, yvar, style=style,
maxPoints=maxPoints, method=method, extract=extract)
return stream
######################
# VISUALIZATION TOOLS
######################
def add_grid_plot(self, target=None, loc=111, do_label=True,
show_nums=False, show_borders=True, cmap='jet_r',
title='BATS-R-US Grid Layout'):
'''
Create a plot of the grid resolution by coloring regions of constant
resolution. Kwarg "target" specifies where to place the plot and can
be a figure, an axis, or None. If target is a figure, a new subplot
is created. The subplot location is set by the kwarg "loc", which
defaults to 111. If target is an axis, the plot is placed into that
axis object. If target is None, a new figure and axis are created
and used to display the plot.
Resolution labels can be disabled by setting kwarg do_label to False.
Plot title is set using the 'title' kwarg, defaults to 'BATS-R-US
Grid Layout'.
Note that if target is not an axis object, the axis will automatically
flip the direction of positive X GSM and turn on equal aspect ratio.
In other words, if you want a customized axis, it's best to create
one yourself.
Figure and axis, even if none given, are returned.
Parameters
----------
target : Matplotlib Figure or Axes object
Set plot destination. Defaults to new figure.
loc : 3-digit integer
Set subplot location. Defaults to 111.
do_label : boolean
Adds resolution legend to righthand margin. Defaults to True.
show_nums : boolean
Adds quadtree values to each region. Useful for debugging.
Defaults to False.
show_borders : True
Show black borders around each quadtree leaf. Defaults to True.
cmap : string
Sets the color map used to color each region. Must be a Matplotlib
named colormap. Defaults to 'jet_r'.
title : string
Sets the title at the top of the plot. Defaults to
'BATS-R-US Grid Layout'.
'''
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.ticker import MultipleLocator
from numpy import linspace
if self['grid'].attrs['gtype'] == 'Regular':
raise ValueError('Function not compatable with regular grids')
# Get dimensions over which we shall plot.
xdim, ydim = self['grid'].attrs['dims'][0:2]
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,10), loc=loc)
ax.set_aspect('equal')
# Set plot range based on quadtree.
ax.set_xlim([self.qtree[1].lim[0],self.qtree[1].lim[1]])
ax.set_ylim([self.qtree[1].lim[2],self.qtree[1].lim[3]])
# Plot.
if(show_borders):
for key in list(self.qtree.keys()):
self.qtree[key].plotbox(ax)
self.qtree.plot_res(ax, tag_leafs=show_nums, do_label=do_label,
cmap=cmap)
# Customize plot.
ax.set_xlabel('GSM %s' % xdim.upper())
ax.set_ylabel('GSM %s' % ydim.upper())
ax.set_title(title)
#if xdim=='x':
# ax.invert_xaxis()
#if ydim=='y':
# ax.invert_yaxis()
self.add_body(ax)
return fig, ax
def add_stream_scatter(self, xcomp, ycomp, nlines=100, target=None, loc=111,
method='rk4', xlim=None, ylim=None, narrow=0,
arrsize=12, arrstyle='->', start_points=None,
**kwargs):
'''
Add a set of stream traces to a figure or axes that are distributed
evenly but randomly throughout the plot domain.
Lines will be seeded randomly over a given spatial range given by
*xlim* and *ylim* **OR** the range of the axes (if *target* is set to
a non-empty axes object) **OR** over the entire object domain (in that
order).
Extra keyword args are handed to matplotlib's LineCollection object:
:class:`matplotlib.collections.LineCollection`.
Parameters
==========
xcomp : string
The first component of the vector field to trace (e.g., 'bx').
ycomp : string
The second component of the vector field to trace (e.g., 'bz').
Other Parameters
================
target : Matplotlib Figure or Axes object
Set plot destination. Defaults to new figure.
loc : 3-digit integer
Set subplot location. Defaults to 111.
xlim : Two-element list/tuple
Set the range in 1st dimension over which lines will be seeded.
ylim : Two-element list/tuple
Set the range in 2nd dimension over which lines will be seeded.
nlines : int
Number of stream lines to create; default is 100.
start_points : nlinesx2 array
Set start_points to define starting location of traces instead of
using random points. This is useful for creating timeseries of
plots.
narrow : int
Add "n" arrows to each line to indicate direction. Default is
zero, or no lines. If narrow=1, arrows will be placed at
*start_points*.
arrstyle : string
Set the arrow style in the same manner as Matplotlib's
annotate function. Default is '->'.
arrsize : int
Set the size, in points, of each directional arrow. Default is 12.
Returns:
========
fig : matplotlib Figure object
ax : matplotlib Axes object
collect : matplotlib Collection object of trace results
start_points : nlines x 2 array of line starting points
'''
from numpy import array
from numpy.random import sample
from matplotlib.collections import LineCollection
from spacepy.plot import add_arrows
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,10), loc=loc)
# Try to determine the order of the dimensions used:
dims = self['grid'].attrs['dims'] # Default to standard order.
flip=False
for letters in zip(xcomp, ycomp): # Check for reverse order.
if letters == dims[::-1]: flip=True # flip it!
# Set if using the axes limits is a viable option for setting limits
# of region over which to seed lines:
use_ax_lims = False
if type(target) == type(ax):
use_ax_lims = bool(ax.artists) and \
not(ax.get_xlim() == ax.get_ylim() == (0,1))
# Set range over which to place lines. Use keyword values if provided
# OR subset of axes ranges that fit in domain (if axes are reasonable)
if xlim == None:
xlim = [self[dims[0]].min(), self[dims[0]].max()]
if use_ax_lims:
axlim = ax.get_xlim()
xlim = [max(xlim[0], axlim[0]), min(xlim[1], axlim[1])]
if ylim == None:
ylim = [self[dims[1]].min(), self[dims[1]].max()]
if use_ax_lims:
axlim = ax.get_ylim()
ylim = [max(ylim[0], axlim[0]), min(ylim[1], axlim[1])]
# If initial source points not given, create a random set:
if not start_points:
# Get random points.
start_points = sample( [nlines, 2] )
# Scale to limits:
start_points[:,0] = (xlim[1]-xlim[0])*start_points[:,0]+xlim[0]
start_points[:,1] = (ylim[1]-ylim[0])*start_points[:,1]+ylim[0]
else:
nlines = start_points.shape[-1]
# Extract stream traces, organize x and y coords:
lines = []
for xstart, ystart in start_points:
# Some index errors crop up from time to time.
# While a better solution is dug up, we use a try
# block for the time being.
try:
stream = self.get_stream(xstart, ystart, xcomp, ycomp,
method=method)
except IndexError:
continue
lines.append(array([stream.x, stream.y][::1-2*flip]).transpose())
# Create line collection & plot.
collect = LineCollection(lines, **kwargs)
ax.add_collection(collect)
# Set the plot limits to match
ax.set_xlim(xlim)
ax.set_ylim(ylim)
# Add arrows if requested. If one point per line, use starting points.
# Otherwise, distribute arrows along lines.
if narrow == 1:
add_arrows(collect, n=narrow, size=arrsize,
positions=start_points, style=arrstyle)
elif narrow > 1:
add_arrows(collect, n=narrow, size=arrsize, style=arrstyle)
return fig, ax, collect, start_points
def find_earth_lastclosed(self, tol=np.pi/360., method='rk4',
max_iter=100, debug=False):
'''
For Y=0 cuts, attempt to locate the last-closed magnetic field line
for both day- and night-sides. This is done using a bisection
approach to precisely locate the transition between open and closed
geometries. The method stops once this transition is found within
a latitudinal tolerance of *tol*, which defaults to $\pi/360.$, or
one-half degree. The tracing *method* can be set via keyword and
defaults to 'rk4' (4th order Runge Kutta, see
:class:`~spacepy.pybats.bats.Stream` for more information).
The maximum number of iterations the algorithm will take is set
by *max_iter*, which defaults to 100. Latitudinal footprints of the
last closed field lines at the inner boundary (not the ionosphere!)
are also returned.
This method returns 5 objects:
* The dipole tilt in radians
* A tuple of the northern/southern hemisphere polar angle of
footpoints for the dayside last-closed field line.
* A tuple of the northern/southern hemisphere polar angle of
footpoints for the nightside last-closed field line.
* The dayside last-closed field line as a
:class:`~spacepy.pybats.bats.Stream` object.
* The nightside last-closed field line as a
:class:`~spacepy.pybats.bats.Stream` object.
In each case, the angle is defined as elevation from the positive
x-axis, in radians.
'''
import matplotlib.pyplot as plt
from numpy import (cos, sin, pi, arctan, arctan2, sqrt)
# Get the dipole tilt by tracing a field line near the inner
# boundary. Find the max radial distance; tilt angle == angle off
# equator of point of min R (~=max |B|).
x_small = self.attrs['rbody']*-1.2 # check nightside.
stream = self.get_stream(x_small, 0, 'bx', 'bz', method=method)
r = stream.x**2 + stream.y**2
loc = r==r.max()
tilt = arctan(stream.y[loc]/stream.x[loc])[0]
if debug:
print('Dipole is tilted {} degress above the z=0 plane.'.format(
tilt*180./pi))
# Dayside- start by tracing from plane of min |B| and perp. to that:
R = self.attrs['rbody']*1.15
s1 = self.get_stream(R*cos(tilt), R*sin(tilt), 'bx','bz', method=method)
# Get initial angle and step.
theta = tilt
dTheta=np.pi/4. # Initially, search 90 degrees.
nIter = 0
while (dTheta>tol)or(s1.open):
nIter += 1
# Are we closed or open? Day or nightside?
closed = not(s1.open) # open or closed?
isNig = s1.x.mean()<0 # line on day or night side?
isDay = not isNig
# Adjust the angle towards the open-closed boundary.
theta += (closed and isDay)*dTheta # adjust nightwards.
theta -= (s1.open or isNig)*dTheta # adjust daywards.
# Trace at the new theta to further restrict angular range:
s1 = self.get_stream(R*cos(theta), R*sin(theta), 'bx', 'bz',
method=method)
# Reduce angular step:
dTheta /= 2.
if nIter>max_iter:
if debug: print('Did not converge before reaching max_iter')
break
# Possible to land on open or nightside line.
# If this happens, inch back to dayside.
isNig = s1.x.mean()<0
while (s1.open or isNig):
theta-=tol/2 # inch daywards.
s1 = self.get_stream(R*cos(theta), R*sin(theta), 'bx', 'bz',
method=method)
isNig = s1.x.mean()<0
# Use last line to get southern hemisphere theta:
npts = int(s1.x.size/2)
r = sqrt(s1.x**2+s1.y**2) # Distance from origin.
# This loc finds the point(s) nearest to Rbody.
loc = np.abs(r-self.attrs['rbody'])==np.min(np.abs(
r[:npts]-self.attrs['rbody'])) #point closest to IB.
xSouth, ySouth = s1.x[loc], s1.y[loc]
# "+ 0" syntax is to quick-copy object.
theta_day = [theta+0, 2*np.pi+arctan(ySouth/xSouth)[0]+0]
day = s1
# Nightside: Use more points in tracing (lines are long!)
theta+=tol/2.0 # Nudge nightwards.
# Set dTheta to half way between equator and dayside last-closed:
dTheta=(pi+tilt-theta)/2.
nIter = 0
while (dTheta>tol)or(s1.open):
nIter += 1
s1 = self.get_stream(R*cos(theta),R*sin(theta), 'bx','bz',
method=method, maxPoints=1E6)
# Closed? Nightside?
closed = not(s1.open)
isNig = s1.x.mean()<0
isDay = not isNig
theta -= (closed and isNig) *dTheta # closed? move poleward.
theta += (s1.open or isDay) *dTheta # open? move equatorward.
# Don't cross over into dayside territory.
if theta < theta_day[0]:
theta = theta_day[0]+tol
s1 = self.get_stream(R*cos(theta),R*sin(theta), 'bx','bz',
method=method, maxPoints=1E6)
if debug: print('No open flux over polar cap.')
break
dTheta /= 2.
if nIter>max_iter:
if debug: print('Did not converge before reaching max_iter')
break
# Use last line to get southern hemisphere theta:
npts = int(s1.x.size/2) # Similar to above for dayside.
r = sqrt(s1.x**2+s1.y**2)
loc = np.abs(r-self.attrs['rbody'])==np.min(np.abs(
r[:npts]-self.attrs['rbody']))
xSouth, ySouth = s1.x[loc], s1.y[loc]
theta_night = [theta+0, 2*np.pi+arctan2(ySouth,xSouth)[0]+0]
night = s1
#plt.plot(s1.x, s1.y, 'r-')
return tilt, theta_day, theta_night, day, night
def add_b_magsphere(self, target=None, loc=111, style='mag',
DoLast=True, DoOpen=True,
compX='bx',compY='bz', narrow=0, arrsize=12,
method='rk4', tol=np.pi/720., DoClosed=True,
colors=None, linestyles=None,
nOpen=5, nClosed=15, arrstyle='->', **kwargs):
'''
Create an array of field lines closed to the central body in the
domain. Add these lines to Matplotlib target object *target*.
If no *target* is specified, a new figure and axis are created.
Note that this should currently only be used for GSM y=0 cuts
of the magnetosphere.
A tuple containing the figure, axes, and LineCollection object
is returned.
Basic styling (color and linestyle) can be handled with the
*style*, *colors*, and *linestyles* kwargs. *style* can accept
style names as defined in :class:`~spacepy.pybats.bats.Stream`, which
colors and styles lines based on characteristics (e.g., open, closed).
The default is 'mag', which colors open lines black and closed lines
white. Alternatively, this kwarg works in a similar manner as
it does in :function:`~matplotlib.pyplot.plot`,
i.e., a string code such as "b-" (a solid blue line) or 'r:' (a
dotted red line), etc. Both *colors* and *linestyles* work much
as they do for :class:`~matplotlib.collections.LineCollection`, but
only a single value (not a list or tuple) should be provided.
*colors* can be a CSS4 color name, an RGB tuple, or a string hex code.
*linestyles* can be the name of the style (e.g., "dashed") or a
shortcut compatable with the *style* kwarg (e.g., "--"). See the
documentation for the associated Matplotlib classes & functions to
see all options. Note that *linestyles* and *colors* override
*style*.
If the styling kwargs are used, they will set the colors for all
lines except last-closed boundaries. Users may control groups
individually using multiple calls and plotting one group at a time.
Note that *colors* and *linestyles* kwargs will override *style*;
*colors* allows for more flexibility concerning color choice.
Algorithm: This method, unlike its predecessor, starts by finding
the last closed field lines via
:func:`~spacepy.pybats.bats.Bats2d.find_earth_lastclosed`. It then
fills the regions between the open and closed regions. Currently, it
does not treat purely IMF field lines.
========== ===========================================================
Kwarg Description
========== ===========================================================
target The figure or axes to place the resulting lines.
style The color coding system for field lines. Defaults to 'mag'.
See :class:`spacepy.pybats.bats.Stream`. Because lines are
added as a :class:`~matplotlib.collections.LineCollection`,
only certain styles are allowed (i.e., line styles only,
no marker styles).
loc The location of the subplot on which to place the lines.
DoLast Plot last-closed lines as red lines. Defaults to **True**.
DoOpen Plot open field lines. Defaults to **True**.
DoClosed Plot closed field lines. Defaults to **True**.
nOpen Number of closed field lines to trace per hemisphere.
Defaults to 5.
nClosed Number of open field lines to trace per hemisphere.
Defaults to 15.
narrow Add "n" arrows to each line to indicate direction.
Default is zero, or no arrows.
arrstyle Set the arrow style in the same manner as Matplotlib's
annotate function. Default is '->'.
arrsize Set the size, in points, of each directional arrow.
Default is 12.
method The tracing method; defaults to 'rk4'. See
:class:`spacepy.pybats.bats.Stream`.
tol Tolerance for finding open-closed boundary; see
:func:`~spacepy.pybats.bats.Bats2d.find_earth_lastclosed`.
compX Name of x-variable through which to trace, defaults to 'bx'.
compY Name of y-variable through which to trace, defaults to 'bz'.
colors Matplotlib-compatable color name (single) to apply to lines.
linestyles A single line style indicator, defaults to '-';
see :class:`~matplotlib.collections.LineCollection` for
possible options.
========== ===========================================================
Extra kwargs are passed to Matplotlib's LineCollection class as
described above.
Returns
=======
fig : matplotlib Figure object
ax : matplotlib Axes object
collect : matplotlib Collection object of trace results
Examples
========
>>> import matplotlib.pyplot as plt
>>> from spacepy.pybats import bats
>>> # Open a 2D slice, add a pressure contour.
>>> # Example file in spacepy/tests/data/pybats_test/:
>>> mhd = bats.Bats2d('./y0_binary.out')
>>> mhd.add_contour('x','z','p')
>>> # Add field lines using default styling:
>>> mhd.add_b_magsphere(target=plt.gca())
>>> # Add a subset of lines using custom styling:
>>> mhd.add_b_magsphere(target=plt.gca(), DoLast=False, DoOpen=False, style='g--')
'''
import re
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from numpy import (arctan, cos, sin, where, pi, log,
arange, sqrt, linspace, array)
from spacepy.plot import add_arrows
# Set ax and fig based on given target.
adj_lims = not(target) # If no target set, adjust axes limits.
fig, ax = set_target(target, figsize=(10,10), loc=111)
self.add_body(ax)
# Lines, colors, and styles:
lines = []
cols = []
# Try to get line style from "style" string.
# Default to regular line if not successful.
if not linestyles:
lstyle = re.sub('\w','',style)
if not lstyle:
linestyles = '-'
else:
linestyles=lstyle
# Start by finding open/closed boundary.
tilt, thetaD, thetaN, last1, last2 = self.find_earth_lastclosed(
method=method, tol=tol)
# Useful parameters for the following traces:
R = self.attrs['rbody']
dTheta = 1.5*np.pi/180.
dThetaN = .05*np.abs(thetaN[0]-thetaD[0])
dThetaS = .05*np.abs(thetaN[1]-thetaD[1])
## Do closed field lines ##
if DoClosed:
for tDay, tNit in zip(
linspace(0, thetaD[0]-dTheta, nClosed),
linspace(np.pi, thetaN[1]-dTheta, nClosed)):
x, y = R*cos(tDay), R*sin(tDay)
sD = self.get_stream(x,y,compX,compY,method=method,
style=style)
x, y = R*cos(tNit), R*sin(tNit)
sN = self.get_stream(x,y,compX,compY,method=method,
maxPoints=1E6,style=style)
# Append to lines, colors.
lines.append(array([sD.x, sD.y]).transpose())
lines.append(array([sN.x, sN.y]).transpose())
cols.append(sD.color)
cols.append(sN.color)
## Do open field lines ##
if DoOpen:
for tNorth, tSouth in zip(
linspace(thetaD[0]+dThetaN, thetaN[0]-dThetaN, nOpen),
linspace(thetaN[1]+dThetaS, thetaD[1]-dThetaS, nOpen)):
x, y = R*cos(tNorth), R*sin(tNorth)
sD = self.get_stream(x,y,compX,compY,method=method,
style=style)
x, y = R*cos(tSouth), R*sin(tSouth)
sN = self.get_stream(x,y,compX,compY,method=method,
style=style)
# Append to lines, colors.
lines.append(array([sD.x, sD.y]).transpose())
lines.append(array([sN.x, sN.y]).transpose())
cols.append(sD.color)
cols.append(sN.color)
## Finalize Collection ##
# If colors is given, replace what is given from
# individual lines. Keep the list-approach, however.
if colors: cols = [colors]*len(cols)
# Add last-closed field lines at end so they are plotted "on top".
if DoLast:
lines+=[array([last1.x,last1.y]).transpose(),
array([last2.x,last2.y]).transpose()]
cols+=2*['r']
# Create line collection & plot.
collect = LineCollection(lines, colors=cols, linestyles=linestyles,
**kwargs)
ax.add_collection(collect)
# Add lines if required:
if narrow>0:
add_arrows(collect, n=narrow, size=arrsize, style=arrstyle)
# On fresh axes, adjust limits from default ([0,1]):
if adj_lims:
# Set defaults:
xlim, ylim = [0,1], [0,1]
# Get x,y locations along each line:
points = [path.vertices for path in collect.get_paths()]
# Get max/min from each line, update lims:
for p in points:
xlim = min(xlim[0], p.min(0)[0]), max(xlim[1], p.max(0)[0])
ylim = min(ylim[0], p.min(0)[1]), max(ylim[1], p.max(0)[1])
# Convert to arrays for element arithmatic:
xlim, ylim = np.array(xlim), np.array(ylim)
# Add a buffer:
dX, dY = min(5,xlim[1]-xlim[0]), min(5,ylim[1]-ylim[0])
xlim+=(-dX, dX)
ylim+=(-dY, dY)
# Set new axes limits:
ax.set_xlim(xlim)
ax.set_ylim(ylim)
return fig, ax, collect
def add_b_magsphere_new(self, *args, **kwargs):
'''
This method is included to ease the transistion from the legacy
version of *add_b_magsphere* to *add_b_magsphere_new*, whose suffix
has been dropped now that it is the standard method.
'''
import warnings
print('ATTENTION: add_b_magsphere_new is now simply add_b_magsphere')
warnings.warn('add_b_magsphere_new is a candidate for removal',
category=DeprecationWarning)
return add_b_magsphere(*args, **kwargs)
def add_b_magsphere_legacy(self, target=None, loc=111, style='mag',
DoImf=False, DoOpen=False, DoTail=False,
DoDay=True, method='rk4', **kwargs):
'''
This object method is considered LEGACY and is a candidate for
removal. An updated algorithm using the same name (add_b_magsphere)
is the replacement.
Create an array of field lines closed to the central body in the
domain. Add these lines to Matplotlib target object *target*.
If no *target* is specified, a new figure and axis are created.
A tuple containing the figure, axes, and LineCollection object
is returned. Any additional kwargs are handed to the LineCollection
object.
Note that this should currently only be used for GSM y=0 cuts
of the magnetosphere.
Method:
First, the title angle is approximated by tracing a dipole-like
field line and finding the point of maximum radial distance on
that line. This is used as the magnetic equator. From this
equator, many lines are traced starting at the central body
radius. More lines are grouped together at higher magnetic
latitude to give good coverage at larger L-shells. Once an
open field line is detected, no more lines are traced. This
is done on the day and night side of the central body.
Because this method rarely captures the position of the night
side x-line, more field lines are traced by marching radially
outwards away from the furthest point from the last traced and
closed field line. This is repeated until open lines are found.
'''
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from numpy import (arctan, cos, sin, where, pi, log,
arange, sqrt, linspace, array)
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,10), loc=loc)
lines = []
colors= []
# Approximate the dipole tilt of the central body.
stream = self.get_stream(3.0, 0, 'bx', 'bz', method=method)
r = stream.x**2 + stream.y**2
loc, = where(r==r.max())
tilt = arctan(stream.y[loc[0]]/stream.x[loc[0]])
# Initial values:
daymax = tilt + pi/2.0
nightmax = tilt + 3.0*pi/2.0
# Day side:
n = arange(25)+1.0
angle = tilt + 5.0*pi*log(n)/(12.0*log(n[-1]))
for theta in angle:
x = self.attrs['rbody'] * cos(theta)
y = self.attrs['rbody'] * sin(theta)
stream = self.get_stream(x,y,'bx','bz', style=style, method=method)
if (stream.y[0] > self.attrs['rbody']) or (stream.style[0]=='k'):
daymax = theta
break
savestream = stream
if DoDay:
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
# Add IMF field lines.
if DoImf:
stream = savestream
r = sqrt(stream.x**2 + stream.y**2)
loc, = where(r==r.max())
x_mp = stream.x[loc[0]]+0.15
y_mp = stream.y[loc[0]]
delx = 2.0
for i, x in enumerate(arange(x_mp, 15.0, delx)):
# From dayside x-line out and up:
y =y_mp-x_mp+x
stream = self.get_stream(x, y, 'bx', 'bz', style=style,
method=method)
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
# From top of magnetosphere down:
y =x_mp+15.0-x+delx/3.0
stream = self.get_stream(x-delx/3.0, y, 'bx', 'bz',
method=method, style=style)
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
# From bottom of mag'sphere down:
y =x_mp-10.0-x+2.0*delx/3.0
stream = self.get_stream(x-2.0*delx/3.0, y, 'bx',
'bz', style=style, method=method)
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
# Night side:
angle = pi + tilt + pi*log(n)/(2.5*log(n[-1]))
for theta in angle:
x = self.attrs['rbody'] * cos(theta)
y = self.attrs['rbody'] * sin(theta)
stream = self.get_stream(x,y,'bx','bz', style=style, method=method)
if stream.open:
nightmax = theta
break
savestream = stream
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
# March down tail.
stream = savestream
r = sqrt(stream.x**2 + stream.y**2)
loc, = where(r==r.max())
x1 = stream.x[loc[0]]
y1 = stream.y[loc[0]]
x = x1
y = y1
while (x-1.5)>self['x'].min():
#print "Closed extension at ", x-1.5, y
#ax.plot(x-1.5, y, 'g^', ms=10)
stream = self.get_stream(x-1.5, y, 'bx', 'bz', style=style,
method=method)
r = sqrt(stream.x**2 + stream.y**2)
if stream.open:
break
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
loc, = where(r==r.max())
x = stream.x[loc[0]]
y = stream.y[loc[0]]
if x1 == x:
stream = self.get_stream(x1+1.0, y1, 'bx', 'bz', method=method)
r = sqrt(stream.x**2 + stream.y**2)
loc, = where(r==r.max())
x1 = stream.x[loc[0]]
y1 = stream.y[loc[0]]
## Add more along neutral sheet.
if DoTail:
m = (y-y1)/(x-x1)
#print "Slope = ", m
#print "From y, y1 = ", y, y1
#print "and x, x1 = ", x, x1
#ax.plot([x,x1],[y,y1], 'wo', ms=10)
xmore = arange(x, -100, -3.0)
ymore = m*(xmore-x)+y
for x, y in zip(xmore[1:], ymore[1:]):
stream = self.get_stream(x, y, 'bx', 'bz', style=style,
method=method)
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
#ax.plot(xmore, ymore, 'r+')
# Add open field lines.
if DoOpen:
for theta in linspace(daymax,0.99*(2.0*(pi+tilt))-nightmax, 15):
x = self.attrs['rbody'] * cos(theta)
y = self.attrs['rbody'] * sin(theta)
stream = self.get_stream(x,y,'bx','bz', method=method)
if stream.open:
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
x = self.attrs['rbody'] * cos(theta+pi)
y = self.attrs['rbody'] * sin(theta+pi)
stream = self.get_stream(x,y,'bx','bz', method=method)
if stream.open:
lines.append(array([stream.x, stream.y]).transpose())
colors.append(stream.style[0])
if 'colors' in kwargs:
colors=kwargs['colors']
kwargs.pop('colors')
collect = LineCollection(lines, colors=colors, **kwargs)
ax.add_collection(collect)
return fig, ax, collect
def add_planet(self, ax=None, rad=1.0, ang=0.0, **extra_kwargs):
'''
Creates a circle of radius=self.attrs['rbody'] and returns the
MatPlotLib Ellipse patch object for plotting. If an axis is specified
using the "ax" keyword, the patch is added to the plot.
Unlike the add_body method, the circle is colored half white (dayside)
and half black (nightside) to coincide with the direction of the
sun. Additionally, because the size of the planet is not intrinsically
known to the MHD file, the kwarg "rad", defaulting to 1.0, sets the
size of the planet.
Extra keywords are handed to the Ellipse generator function.
'''
from matplotlib.patches import Circle, Wedge
if 'rbody' not in self.attrs:
raise KeyError('rbody not found in self.attrs!')
body = Circle((0,0), rad, fc='w', zorder=1000, **extra_kwargs)
arch = Wedge((0,0), rad, 90.+ang, -90.+ang, fc='k',
zorder=1001, **extra_kwargs)
if ax != None:
ax.add_artist(body)
ax.add_artist(arch)
return body, arch
def add_body(self, ax=None, facecolor='lightgrey', DoPlanet=True, ang=0.0,
**extra_kwargs):
'''
Creates a circle of radius=self.attrs['rbody'] and returns the
MatPlotLib Ellipse patch object for plotting. If an axis is specified
using the "ax" keyword, the patch is added to the plot.
Default color is light grey; extra keywords are handed to the Ellipse
generator function.
Because the body is rarely the size of the planet at the center of
the modeling domain, add_planet is automatically called. This can
be negated by using the DoPlanet kwarg.
'''
from matplotlib.patches import Ellipse
if 'rbody' not in self.attrs:
raise KeyError('rbody not found in self.attrs!')
dbody = 2.0 * self.attrs['rbody']
body = Ellipse((0,0),dbody,dbody,facecolor=facecolor, zorder=999,
**extra_kwargs)
if DoPlanet:
self.add_planet(ax, ang=ang)
if ax != None:
ax.add_artist(body)
def add_pcolor(self, dim1, dim2, value, zlim=None, target=None, loc=111,
title=None, xlabel=None, ylabel=None,
ylim=None, xlim=None, add_cbar=False, clabel=None,
add_body=True, dolog=False, *args, **kwargs):
'''
Create a pcolor plot of variable **value** against **dim1** on the
x-axis and **dim2** on the y-axis. Pcolor plots shade each
computational cell with the value at the cell center. Because no
interpolation or smoothing is used in the visualization, pcolor plots
are excellent for examining the raw output.
Simple example:
>>> self.add_pcolor('x', 'y', 'rho')
If kwarg **target** is None (default), a new figure is
generated from scratch. If target is a matplotlib Figure
object, a new axis is created to fill that figure at subplot
location **loc**. If **target** is a matplotlib Axes object,
the plot is placed into that axis.
Four values are returned: the matplotlib Figure and Axes objects,
the matplotlib contour object, and the matplotlib colorbar object
(defaults to *False* if not used.)
=========== ==========================================================
Kwarg Description
=========== ==========================================================
target Set plot destination. Defaults to new figure.
loc Set subplot location. Defaults to 111.
title Sets title of axes. Default is no title.
xlabel Sets x label of axes. Defaults to **dim1** and units.
ylabel Sets y label of axes. Defaults to **dim2** and units.
xlim Sets limits of x-axes. Defaults to whole domain.
ylim Sets limits of y-axes. Defaults to whole domain.
zlim Sets color bar range. Defaults to variable max/min.
add_cbar Adds colorbar to plot. Defaults to *False*.
clabel Sets colorbar label. Defaults to **var** and units.
add_body Places planetary body in plot. Defaults to **True**.
dolog Sets use of logarithmic scale. Defaults to **False**.
=========== ==========================================================
'''
import numbers
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,10), loc=loc)
# Get max/min if none given.
if zlim is None:
zlim=[0,0]
zlim[0]=self[value].min(); zlim[1]=self[value].max()
if dolog and zlim[0]<=0:
zlim[0] = np.min( [0.0001, zlim[1]/1000.0] )
# Logarithmic scale?
if dolog:
z=np.where(self[value]>zlim[0], self[value], 1.01*zlim[0])
norm=LogNorm()
else:
z=self[value]
norm=None
if self['grid'].attrs['gtype']=='Regular':
pass
else:
# Indices corresponding to QTree dimensions:
ix=self['grid'].attrs['dims'].index(dim1)
iy=self['grid'].attrs['dims'].index(dim2)
for k in self.qtree:
# Plot only leafs of the tree.
if not self.qtree[k].isLeaf: continue
leaf=self.qtree[k]
x=leaf.cells[ix]
y=leaf.cells[iy]
z_local=z[leaf.locs]
pcol=ax.pcolormesh(x,y,z_local,vmin=zlim[0],vmax=zlim[1],
norm=norm,**kwargs)
# Add cbar if necessary.
if add_cbar:
cbar=plt.colorbar(pcol, ax=ax, pad=0.01)
if clabel==None:
clabel="%s (%s)" % (value, self[value].attrs['units'])
cbar.set_label(clabel)
else:
cbar=None # Need to return something, even if none.
# Set title, labels, axis ranges (use defaults where applicable.)
if title: ax.set_title(title)
if ylabel==None: ylabel='%s ($R_{E}$)'%dim2.upper()
if xlabel==None: xlabel='%s ($R_{E}$)'%dim1.upper()
ax.set_ylabel(ylabel); ax.set_xlabel(xlabel)
try:
assert len(xlim)==2
assert isinstance(xlim[0], numbers.Number)
assert isinstance(xlim[1], numbers.Number)
except (TypeError, AssertionError):
if xlim is not None:
raise ValueError('add_pcolor: xlim must be list- or array-like and have 2 elements')
else:
ax.set_xlim(xlim)
try:
assert len(ylim)==2
assert isinstance(ylim[0], numbers.Number)
assert isinstance(ylim[1], numbers.Number)
except (TypeError, AssertionError):
if ylim is not None:
raise ValueError('add_pcolor: ylim must be list- or array-like and have 2 elements')
else:
ax.set_ylim(ylim)
# Add body/planet. Determine where the sun is first.
if dim1=='x':
ang=0.0
elif dim2=='x':
ang=90.0
else: ang=0.0
if add_body: self.add_body(ax, ang=ang)
return fig, ax, pcol, cbar
def add_contour(self, dim1, dim2, value, nlev=30, target=None, loc=111,
title=None, xlabel=None, ylabel=None,
ylim=None, xlim=None, add_cbar=False, clabel=None,
filled=True, add_body=True, dolog=False, zlim=None,
*args, **kwargs):
'''
Create a contour plot of variable **value** against **dim1** on the
x-axis and **dim2** on the y-axis.
Simple example:
>>> self.add_contour('x', 'y', 'rho')
If kwarg **target** is None (default), a new figure is
generated from scratch. If target is a matplotlib Figure
object, a new axis is created to fill that figure at subplot
location **loc**. If **target** is a matplotlib Axes object,
the plot is placed into that axis.
Four values are returned: the matplotlib Figure and Axes objects,
the matplotlib contour object, and the matplotlib colorbar object
(defaults to *False* if not used.)
=========== ==========================================================
Kwarg Description
=========== ==========================================================
target Set plot destination. Defaults to new figure.
loc Set subplot location. Defaults to 111.
nlev Number of contour levels. Defaults to 30.
title Sets title of axes. Default is no title.
xlabel Sets x label of axes. Defaults to **dim1** and units.
ylabel Sets y label of axes. Defaults to **dim2** and units.
xlim Sets limits of x-axes. Defaults to whole domain.
ylim Sets limits of y-axes. Defaults to whole domain.
zlim Sets contour range. Defaults to variable max/min.
add_cbar Adds colorbar to plot. Defaults to *False*.
clabel Sets colorbar label. Defaults to **var** and units.
add_body Places planetary body in plot. Defaults to **True**.
dolog Sets use of logarithmic scale. Defaults to **False**.
=========== ==========================================================
'''
import numbers
import matplotlib.pyplot as plt
from matplotlib.colors import (LogNorm, Normalize)
from matplotlib.ticker import (LogLocator, LogFormatter,
LogFormatterMathtext, MultipleLocator)
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,10), loc=loc)
# Get max/min if none given.
if zlim is None:
zlim=[0,0]
zlim[0]=self[value].min(); zlim[1]=self[value].max()
if dolog and zlim[0]<=0:
zlim[0] = np.min( [0.0001, zlim[1]/1000.0] )
# Set contour command based on grid type.
if self['grid'].attrs['gtype'] != 'Regular': # Non-uniform grids.
if filled:
contour=ax.tricontourf
else:
contour=ax.tricontour
else: # Uniform grids.
if filled:
contour=ax.contourf
else:
contour=ax.contour
# Create levels and set norm based on dolog.
if dolog:
levs = np.power(10, np.linspace(np.log10(zlim[0]),
np.log10(zlim[1]), nlev))
z=np.where(self[value]>zlim[0], self[value], 1.01*zlim[0])
norm=LogNorm()
ticks=LogLocator()
fmt=LogFormatterMathtext()
else:
levs = np.linspace(zlim[0], zlim[1], nlev)
z=self[value]
norm=None
ticks=None
fmt=None
# Plot contour.
cont=contour(self[dim1],self[dim2],np.array(z),
levs,*args, norm=norm, **kwargs)
# Add cbar if necessary.
if add_cbar:
cbar=plt.colorbar(cont, ax=ax, ticks=ticks, format=fmt, pad=0.01)
if clabel==None:
clabel="%s (%s)" % (value, self[value].attrs['units'])
cbar.set_label(clabel)
else:
cbar=None # Need to return something, even if none.
# Set title, labels, axis ranges (use defaults where applicable.)
if title: ax.set_title(title)
if ylabel==None: ylabel='%s ($R_{E}$)'%dim2.upper()
if xlabel==None: xlabel='%s ($R_{E}$)'%dim1.upper()
ax.set_ylabel(ylabel); ax.set_xlabel(xlabel)
try:
assert len(xlim)==2
assert isinstance(xlim[0], numbers.Number)
assert isinstance(xlim[1], numbers.Number)
except (TypeError, AssertionError):
if xlim is not None:
raise ValueError('add_contour: xlim must be list- or array-like and have 2 elements')
else:
ax.set_xlim(xlim)
try:
assert len(ylim)==2
assert isinstance(ylim[0], numbers.Number)
assert isinstance(ylim[1], numbers.Number)
except (TypeError, AssertionError):
if ylim is not None:
raise ValueError('add_contour: ylim must be list- or array-like and have 2 elements')
else:
ax.set_ylim(ylim)
# Add body/planet. Determine where the sun is first.
if dim1=='x':
ang=0.0
elif dim2=='x':
ang=90.0
else: ang=0.0
if add_body: self.add_body(ax, ang=ang)
return fig, ax, cont, cbar
class ShellSlice(IdlFile):
'''
Shell slices are special MHD outputs where the domain is interpolated
onto a spherical slice in 1, 2, or 3 dimensions. Some examples
include radial or azimuthal lines, spherical shells, or 3D wedges.
The *Shell* class reads and handles these output types.
'''
def __init__(self, filename, format='binary', *args, **kwargs):
from spacepy.pybats import parse_filename_time
IdlFile.__init__(self, filename, format=format, keep_case=False)
# Extract time from file name:
i_iter, runtime, time = parse_filename_time(self.attrs['file'])
if 'time' not in self.attrs: self.attrs['time'] = time
if 'iter' not in self.attrs: self.attrs['iter'] = i_iter
### Create some helper variables for plotting and calculations
d2r = np.pi/180. # Convert degrees to radians.
# Get grid spacing. If npoints ==1, set to 1 to avoid math errors.
self.drad = (self['r'][ -1] - self['r'][ 0])/max(self['grid'][0]-1,1)
self.dlon = (self['lon'][-1] - self['lon'][0])/max(self['grid'][1]-1,1)
self.dlat = (self['lat'][-1] - self['lat'][0])/max(self['grid'][2]-1,1)
self.dphi = d2r*self.dlon
self.dtheta = d2r*self.dlat
# Get spherical, uniform grid in units of r_body/radians:
self.lon, self.r, self.lat = np.meshgrid(
np.array(self['lon']), np.array(self['r']), np.array(self['lat']))
self.phi = d2r*self.lon
self.theta = d2r*(90-self.lat)
@calc_wrapper
def calc_urad(self):
'''
Calculate radial velocity.
'''
ur = self['ux']*np.sin(self.theta)*np.cos(self.phi) + \
self['uy']*np.sin(self.theta)*np.sin(self.phi) + \
self['uz']*np.cos(self.theta)
self['ur'] = dmarray(ur, {'units':self['ux'].attrs['units']})
@calc_wrapper
def calc_radflux(self, var, conv=1000. * (100.0)**3):
'''
For variable *var*, calculate the radial flux of *var* through each
grid point in the slice as self[var]*self['ur'].
Resulting value stored as "var_rflx".
'''
if var+'_rflx' in self: return
# Make sure we have radial velocity.
if 'ur' not in self: self.calc_urad()
# Calc flux:
self[var+'_rflx'] = self[var] * self['ur'] * conv
@calc_wrapper
def calc_radflu(self, var):
'''
For variable *var*, calculate the radial fluence, or the
spatially integrated radial flux through 2D surfaces of
constant radius.
Resulting variable stored as "var_rflu". Result will be an array
with one value for each radial distance within the object.
'''
# Need at least 2D in angle space:
if self.dphi==0 or self.dtheta==0:
raise ValueError('Fluence can only be calculated for 2D+ surfaces.')
# Trim flux off of val name:
if '_rflx' in var: var = var[:-5]
# Convenience:
flux = var + '_rflx'
flu = var + '_rflu'
# Make sure flux exists:
if flux not in self: self.calc_radflux(var)
if flu in self: return
# Create output container, one point per radial distance:
self[flu] = np.zeros( self['grid'][0] )
# Integrate over all radii.
# Units: convert R to km and cm-3 to km.
for i, R in enumerate(self['r']):
self[flu] = np.sum(
(R* 6371.0)**2 * self[flux][i,:,:] * \
np.sin(self.theta[i,:,:]) * \
self.dtheta * self.dphi * 1000.**2 )
def add_cont_shell(self, value, irad=0, target=None, loc=111,
zlim=None, dolabel=True, add_cbar=False,
dofill=True, nlev=51, colat_max=90,
rotate=np.pi/2, dolog=False, clabel=None,
latticks=15., yticksize=14, extend='both', **kwargs):
'''
For slices that cover a full hemisphere or more, create a polar
plot of variable *value*.
Extra keywords are sent to the matplotlib contour command.
'''
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib.colors import (LogNorm, Normalize)
from matplotlib.ticker import (LogLocator, LogFormatter,
LogFormatterMathtext, MultipleLocator)
fig, ax = set_target(target, figsize=(10,10), loc=loc, polar=True)
# Get max/min if none given.
if zlim is None:
zlim=[0,0]
zlim[0]=self[value][irad,:,:].min()
zlim[1]=self[value][irad,:,:].max()
# For log space, no negative zlimits.
if dolog and zlim[0]<=0:
zlim[0] = np.min( [0.0001, zlim[1]/1000.0] )
# Create levels and set norm based on dolog.
if dolog: # Log space!
levs = np.power(10, np.linspace(np.log10(zlim[0]),
np.log10(zlim[1]), nlev))
z=np.where(self[value]>zlim[0], self[value], 1.01*zlim[0])
norm=LogNorm()
ticks=LogLocator()
fmt=LogFormatterMathtext()
else:
levs = np.linspace(zlim[0], zlim[1], nlev)
z=self[value]
norm=None
ticks=None
fmt=None
# Select proper contour function based on fill/don't fill.
if dofill:
func = ax.contourf
else:
func = ax.contour
# Plot result. Rotate "rotate" radians to get sun in right spot.
# Plot against colatitude to arrange results correctly.
cnt = func(self.phi[irad,:,:]+rotate, 90-self.lat[irad,:,:],
np.array(self[value][irad,:,:]), levs, norm=norm,
extend=extend, **kwargs)
# Add cbar if necessary.
if add_cbar:
cbar=plt.colorbar(cnt, ax=ax, ticks=ticks, format=fmt, shrink=.85)
if clabel==None:
clabel="{} ({})".format(value, self[value].attrs['units'])
cbar.set_label(clabel)
else:
cbar=None # Need to return something, even if none.
# Adjust latitude
ax.set_ylim( [0, colat_max] )
# Adjust atitude and add better labels:
ax.yaxis.set_major_locator(MultipleLocator(latticks))
ax.set_ylim([0,colat_max])
ax.set_yticklabels('')
opts = {'size':yticksize, 'rotation':-45, 'ha':'center', 'va':'center'}
for theta in np.arange(90-latticks, 90-colat_max, -latticks):
txt = '{:02.0f}'.format(theta)+r'$^{\circ}$'
ax.text(pi/4., 90.-theta, txt, color='w', weight='extra bold',**opts)
ax.text(pi/4., 90.-theta, txt, color='k', weight='light', **opts)
# Use MLT-type labels.
lt_labels = ['Noon', '18', '00', '06']
xticks = [ 0, pi/2, pi, 3*pi/2]
xticks = np.array(xticks) + rotate
# Apply x-labels:
ax.set_xticks(xticks)
ax.set_xticklabels(lt_labels)
return fig, ax, cnt, cbar
[docs]class Mag(PbData):
'''
A container for data from a single BATS-R-US virtual magnetometer. These
work just like a typical :class:`spacepy.pybats.PbData` object. Beyond
raw magnetometer data, additional values are calculated and stored,
including total pertubations (the sum of all global and ionospheric
pertubations as measured by the magnetometer). Users will be interested
in methods :meth:`~spacepy.pybats.bats.Mag.add_comp_plot` and
:meth:`~spacepy.pybats.bats.Mag.calc_dbdt`.
Instantiation is best done through :class: `spacepy.pybats.MagFile`
objects, which load and parse organize many virtual magnetometers from a
single output file into a single object. However, they can be created
manually, though painfully. Users must instantiate by handing the
new object the number of lines that will be parsed (rather, the number
of data points that will be needed), a time vector, and (optionally)
the list of variables coming from the GM and IE module. While the
latter two are keyword arguments, at least one should be provided.
Next, the arrays whose keys were given by the *gmvars* and *ievars*
keyword arguments in the instantiation step can either be filled manually
or by using the :meth:`~spacepy.pybats.bats.Mag.parse_gmline` and
:meth:`~spacepy.pybats.bats.Mag.parse_ieline` methods to parse lines of
ascii data from a magnetometer output file. Finally, the
:meth:`~spacepy.pybats.bats.Mag.recalc` method should be called to
calculate total perturbation.
'''
def __init__(self, nlines, time, gmvars=(), ievars=(), *args, **kwargs):
from numpy import zeros
super(Mag, self).__init__(*args, **kwargs) # Init as PbData.
self['time']=time
self.attrs['nlines']=nlines
self['x']=np.zeros(nlines)
self['y']=np.zeros(nlines)
self['z']=np.zeros(nlines)
# Create IE and GM specific containers.
for key in gmvars:
self[key]=np.zeros(nlines)
for key in ievars:
self['ie_'+key]=np.zeros(nlines)
def parse_gmline(self, i, line, namevar):
'''
Parse a single line from a GM_mag*.out file and put into
the proper place in the magnetometer arrays. The line should
have the same number of variables as was initially given to the
:class:`~spacepy.pybats.bats.Mag` object. This method is best
used through the :class:`~spacepy.pybats.bats.MagFile` class interface.
Usage:
>>> self.parse_gmline(i, line, namevar)
where i is the entry number, line is the raw ascii line, and
namevar is the list of variable names.
'''
parts=line.split()
self['x'][i]=float(parts[9])
self['y'][i]=float(parts[10])
self['z'][i]=float(parts[11])
for j, key in enumerate(namevar):
self[key][i]=float(parts[j+12])
def parse_ieline(self, i, line, namevar):
'''
Parse a single line from a IE_mag*.out file and put into
the proper place in the magnetometer arrays. The line should
have the same number of variables as was initially given to the
:class:`~spacepy.pybats.bats.Mag` object. This method is best
used through the :class:`~spacepy.pybats.bats.MagFile` class interface.
Usage:
>>> self.parse_gmline(i, line, namevar)
where i is the entry number, line is the raw ascii line, and
namevar is the list of variable names.
'''
parts=line.split()
for j, key in enumerate(namevar):
self['ie_'+key][i]=float(parts[j+11])
def _recalc(self):
'''
Calculate total :math:`\Delta B` from GM and IE; store under object keys
*totaln*, *totale*, and *totald* (one for each component of the HEZ
coordinate system).
This function should only be called to correct legacy versions of
magnetometer files.
'''
from numpy import sqrt, zeros
# If values already exist, do not overwrite.
if 'dBn' in self: return
# New containers:
self['totaln']=np.zeros(self.attrs['nlines'])
self['totale']=np.zeros(self.attrs['nlines'])
self['totald']=np.zeros(self.attrs['nlines'])
for key in list(self.keys()):
if key[-2:]=='Bn':
self['totaln']=self['totaln']+self[key]
if key[-2:]=='Be':
self['totale']=self['totale']+self[key]
if key[-2:]=='Bd':
self['totald']=self['totald']+self[key]
# Old names -> new names:
varmap={'totaln':'dBn', 'totale':'dBe', 'totald':'dBd',
'gm_dBn':'dBnMhd', 'gm_dBe':'dBeMhd', 'gm_dBd':'dBdMhd',
'gm_facdBn':'dBnFac', 'gm_facdBe':'dBeFac','gm_facdBd':'dBdFac',
'ie_JhdBn':'dBnHal', 'ie_JhdBe':'dBeHal', 'ie_JhdBd':'dBdHal',
'ie_JpBn':'dBnPed', 'ie_JpBe':'dBePed', 'ie_JpBd':'dBdPed'}
# Replace variable names.
for key in list(self.keys()):
if key in varmap:
self[ varmap[key] ] = self.pop(key)
def calc_h(self):
'''
Calculate the total horizontal perturbation, 'H', using the pythagorean
sum of the two horizontal components (north-south and east-west
components):
$\Delta B_H = \sqrt{\Delta B_N^2 + \Delta B_E^2}$
'''
allvars = list(self.keys())
for v in allvars:
# Find all dB-north variables:
if v[:3] == 'dBn':
v_east = v.replace('dBn', 'dBe')
self[v.replace('dBn', 'dBh')] = dmarray(
np.sqrt(self[v]**2+self[v_east]**2), {'units':'nT'})
def calc_dbdt(self):
'''
Calculate the time derivative of all dB-like variables and save as
'dBdt[direction][component]. For example, the time derivative of
dBeMhd will be saved as dBdteMhd.
|dB/dt|_h is also calculated following the convention of
Pulkkinen et al, 2013:
$|dB/dt|_H = \sqrt{(\dB_N/dt)^2 + (dB_E/dt)^2}$
A 2nd-order accurate centeral difference method is used to
calculate the time derivative. For the first and last points,
2nd-order accurate forward and backward differences are taken,
respectively.
'''
# Do not calculate twice.
if 'dBdtn' in self: return
# Get dt values:
dt = np.array([x.total_seconds() for x in np.diff(self['time'])])
# Loop through variables:
oldvars = list(self.keys())
for k in oldvars:
if 'dB' not in k: continue
# Create new variable name and container:
new = k.replace('dB','dBdt')
self[new] = dmarray(np.zeros(self.attrs['nlines']),{'units':'nT/s'})
# Central diff:
self[new][1:-1] = (self[k][2:]-self[k][:-2])/(dt[1:]+dt[:-1])
# Forward diff:
self[new][0] = (-self[k][2] + 4*self[k][1] - 3*self[k][0]) \
/ (dt[1]+dt[0])
# Backward diff:
self[new][-1] = (3*self[k][-1] - 4*self[k][-2] + self[k][-3]) \
/ (dt[-1]+dt[-2])
self['dBdth'] = np.sqrt(self['dBdtn']**2+self['dBdte']**2)
def add_plot(self, value, style='-', target=None, loc=111, label=None,
**kwargs):
'''
Plot **value**, which should be a key corresponding to a data vector
stored in the :class:`~spacepy.pybats.bats.Mag` object,
against the object's *time*. The **target** kwarg specifies the
destination of the plot. If not set, **target** defaults to None and
a new figure and axis will be created. If target is a matplotlib
figure, a new axis is created at subplot location 111 (which can be
changed
using kwarg **loc**). If target is a matplotlib Axes object, the line
is added to the plot as if ``Axes.plot()`` was used. The line label,
which is used for setting labels on figure legends,
defaults to the value key but can be customized with the "label"
kwarg. The **style** keyword accepts basic Matplotlib line style
strings such as '--r' or '+g'. This string will be passed on to the
plot command to customize the line.
All extra kwargs are handed to ``Axes.plot``, allowing the user to set
any additional options (e.g., line color and style, etc.).
Three values are returned: the Figure object, Axis object, and
newly created line object. These can be used to further customize
the figure, axis, and line as necessary.
Example: Plot total :math:`\Delta B_N` onto an existing axis with line
color blue, line style dashed, and line label "Wow!":
>>> self.plot('dBn', target='ax', label='Wow!', lc='b', ls='--')
Example: Plot total :math:`\Delta B_N` on a new figure, save returned
values and overplot additional values on the returned axis. Default
labels and line styles are used in this example.
>>> fig, ax, line = self.plot('n')
>>> self.plot('dBe', target = ax)
'''
import matplotlib.pyplot as plt
if not label:
label=value
# Set figure and axes based on target:
fig, ax = set_target(target, figsize=(10,4), loc=loc)
line=ax.plot(self['time'], self[value], style, label=label, **kwargs)
applySmartTimeTicks(ax, self['time'], dolabel=True)
return fig, ax
def add_comp_plot(self, direc, target=None, add_legend=True,
loc=111, lw=2.0):
'''
Create a plot with, or add to an existing plot, an illustration of
how the
separate components sum together to make the total disturbance in a
given orthongal direction (arg **direc**). The three possible
components are 'n' (northwards, towards the magnetic pole), 'e'
(eastwards), or 'd' (downwards towards the center of the Earth.) The
components of the total disturbance in any on direction are
magnetospheric currents ('gm_dB'), gap-region field-aligned currents
('gm_facdB'), and ionospheric Hall and Pederson currents ('ie_Jp' and
'ie_Jh').
Example usage:
>>> self.add_comp_plot('n')
This will create a new plot with the total disturbance in the 'n'
direction along with line plots of each component that builds this
total. This method uses the familiar PyBats **target** kwarg system
to allow users to add these plots to existing figures or axes.
Parameters
==========
direc : string
Indicate the direction to plot: either 'n', 'e', 'd', or
'h' if calculated.
Other Parameters
================
target : Matplotlib Figure or Axes object
Set plot destination. Defaults to new figure.
loc : 3-digit integer
Set subplot location. Defaults to 111.
add_legend : bool
Add legend to plot. Defaults to True.
lw : float
Set the width of the lines. Defaults to 2.0 Total field
is always 1.5 times thicker.
'''
import matplotlib.pyplot as plt
# Set plot targets.
fig, ax = set_target(target, figsize=(10,4), loc=loc)
prefix = 'dB'+direc
# Use a dictionary to assign line styles, widths.
styles={prefix+'Mhd':'--', prefix+'Fac':'--',
prefix+'Hal':'-.', prefix+'Ped':'-.',
prefix:'-'}
widths={prefix+'Mhd':lw, prefix+'Fac':lw,
prefix+'Hal':lw, prefix+'Ped':lw,
prefix:1.5*lw}
colors={prefix+'Mhd':'#FF6600', prefix+'Fac':'r',
prefix+'Hal':'b', prefix+'Ped':'c',
prefix:'k'}
# Labels:
labels={prefix+'Mhd':r'$J_{Mag}$', prefix+'Fac':r'$J_{Gap}$',
prefix+'Hal':r'$J_{Hall}$', prefix+'Ped':r'$J_{Peder}$',
prefix:r'Total $\Delta B'+'_{}$'.format(direc)}
# Plot.
for k in sorted(self):
if ('dB'+direc not in k) or (k=='time') or (k[:4]=='dBdt'):continue
ax.plot(self['time'], self[k], label=labels[k],
lw=widths[k], c=colors[k])#,ls=styles[k]
# Ticks, zero-line, and legend:
applySmartTimeTicks(ax, self['time'], True, True)
ax.hlines(0.0, self['time'][0], self['time'][-1],
linestyles=':', lw=2.0, colors='k')
if add_legend: ax.legend(ncol=3, loc='best')
# Axis labels:
ax.set_ylabel(r'$\Delta B_{%s}$ ($nT$)'%(direc.upper()))
if target==None: fig.tight_layout()
return fig, ax
[docs]class MagFile(PbData):
'''
BATS-R-US magnetometer files are powerful tools for both research and
operations. :class:`~spacepy.pybats.bats.MagFile` objects open, parse,
and visualize such output.
The $\delta B$ calculated by the SWMF requires two components: GM (BATSRUS)
and IE (Ridley_serial). The data is spread across two files: GM_mag*.dat
and IE_mag*.dat. The former contains $\delta B$ caused by gap-region
(i.e., inside the inner boundary) FACs and the changing global field.
The latter contains the $\delta B$ caused by Pederson and Hall
currents in the ionosphere. :class:`~spacepy.pybats.bats.MagFile` objects
can open one or both of these files at a time; when both are opened, the
total $\delta B$ is calculated and made available to the user.
Usage:
>>> # Open up the GM magnetometer file only.
>>> obj = spacepy.pybats.bats.MagFile('GM_file.mag')
>>>
>>> # Open up both the GM and IE file [LEGACY SWMF ONLY]
>>> obj = spacepy.pybats.bats.MagFile('GM_file.mag', 'IE_file.mag')
>>>
>>> # Open up the GM magnetometer file; search for the IE file.
>>> obj = spacepy.pybats.bats.MagFile('GM_file.mag', find_ie=True)
Note that the **find_ie** kwarg uses a simple search assuming the data
remain in a typical SWMF-output organizational tree (i.e., if the results
of a simulation are in folder *results*, the GM magnetometer file can be
found in *results/GM/* or *results/GM/IO2/* while the IE file can be found
in *results/IE/* or *results/IE/ionosphere/*). It will also search the
present working directory. This method is not robust; the user must take
care to ensure that the two files correspond to each other.
'''
def __init__(self, filename, ie_name=None, find_ie=False, *args, **kwargs):
from glob import glob
super(MagFile, self).__init__(*args, **kwargs) # Init as PbData.
self.attrs['gmfile']=filename
# Try to find the IE file based on our current location.
if(find_ie and not ie_name):
basedir = filename[0:filename.rfind('/')+1]
if glob(basedir + '../IE/IE_mag_*.mag'):
self.attrs['iefile']=glob(basedir + '../IE/IE_mag_*.mag')[-1]
elif glob(basedir + '../../IE/ionosphere/IE_mag_*.mag'):
self.attrs['iefile'] = \
glob(basedir + '../../IE/ionosphere/IE_mag_*.mag')[-1]
elif glob(basedir + '/IE_mag_*.mag'):
self.attrs['iefile']= glob(basedir + '/IE_mag_*.mag')[-1]
else:
self.attrs['iefile']=None
else:
self.attrs['iefile']=ie_name
# Set legacy mode to handle old variable names:
self.legacy = find_ie or bool(ie_name)
self.readfiles()
def readfiles(self):
'''
Read and parse GM file and IE file (if name given.)
'''
import datetime as dt
from numpy import zeros
# Slurp lines.
infile = open(self.attrs['gmfile'], 'r')
lines=infile.readlines()
infile.close()
# Parse header.
# Get number of stations.
nmags=int((lines[0].split(':')[0]).split()[0])
# Get station names.
names=lines[0].split(':')[1]
namemag = names.split()
self.attrs['namemag']=namemag
# Check nmags vs number of mags in header.
if nmags != len(namemag):
raise BaseException(
'ERROR: GM file claims %i magnetomers, lists %i'
% (nmags, len(namemag)))
# Grab variable names. Use legacy mode if necessary:
prefix = 'gm_'*self.legacy
# skip time, iter, and loc; add prefix to var names:
gm_namevar = [prefix+x for x in lines[1].split()[12:]]
# Set number of mags and records.
self.attrs['nmag']=len(namemag)
nrecords = (len(lines)-2)//nmags
# If there is an IE file, Parse that header, too.
if self.attrs['iefile']:
infile=open(self.attrs['iefile'], 'r')
ielns =infile.readlines()
infile.close()
nmags=int((ielns[0].split(':')[0]).split()[0])
iestats=(ielns[0].split(':')[1]).split()
# Check nmags vs number of mags in header.
if nmags != len(iestats):
raise BaseException(
'ERROR: IE file claims %i magnetomers, lists %i'
% (nmags, len(namemag)))
if iestats != self.attrs['namemag']:
raise RuntimeError("Files do not have matching stations.")
ie_namevar=ielns[1].split()[11:]
self.attrs['ie_namevar']=ie_namevar
if (len(ielns)/self.attrs['nmag']) != (nrecords-1):
print('Number of lines do not match: GM=%d, IE=%d!' %
(nrecords-1, len(ielns)/self.attrs['nmag']))
nrecords=min(ielns, nrecords-1)
else:
ie_namevar=()
self.attrs['ie_namevar']=()
# Build containers.
self['time']=np.zeros(nrecords, dtype=object)
self['iter']=np.zeros(nrecords, dtype=float)
for name in namemag:
self[name]=Mag(nrecords, self['time'], gm_namevar, ie_namevar)
data_buffer=np.zeros((nrecords,nmags,(len(gm_namevar)+3)))
# Read file data.
for i in range(nrecords):
line = lines[i*nmags+2]
parts=line.split()
self['iter'][i]=parts[0]
self['time'][i]=dt.datetime(
int(parts[1]), #year
int(parts[2]), #month
int(parts[3]), #day
int(parts[4]), #hour
int(parts[5]), #minute
int(parts[6]), #second
int(parts[7])*1000 #microsec
)
for j in range(nmags):
line=lines[i*nmags+j+2]
if j>0:
parts=line.split()
values=[float(part) for part in parts[9:]]
data_buffer[i,j]=values
if self.attrs['iefile'] and i>0:
line=ielns[i*nmags+2]
self[namemag[0]].parse_ieline(i, line, ie_namevar)
for j in range(1, nmags):
self[namemag[j]].parse_ieline(i, ielns[i*nmags+j+2],
ie_namevar)
for j in range(nmags):
mag=self[namemag[j]]
mag['x']=data_buffer[:,j,0]
mag['y']=data_buffer[:,j,1]
mag['z']=data_buffer[:,j,2]
for k, key in enumerate(gm_namevar):
mag[key]=data_buffer[:,j,k+3]
# Sum up IE/GM components if necessary (legacy only):
if self.legacy: self._recalc()
# Get time res.
self.attrs['dt']=(self['time'][1]-self['time'][0]).seconds/60.0
def _recalc(self):
'''
Old magnetometer files had different variable names and did not
contain the total perturbation. This function updates variable names
and sums all contributions from all models/regions to get total
:math:`\Delta B`.
This function is only required for legacy results. New versions of
the SWMF include both GM, IE, and total perturbations in a single
file.
'''
for mag in self.attrs['namemag']:
self[mag]._recalc()
def calc_h(self):
'''
For each magnetometer object, calculate the horizontal component of
the perturbations using the pythagorean sum of the two horizontal
components (north-south and east-west components):
$\Delta B_H = \sqrt{\Delta B_N^2 + \Delta B_E^2}$
'''
for k in self:
if k=='time' or k=='iter': continue
self[k].calc_h()
def calc_dbdt(self):
'''
For each magnetometer object, calculate the horizontal component of
the perturbations.
|dB/dt|_h is also calculated following the convention of
Pulkkinen et al, 2013:
$|dB/dt|_H = \sqrt{(\dB_N/dt)^2 + (dB_E/dt)^2}$
'''
for k in self:
if k=='time' or k=='iter': continue
self[k].calc_dbdt()
class MagGridFile(IdlFile):
'''
Magnetometer grids are a recent addition to BATS-R-US: instead of
specifying a small set of individual stations, the user can specify a
grid of many stations spanning a latitude/longitude range. The files
are output in the usual :class:`spacepy.pybats.IdlFile` format. This
class handles the reading, manipulating, and visualization of these files.
'''
def __init__(self, *args, **kwargs):
import re
from spacepy.pybats import parse_filename_time
from spacepy.coordinates import Coords
# Initialize as an IdlFile.
super(MagGridFile, self).__init__(header=None, *args, **kwargs)
# Additional header parsing:
head = self.attrs['header']
match = re.search('\((\w+)\).*\[(\w+)\].*\[(\w+)\]', head)
coord, unit1, unit2 = match.groups()
self['grid'].attrs['coord']=coord
# Extract time from file name:
i_iter, runtime, time = parse_filename_time(self.attrs['file'])
if 'time' not in self.attrs: self.attrs['time'] = time
if 'iter' not in self.attrs: self.attrs['iter'] = i_iter
# Set units based on header parsing:
for v in self:
if v == 'grid':
continue
elif v in self['grid'].attrs['dims']:
self[v].attrs['units']=unit1
else:
self[v].attrs['units']=unit2
# Get cooridnates in both SM and GEO.
#if coord == 'GEO':
# self['Lat_geo']=self['Lat']
# self['Lon_geo']=self['Lon']
@calc_wrapper
def calc_h(self):
'''
Calculate the total horizontal perturbation, 'h', using the pythagorean
sum of the two horizontal components (east and west components).
'''
allvars = list(self.keys())
for v in allvars:
# Find all dB-north variables:
if v[:3] == 'dBn':
v_east = v.replace('dBn', 'dBe')
self[v.replace('dBn', 'dBh')] = dmarray(
np.sqrt(self[v]**2+self[v_east]**2), {'units':'nT'})
def add_contour(self, value, nlev=30, target=None, loc=111,
title=None, xlabel=None, ylabel=None,
ylim=None, xlim=None, add_cbar=False, clabel=None,
filled=True, dolog=False, zlim=None,
coords=None, add_conts=False,
*args, **kwargs):
'''
Add a lat-lon contour plot of variable *value*.
'''
from spacepy.pybats import mhdname_to_tex
import matplotlib.pyplot as plt
from matplotlib.colors import (LogNorm, Normalize)
from matplotlib.ticker import (LogLocator, LogFormatter, FuncFormatter,
LogFormatterMathtext, MultipleLocator)
# Set ax and fig based on given target.
fig, ax = set_target(target, figsize=(10,7), loc=loc)
# Get max/min if none given.
if zlim is None:
zlim = [self[value].min(), self[value].max()]
if zlim[1]-zlim[0]>np.max(zlim):
maxval = max(abs(zlim[0]), abs(zlim[1]))
zlim=[-maxval, maxval]
# No zero-level for log scales:
if dolog and zlim[0]<=0:
zlim[0] = np.min( [0.0001, zlim[1]/1000.0] )
# Better default color maps:
if 'cmap' not in kwargs:
# If zlim spans positive and negative:
if zlim[1]-zlim[0]>np.max(zlim):
kwargs['cmap']='bwr'
else:
kwargs['cmap']='Reds'
# Set contour command based on filled/unfilled contours:
if filled:
contour=ax.contourf
else:
contour=ax.contour
# Create levels and set norm based on dolog.
if dolog:
levs = np.power(10, np.linspace(np.log10(zlim[0]),
np.log10(zlim[1]), nlev))
z=np.where(self[value]>zlim[0], self[value], 1.01*zlim[0])
norm=LogNorm()
ticks=LogLocator()
fmt=LogFormatterMathtext()
else:
levs = np.linspace(zlim[0], zlim[1], nlev)
z=self[value]
norm=None
ticks=MultipleLocator((zlim[1]-zlim[0])/10) ### fix this
fmt=None
# Add Contour to plot:
cont = contour(self['Lon'], self['Lat'], np.array(np.transpose(z)),
levs, *args, norm=norm, **kwargs)
# Add cbar if necessary.
if add_cbar:
cbar=plt.colorbar(cont, ax=ax, ticks=ticks, format=fmt, pad=0.01)
if clabel==None:
varname = mhdname_to_tex(value)
units = mhdname_to_tex(self[value].attrs['units'])
clabel="%s (%s)" % (varname, units)
cbar.set_label(clabel)
else:
cbar=None # Need to return something, even if none.
# Set title, labels, axis ranges (use defaults where applicable.)
if title: ax.set_title(title)
coord_sys = self['grid'].attrs['coord']
if ylabel==None: ylabel='Latitude ({})'.format(coord_sys)
if xlabel==None: xlabel='Longitude ({})'.format(coord_sys)
ax.set_ylabel(ylabel); ax.set_xlabel(xlabel)
if type(xlim) != None: ax.set_xlim(xlim)
if type(ylim) != None: ax.set_ylim(ylim)
# If a brand-new axes was created, use custom ticks:
if not issubclass(type(target), plt.Axes):
fmttr = FuncFormatter(lambda x, pos:
'{:4.1f}'.format(x)+r'$^{\circ}$')
ax.xaxis.set_major_formatter(fmttr)
ax.yaxis.set_major_formatter(fmttr)
# If a brand-new figure was created, use tight-layout.
if target==None: fig.tight_layout()
return fig, ax, cont, cbar
def interp(self, var, lons, lats):
'''
Interpolate variable *var* to point(s) *lats*, *lons*. *lats* and
*lons* can either be scalar or a numpy-like array. Values interpolated
to these positions is returned to the caller either as a scalar or
array.
Simple bilinear interpolation is used.
'''
from spacepy.pybats.batsmath import interp_2d_reg as intp
return intp(lats, lons, self['Lat'], self['Lon'], self[var])
[docs]class GeoIndexFile(LogFile):
'''
Geomagnetic Index files are a specialized BATS-R-US output that contain
geomagnetic indices calculated from simulated ground-based magnetometers.
Currently, the only index instituted is Kp through the faKe_p setup.
Future work will expand the system to include Dst, AE, etc.
GeoIndFiles are a specialized subclass of pybats.LogFile. It includes
additional methods to quickly visualize the output, perform data-model
comparisons, and more.
'''
def __repr__(self):
return 'GeoIndexFile object at %s' % (self.attrs['file'])
def __init__(self, filename, keep_case=True, *args, **kwargs):
'''
Load ascii file located at self.attrs['file'].
'''
# Call super __init__:
super(GeoIndexFile, self).__init__(filename, *args, **kwargs)
# Re-parse the file header; look for key info.
head = (self.attrs['descrip']).replace('=',' ')
parts = head.split()
if 'DtOutput=' in head:
self.attrs['dt'] = float(parts[parts.index('DtOutput=')+1])
if 'SizeKpWindow' in head:
self.attrs['window'] = float(parts[parts.index(
'SizeKpWindow(Mins)')+1])
if 'Lat' in head:
self.attrs['lat'] = float(parts[parts.index('Lat')+1])
if 'K9' in head:
self.attrs['k9'] = float(parts[parts.index('K9') +1])
def fetch_obs_kp(self):
'''
Fetch the observed Kp index for the time period covered in the
logfile. Return *True* on success.
Observed Kp is automatically fetched from the Kyoto World Data Center
via the :mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoKp` object, which holds the observed
kp, is stored as *self.obs_kp* for future use.
'''
import spacepy.pybats.kyoto as kt
# Return if already obtained:
if hasattr(self, 'obs_kp'): return True
# Start and end time to collect observations:
stime = self['time'][0]; etime = self['time'][-1]
# Attempt to fetch from Kyoto website:
try:
self.obs_kp = kt.fetch('kp', stime, etime)
# Warn on failure:
except BaseException as args:
raise Warning('Failed to fetch Kyoto Kp: ', args)
return False
return True
def fetch_obs_ae(self):
'''
Fetch the observed AE index for the time period covered in the
logfile. Return *True* on success.
Observed AE is automatically fetched from the Kyoto World Data Center
via the :mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoAe` object, which holds the observed
AE, is stored as *self.obs_ae* for future use.
'''
import spacepy.pybats.kyoto as kt
# Return if already obtained:
if hasattr(self, 'obs_ae'): return True
# Start and end time to collect observations:
stime = self['time'][0]; etime = self['time'][-1]
# Attempt to fetch from Kyoto website:
try:
self.obs_ae = kt.fetch('ae', stime, etime)
# Warn on failure:
except BaseException as args:
raise Warning('Failed to fetch Kyoto AE: ', args)
return False
return True
def add_kp_quicklook(self, target=None, loc=111, label=None,
plot_obs=False, add_legend=True,
obs_kwargs={'c':'k', 'ls':'--', 'lw':2}, **kwargs):
'''
Similar to "dst_quicklook"-type functions, this method fetches observed
Kp from the web and plots it alongside the Kp read from the GeoInd file.
Usage:
>>> obj.kp_quicklook(target=SomeMplTarget)
The target kwarg works like in other PyBats plot functions: it can be
a figure, an axes, or None, and it determines where the plot is placed.
Other kwargs customize the line. Label defaults to fa$K$e$_{P}$, extra
kwargs are passed to pyplot.plot.
Observed Kp can be added via the *plot_obs* kwarg. Kp is automatically
fetched from the Kyoto World Data Center via the
:mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoKp` object, which holds the observed
Kp, is stored as *self.obs_kp* for future use.
The observed line can be customized via the *obs_kwargs* kwarg, which
is a dictionary of plotting keyword arguments.
'''
import matplotlib.pyplot as plt
# Set up plot target.
fig, ax = set_target(target, figsize=(10,4), loc=loc)
# Create label:
if not(label):
label = 'fa$K$e$_{P}$'
if ('lat' in self.attrs) and ('k9' in self.attrs):
label += ' (Lat=%04.1f$^{\circ}$, K9=%03i)' % \
(self.attrs['lat'], self.attrs['k9'])
# Sometimes, the "Kp" varname is caps, sometimes not.
kp = 'Kp'
if kp not in self.keys(): kp='kp'
ax.plot(self['time'], self['Kp'], label=label,**kwargs)
ax.set_ylabel('$K_{P}$')
ax.set_xlabel('Time from '+ self['time'][0].isoformat()+' UTC')
applySmartTimeTicks(ax, self['time'])
if target==None: fig.tight_layout()
if plot_obs:
# Attempt to fetch Kp:
if self.fetch_obs_kp():
# If successful, add to plot:
self.obs_kp.add_histplot(target=ax, **obs_kwargs)
applySmartTimeTicks(ax, self['time'])
if add_legend: ax.legend(loc='best')
return fig, ax
def add_ae_quicklook(self, target=None, loc=111, label=None,
plot_obs=False, val='AE', add_legend=True,
obs_kwargs={'c':'k', 'ls':'--', 'lw':1.5}, **kwargs):
'''
Similar to "dst_quicklook"-type functions, this method fetches observed
AE indices from the web and plots it alongside the corresponding
AE read from the GeoInd file. Because there are four AE-like indices
(AL, AU, AE, and AO), the kwarg *val* specifies which to plot
(default is AE).
Usage:
>>> obj.ae_quicklook(target=SomeMplTarget)
The target kwarg works like in other PyBats plot functions: it can be
a figure, an axes, or None, and it determines where the plot is placed.
Other kwargs customize the line. Extra kwargs are passed to pyplot.plot
Observed AE can be added via the *plot_obs* kwarg. AE is automatically
fetched from the Kyoto World Data Center via the
:mod:`spacepy.pybats.kyoto` module. The associated
:class:`spacepy.pybats.kyoto.KyotoAe` object, which holds the observed
AE, is stored as *self.obs_kp* for future use.
The observed line can be customized via the *obs_kwargs* kwarg, which
is a dictionary of plotting keyword arguments.
'''
import matplotlib.pyplot as plt
# Set up plot target.
fig, ax = set_target(target, figsize=(10,4), loc=loc)
if not(label):
label = 'Virtual {}'.format(val)
ax.plot(self['time'], self[val], label=label,**kwargs)
ax.set_ylabel('{} ($nT$)'.format(val))
ax.set_xlabel('Time from '+ self['time'][0].isoformat()+' UTC')
applySmartTimeTicks(ax, self['time'])
if target==None: fig.tight_layout()
if plot_obs:
if self.fetch_obs_ae():
ax.plot(self.obs_ae['time'], self.obs_ae[val.lower()],
label='Obs.', **obs_kwargs)
applySmartTimeTicks(ax, self['time'])
if add_legend: ax.legend(loc='best')
return fig, ax
[docs]class VirtSat(LogFile):
'''
A :class:`spacepy.pybats.LogFile` object tailored to virtual satellite
output; includes special satellite-specific plotting methods.
'''
def __init__(self, *args, **kwargs):
from re import findall
from scipy.interpolate import interp1d
from matplotlib.dates import date2num
super(VirtSat, self).__init__(keep_case=False, *args, **kwargs)
# Attempt to extract the satellite's name and save it.
try:
s = self.attrs['file']
a = findall('sat_([\-\w]+)_\d+_n\d+\.sat|sat_(\w+)_n\d+\.sat', s)[0]
name = list(filter(None, a))[0]
except IndexError:
name = None
self.attrs['name']=name
# Create interpolation functions for position.
self._interp={}
tnum=date2num(self['time'])
for x in 'xyz':
self._interp[x] = interp1d(tnum, self[x])
def __repr__(self):
return 'Satellite object {} at {}'.format(self.attrs['name'],
self.attrs['file'])
def calc_ndens(self):
'''
Calculate number densities for each fluid. Species mass is ascertained
via recognition of fluid name (e.g. OpRho is clearly oxygen). A full
list of recognized fluids/species can be found by exploring the
dictionary *mass* found in :mod:`~spacepy.pybats.bats`. Composition is
also calculated as percent of total number density.
New values are saved using the keys *speciesN* (e.g. *opN*) and
*speciesFrac* (e.g. *opFrac*).
'''
_calc_ndens(self)
def calc_temp(self, units='eV'):
'''
Calculate plasma temperature for each fluid. Number density is
calculated using *calc_ndens* if it hasn't been done so already.
Temperature is obtained via density and pressure through the simple
relationship P=nkT.
Use the units kwarg to set output units. Current choices are
KeV, eV, and K. Default is eV.
'''
from spacepy.datamodel import dmarray
units = units.lower()
# Create dictionary of unit conversions.
conv = {'ev' : 6241.50935, # nPa/cm^3 --> eV.
'kev': 6.24150935, # nPa/cm^3 --> KeV.
'k' : 72429626.47} # nPa/cm^3 --> K.
# Calculate number density if not done already.
if not 'N' in self:
self.calc_ndens()
# Find all number density variables.
for key in list(self.keys()):
# Next variable if not number density:
if key[-1] != 'N':
continue
# Next variable if no matching pressure:
if not key[:-1]+'p' in self:
continue
self[key[:-1]+'t'] = dmarray(
conv[units] * self[key[:-1]+'p']/self[key],
attrs = {'units':units})
def calc_bmag(self):
'''
Calculates total magnetic field magnitude via:
$|B|=\sqrt{B_X^2+B_Y^2+B_Z^2}$. Results stored as variable name *b*.
'''
if 'b' not in self:
self['b']=np.sqrt(self['bx']**2+
self['by']**2+
self['bz']**2)
return True
def calc_magincl(self, units='deg'):
'''
Magnetic inclination angle (a.k.a. inclination angle) is defined:
$\Theta = sin^{-1}(B_Z/B)$. It is a crucial value when examining
magnetic dynamics about geosychronous orbit, the tail, and other
locations.
This function calculates the magnetic inclination for the Virtual
Satellite object and saves it as *b_incl*. Units default to degrees;
the keyword **units** can be changed to 'rad' to change this.
'''
if 'b' not in self: self.calc_bmag()
incl = np.arcsin(self['bz']/self['b'])
if units=='deg':
self['b_incl'] = dmarray(incl*180./np.pi, {'units':'degrees'})
elif units=='rad':
self['b_incl'] = dmarray(incl, {'units':'radians'})
else:
raise ValueError('Unrecognized units. Use "deg" or "rad"')
return True
def get_position(self, time):
'''
For an arbitrary time, *time*, return a tuple of coordinates for
the satellite's position at that time in GSM coordinates.
*time* should be type datetime.datetime, matplotlib.dates.date2num
output (i.e., number of days (fraction part represents hours,
minutes, seconds) since 0001-01-01 00:00:00 UTC, *plus* *one*), or
a sequence of either.
The satellite's position is interpolated (linearly) to *time*.
'''
from matplotlib.dates import date2num
# Test if sequence.
if isinstance(time, (list, np.ndarray)):
testval=time[0]
else:
testval=time
# Test if datetime or not.
if type(testval) == type(self['time'][0]):
time=date2num(time)
# Interpolate, using "try" as to not pass time limits and extrapolate.
try:
loc = (self._interp['x'](time),
self._interp['y'](time),
self._interp['z'](time))
except ValueError:
loc = (None, None, None)
return loc
def add_sat_loc(self, time, target, plane='XY', dobox=False,
dolabel=False, size=12, c='k', **kwargs):
'''
For a given axes, *target*, add the location of the satellite at
time *time* as a circle. If kwarg *dolabel* is True, the satellite's
name will be used to label the dot. Optional kwargs are any accepted by
matplotlib.axes.Axes.plot. The kwarg *plane* specifies the plane
of the plot, e.g., 'XY', 'YZ', etc.
'''
# Get Axes' limits:
xlim = target.get_xlim()
ylim = target.get_ylim()
plane=plane.lower()
loc = self.get_position(time)
if None in loc: return
x=loc['xyz'.index(plane[0])]
y=loc['xyz'.index(plane[1])]
# Do not label satellite if outside axes bounds.
if (x<min(xlim))or(x>max(xlim))or \
(y<min(ylim))or(y>max(ylim)): dolabel=False
target.plot(x, y, 'o', **kwargs)
if dolabel:
xoff = 0.03*(xlim[1]-xlim[0])
if dobox:
target.text(x+xoff,y,self.attrs['name'],
bbox={'fc':'w','ec':'k'},size=size, va='center')
else:
target.text(x+xoff,y,self.attrs['name'],
size=size, va='center', color=c)
# Restore Axes' limits.
target.set_ylim(ylim)
target.set_xlim(xlim)
def add_orbit_plot(self, plane='XY', target=None, loc=111, rbody=1.0,
title=None, trange=None, add_grid=True, style='g.',
adjust_axes=True, add_arrow=True,
arrow_kwargs={'color':'g', 'width':.05, 'ec':'k'},
**kwargs):
'''
Create a 2D orbit plot in the given plane (e.g. 'XY' or 'ZY').
Extra kwargs are handed to the plot function.
Parameters
==========
None
Other Parameters
================
plane : string
Set the plane in which to plot (XY, XZ, etc.) Defaults to XY.
target : Matplotlib Figure or Axes object
Set plot destination. Defaults to new figure.
loc : 3-digit integer
Set subplot location. Defaults to 111.
adjust_axes : bool
If True, axes will be customized to best display orbit (equal
aspect ratio, grid on, planet drawn, etc.). Defaults to True.
style: string
A matplotlib line style specifier. Defaults to 'g.'.
title : string
Set title of axes.
rbody : real
Set radius of model inner boundary. Defaults to 1.0
trange : list of datetimes
Set the time range to plot. Defaults to None, or entire dataset.
add_grid : boolean
Turn on or off grid style of axes. Default is True.
add_arrow : boolean
Add arrow at end of orbit path to indicate direction.
Default is True.
arrow_kwargs : dict
Dictionary of arrow kwargs. Defaults to match line style.
'''
from spacepy.pybats import add_body
from spacepy.pybats.ram import grid_zeros
import matplotlib.pyplot as plt
fig, ax = set_target(target, figsize=(5,5), loc=loc)
plane=plane.upper()
# Set time range of plot.
if not trange: trange = [self['time'].min(), self['time'].max()]
tloc = (self['time']>=trange[0])&(self['time']<=trange[-1])
# Extract orbit X, Y, or Z.
plane=plane.lower()
if plane[0] in ['x','y','z']:
x = self[plane[0]][tloc]
else:
raise ValueError('Bad dimension specifier: ' + plane[0])
if (plane[1] in ['x','y','z']) and (plane[0]!=plane[1]):
y = self[plane[1]][tloc]
else:
raise ValueError('Bad dimension specifier: ' + plane[1])
# Actually plot:
ax.plot(x, y, style, **kwargs)
# Add arrow to plot:
if add_arrow:
x_arr, y_arr = x[-1], y[-1]
dx, dy = x[-1]-x[-2], y[-1] - y[-2]
ax.arrow(x_arr, y_arr, dx, dy, **arrow_kwargs)
# Finish customizing axis.
if adjust_axes:
ax.axis('equal')
ax.set_xlabel('GSM %s'%(plane[0].upper()))
ax.set_ylabel('GSM %s'%(plane[1].upper()))
if title:
ax.set_title(title)
if add_grid:
ax.grid()
grid_zeros(ax)
add_body(ax, rad=rbody, add_night=('X' in plane.upper()) )
return fig, ax
