Reproducible LMA research with the IPython notebook and brawl4d

This notebook demonstrates how to download and display data from the NSF-sponsored Deep Convective Clouds and Chemistry campaign.

To download data, go to the DC3 data archive and choose one of the LMA datasets. This example will assume use of the LMA VHF source and flash data on June 4, 2012 from 2050-2100 UTC from West Texas, as retrieved using the DC3 data download form.

Note: the hour of data below is 960 MB, so caveat downloader. Try just the 2050--2100 interval.

Running the code below is some basic boilerplate.


In [1]:
%pylab
from brawl4d.brawl4d import B4D_startup

# this function forces a manual redraw / re-flow of the data to the plot.
def redraw():
    from stormdrain.pubsub import get_exchange
    get_exchange('SD_bounds_updated').send(panels.bounds)


Welcome to pylab, a matplotlib-based Python environment [backend: Qt4Agg].
For more information, type 'help(pylab)'.
/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/matplotlib/__init__.py:908: UserWarning:  This call to matplotlib.use() has no effect
because the the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  if warn: warnings.warn(_use_error_msg)

In the cell below, note that the basedate has been set to match the dataset we downloaded 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 [2]:
from datetime import datetime
panels = B4D_startup(basedate=datetime(2012,6,4), ctr_lat=33.5, ctr_lon=-101.5)

Below, set a valid path to lma_file. The IPython notebook will try to tab-complete paths.


In [3]:
from brawl4d.LMA.controller import LMAController
lma_file = '/data/20120604/LYLOUT_120604_205000_0600.dat.flash.h5'
lma_ctrl = LMAController()
d, post_filter_brancher, scatter_ctrl, charge_lasso = lma_ctrl.load_hdf5_to_panels(panels, lma_file)


found flash data

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.


In [4]:
panels.panels['tz'].axis((20*3600 + 50*60, 20*3600 + 52*60, 1, 15))
panels.panels['xy'].axis((-40, 0, 50, 90))


Out[4]:
(-40, 0, 50, 90)

After adjusting the lower bound on station count below, it's necessary to trigger a reflow of the data down the pipeline. You can do this by dragging the plot to adjust the view bounds. Or, simply call the redraw() function we defined above.


In [5]:
lma_ctrl.bounds.stations=(7,99)
redraw()

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.


In [8]:
charge_lasso.charge=1
panels.lasso()

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 [16]:
chg = d.data['charge']
wh = np.where(chg > 0)
print d.data[wh]['time']


[]

Color by...

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


In [21]:
scatter_ctrl.color_field = 'chi2'
redraw()

In [7]:
scatter_ctrl.color_field = 'time'
redraw()

In [22]:
# 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'

redraw()

Animation

Animate the current view


In [18]:
n_seconds = 5
anim=scatter_ctrl.animate(n_seconds, repeat=False)


78 of 134 flashes have > 10 points. Their average area =  30.5 km^2

In [9]:
anim.repeat = False # stops the animation if repeat=True above

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 [6]:
current_events_flashes = lma_ctrl.flash_stats_for_dataset(d, scatter_ctrl.branchpoint)


1 of 3 flashes have > 10 points. Their average area =  26.4 km^2
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2

In [6]:

Flash volume for current sources


In [7]:
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
            stormdrain.support.coords.systems.GeographicSystem instance
        """
        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 
        # https://github.com/scipy/scipy/blob/master/scipy/spatial/tests/test_qhull.py#L106
        # Except the formula needs to be divided by ndim! to get the volume, cf., 
        # http://en.wikipedia.org/wiki/Simplex#Geometric_properties
        # 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
        self.volume=np.sum(np.abs(simplex_volumes))
        
    
    @coroutine
    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
            self._hull_volume()
            print "Volume of hull of points in current view is {0:5.1f}".format(
                        self.volume / 1.0e9) # (1000 m)^3




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

In [9]:
current_events_flashes.targets.add(stat_maker)

In [10]:
print current_events_flashes.targets


set([<generator object flash_stat_printer at 0x21140c88>, <generator object events_flashes_receiver at 0x200f7cd8>])
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 1 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  24.5
1 of 3 flashes have > 10 points. Their average area =  26.4 km^2
Volume of hull of points in current view is  26.3
123 of 262 flashes have > 10 points. Their average area =  22.7 km^2
Volume of hull of points in current view is 10875.9
123 of 262 flashes have > 10 points. Their average area =  22.7 km^2
Volume of hull of points in current view is 10875.9
123 of 262 flashes have > 10 points. Their average area =  22.7 km^2
Volume of hull of points in current view is 10875.9
61 of 111 flashes have > 10 points. Their average area =  31.0 km^2
Volume of hull of points in current view is 1284.7
8 of 9 flashes have > 10 points. Their average area =  29.8 km^2
Volume of hull of points in current view is 412.6
9 of 10 flashes have > 10 points. Their average area =  27.6 km^2
Volume of hull of points in current view is 703.1
8 of 9 flashes have > 10 points. Their average area =  29.8 km^2
Volume of hull of points in current view is 412.6
8 of 9 flashes have > 10 points. Their average area =  29.8 km^2
Volume of hull of points in current view is 412.6
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
2 of 2 flashes have > 10 points. Their average area =  84.2 km^2
Volume of hull of points in current view is 217.3
2 of 2 flashes have > 10 points. Their average area =  84.2 km^2
Volume of hull of points in current view is 217.3
2 of 2 flashes have > 10 points. Their average area =  84.2 km^2
Volume of hull of points in current view is 217.3
2 of 2 flashes have > 10 points. Their average area =  84.2 km^2
Volume of hull of points in current view is 217.3
2 of 2 flashes have > 10 points. Their average area =  84.2 km^2
Volume of hull of points in current view is 217.3
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 107.1
8 of 9 flashes have > 10 points. Their average area =  29.8 km^2
Volume of hull of points in current view is 412.6
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 177.2
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is   1.9
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is   1.0
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is   1.9
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is 177.2
1 of 1 flashes have > 10 points. Their average area = 158.3 km^2
Volume of hull of points in current view is  89.9

What if we take the points on the volumetric hull, and PCA decompose those? Will the envelope have a major axis that is more aligned along the plate area?

Further analysis

If your analysis is hard to explain, maybe it would help to include an equation with LaTeX: $$y=x^2$$


In [ ]: