Reproducible LMA research with the IPython notebook and brawl4d

This notebook demonstrates how to display data included in the lmatools repository.

If you haven't yet, process and view sample data included with lmatools

python ~/code/lmatools/testing/ /path/to/output/files/

Then, edit the data_path in the second cell to include /path/to/output/files/ as defined above. Run all the cells prior to the "Charge Analysis" and try interacting with the plot. You should see some data.

In [5]:
%matplotlib qt4
# import matplotlib
# matplotlib.use('nbagg')
#import matplotlib.pyplot as plt
from brawl4d.brawl4d import B4D_startup, redraw
import os

In [8]:
data_path = '/data/tmp/flash_sort_test/'
lma_file = os.path.join(data_path, 'h5_files/2014/May/26/LYLOUT_140526_094000_0600.dat.flash.h5')

In [9]:

In the cell below, note that the basedate has been set to match the dataset we specified above.

If you are not using data from the WTLMA, then you'll also need to pass ctr_lon=value and ctr_lat=value to B4D_startup.

In [10]:
from datetime import datetime
panels = B4D_startup(basedate=datetime(2014,5,26), ctr_lat=33.5, ctr_lon=-101.5)

In [11]:
import matplotlib.pyplot as plt;

Below, set a valid path to lma_file. IPython will tab-complete paths, like the shell.

In [12]:
from brawl4d.LMA.controller import LMAController
lma_ctrl = LMAController()
d, post_filter_brancher, scatter_ctrl, charge_lasso = lma_ctrl.load_hdf5_to_panels(panels, lma_file)

found flash data

In [16]:
panels.panels['tz'].axis((9*3600 + 40*60, 9*3600 + 42*60, 1, 15))
panels.panels['xy'].axis((-130, 20, -10, 140))

(-130, 20, -10, 140)

In [14]:
from brawl4d.LMA.widgets import LMAwidgetController
from IPython.display import display
from brawl4d.LMA.controller import LMAController

lma_tools = LMAwidgetController(panels, lma_ctrl, scatter_ctrl, charge_lasso, d)

The cell above contains and combines all functions of Brawl4d into a single .py file. This allows all tools to be activiated simultaneously, while also enabling a centralized container for all active widgets for convienence.

LMA Tools Contained:

Number of Stations: Specify accordingly with data file (min=1; max=11)

Max Chi2: Values for chi2 obtained from the data file (min=0.0; max=1.0)

Charge Selection: Selection for Negative (-1), Neutral (0), and Positive (1) charge for charge selection and analyzation in the browser. The draw button activates the lasso tool enabling charge selection; re-clicking the draw button is necessary upon each selection made.

Color By: Allows the display of LMA data by chi2, time, or charge; selecting one will redraw the plot in the browser.

Animation Time: Allows for animation of the LMA data in the browser for charge polarity determination made by Charge Selection. The slider allows the user to select a desired time for total animation duration (min=1s, max=30s). Clicking Animate will then animate the data after a desired time has been selected.

Zoom in on a few cells of interest. The smaller, western and northern cells here are anomalously electrified, while the larger cluster is normally electrified.

Charge analysis

It's possible to use a lasso to classify charge regions inferred from LMA data. Set the polarity and run the code below to start the lasso. On the plot, left click to draw the lasso, and then right click to close the lasso and assign the charge.

If you're using an HDF5-format LMA data file, the analyzed charge is automatically written to the HDF5 file. The results of the operation can be queried by looking for the points that have had their charge set to the value defined above.

In [ ]:
chg =['charge']
wh = np.where(chg > 0)

Color by...

In addition to coloring the scatter plots by time, it's possible to use other values in the LMA data array.

In [ ]:
# A reference to the current data in the view is cached by the charge lasso.
current_data = charge_lasso.cache_segment.cache[-1]
# Manually set the color limits on the flash_id variable
scatter_ctrl.default_color_bounds.flash_id =(current_data['flash_id'].min(), current_data['flash_id'].max())
# Color by flash ID.
scatter_ctrl.color_field = 'flash_id'


Flash statistics

If the LMA controller found flash data, then it's possible to get a live update of flashes in the current view. current_events_flashes is an analysis pipeline branchpoint, which will send events and flashes to another analysis pipeline segment that can be specified with current_events_flashes.targets.add(target). Behind the scenes, it's hooked up to an segment that receives the events and flashes, and prints the average flash area of all flashes that have more than a threshold number of points.

Change the view a few times and you'll see updated flash stats below.

In [ ]:
current_events_flashes = lma_ctrl.flash_stats_for_dataset(d, scatter_ctrl.branchpoint)

In [ ]:

Flash volume for current sources

In [ ]:
from scipy.spatial import Delaunay
from scipy.misc import factorial

from stormdrain.pipeline import coroutine
class LMAEventStats(object):
    def __init__(self, GeoSys):
        """ GeoSys is an instance of
        self.GeoSys = GeoSys
    def ECEF_coords(self, lon, lat, alt):
        x,y,z = self.GeoSys.toECEF(lon, lat, alt)
        return x,y,z
    def _hull_volume(self):
        tri = Delaunay(self.xyzt[:,0:3])
        vertices = tri.points[tri.vertices]
        # This is the volume formula in 
        # Except the formula needs to be divided by ndim! to get the volume, cf., 
        # Credit Pauli Virtanen, Oct 14, 2012, scipy-user list
        q = vertices[:,:-1,:] - vertices[:,-1,None,:]
        simplex_volumes = (1.0 / factorial(q.shape[-1])) * np.fromiter(
                (np.linalg.det(q[k,:,:]) for k in range(tri.nsimplex)) , dtype=float)
        self.tri = tri
        # The simplex volumes have negative values since they are oriented 
        # (think surface normal direction for a triangle
    def events_flashes_receiver(self):
        while True:
            evs, fls = (yield)
            x,y,z = self.ECEF_coords(evs['lon'], evs['lat'], evs['alt'])
            t = evs['time']
            self.xyzt = np.vstack((x,y,z,t)).T
            print "Volume of hull of points in current view is {0:5.1f}".format(
                        self.volume / 1.0e9) # (1000 m)^3

In [ ]:
stats = LMAEventStats(panels.cs.geoProj)
stat_maker = stats.events_flashes_receiver()

In [ ]:

In [ ]:
print current_events_flashes.targets

Leader speed analysis

Run the cell below after adjusting the plot to plot the distance with time from the first point.

  • Negative leaders: $10^5$ m/s (Behnke et al., 2005, GRL)
  • Positive leaders: $10^4$ m/s (Williams 2008 ICAE, after Les Renardières Group 1977, 81)

In [ ]:
import scipy.spatial.distance as distance
all_dist_pairs = distance.pdist(stats.xyzt[:,0:2])

shift_t  = stats.xyzt[:,3]-stats.xyzt[0,3]

fig = plt.figure()
ax.scatter(shift_t, sqd[0,:], cmap='viridis')
t0, t1 = shift_t.min(), shift_t.max()
d0 = 0.0
d_c = 3.0e8*(t1-t0)
d_8 = 1.0e8*(t1-t0)
d_7 = 1.0e7*(t1-t0)
d_6 = 1.0e6*(t1-t0)
d_5 = 1.0e5*(t1-t0)
d_4 = 1.0e4*(t1-t0)
# ax.plot((t0, t1), (d0, d_c), label='c')
# ax.plot((t0, t1), (d0, d_8), label='1e8')
# ax.plot((t0, t1), (d0, d_7), label='1e7')
# ax.plot((t0, t1), (d0, d_6), label='1e6')
ax.plot((t0, t1), (d0, d_5), label='1e5')
ax.plot((t0, t1), (d0, d_4), label='1e4')

3D visualization of current sources with Mayavi

This will receive the events and flashes from the current 2D view, and update an interactive 3D view.

It should be possible to use the scatter_ctrl.branchpoint and the mappablerangeupdater setup by scatter_ctrl instead of the raw events flashes. This would permit synchronization with the same vmin, vmax, which is already being figured out.

In [ ]:
from brawl4d.brawl4d import redraw
import mayavi.mlab as mvlab
from stormdrain.pipeline import coroutine
class MayaviOutlet(object):
    def __init__(self, panels, ev_fl_broadcaster):
        self.ev_fl_broadcaster = ev_fl_broadcaster
        self.p3d = mvlab.points3d([0], [0], [0], [0], scale_factor=5e-5)
        self.scene = self.p3d.scene
        self.scene.background = (0,0,0)
        # Force a reflow of data
        # move camera to see everything after data are plotted
    def rx(self):
        while True:
            ev, fl = (yield)
#             self.ev = ev
#             self.fl = fl
            evx, evy, evz, evt = ev['x'], ev['y'], ev['z'], ev['time']
            self.p3d.mlab_source.reset(x=evx, y=evy, z=evz, scalars=evt)
current_events_flashes = lma_ctrl.flash_stats_for_dataset(d, scatter_ctrl.branchpoint)
mvo = MayaviOutlet(panels, current_events_flashes)

In [ ]:

In [ ]:

In [ ]: