diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/density_to_files.py b/density_to_files.py
new file mode 100644
index 0000000..f1b2cd2
--- /dev/null
+++ b/density_to_files.py
@@ -0,0 +1,269 @@
+import glob
+
+import numpy as np
+
+from density_tools import unique_vectors
+
+
+
+def coroutine(func):
+    def start(*args,**kwargs):
+        cr = func(*args,**kwargs)
+        cr.next()
+        return cr
+    return start
+
+@coroutine
+def broadcast(targets):
+    while True:
+        stuff = (yield)
+        for target in targets:
+            target.send(stuff)
+
+
+# class map_projector(object):
+#     def __init__(self, ctr_lat, ctr_lon, proj_name='eqc'):
+#         self.mapProj = MapProjection(projection=proj_name, ctrLat=ctr_lat, ctrLon=ctr_lon, lat_ts=ctr_lat, lon_0=ctr_lon)
+#         self.geoProj = GeographicSystem()
+# 
+#     def __call__(self, lon, lat, alt):
+#         x,y,z = self.mapProj.fromECEF( 
+#                 *self.geoProj.toECEF(lon, lat, alt)
+#                 )
+#         return x, y, z
+#         
+# @coroutine
+# def map_projector(ctr_lat, ctr_lon, target, proj_name='eqc'):
+#     mapProj = MapProjection(projection=proj_name, ctrLat=ctr_lat, ctrLon=ctr_lon, lat_ts=ctr_lat, lon_0=ctr_lon)
+#     geoProj = GeographicSystem()
+#     while True:
+#         lon, lat, alt = (yield)
+#         x,y,z = self.mapProj.fromECEF( 
+#                 *self.geoProj.toECEF(lon, lat, alt)
+#                 )
+#         target.send((x,y,z))
+
+
+
+
+@coroutine
+def filter_flash(target, min_points=10):
+    """ Filters flash by minimum number of points.
+    """
+    while True:
+        evs, flash = (yield)           # Receive a flash
+        if (flash['n_points'] >= 10):
+            target.send((evs, flash))
+            
+
+@coroutine
+def flashes_to_frames(time_edges, targets, time_key='start'):
+    assert len(time_edges) == (len(targets)+1)
+    while True:
+        events, flashes = (yield)
+        start_times = flashes[time_key]
+        sort_idx = np.argsort(start_times) #, order=[time_key])
+        idx = np.searchsorted(start_times[sort_idx], time_edges)
+        slices = [slice(*i) for i in zip(idx[0:-1], idx[1:])]
+        for target, s in zip(targets, slices):
+            target.send((events, flashes[sort_idx][s]))
+
+def event_yielder(evs, fls):
+    for fl in fls:
+        these_events = evs[evs['flash_id'] == fl['flash_id']]
+        # if len(these_events) <> fl['n_points']:
+            # print 'not giving all ', fl['n_points'], ' events? ', these_events.shape
+        for an_ev in these_events:
+            yield an_ev
+
+
+@coroutine
+def extract_events_for_flashes(target, flashID_key='flash_id'):
+    """ Takes a large table of events and grabs only the events belonging to the flashes.
+    """
+    
+    while True:
+        evs, fls = (yield)
+        # print 'extracting events'
+        event_dtype = evs[0].dtype
+        
+        events = np.fromiter( (event_yielder(evs, fls)) , dtype=event_dtype)
+        
+        # The line below (maybe maybe maybe)
+        # events = np.fromiter((evs[evs['flash_id'] == fl['flash_id']] for fl in fls), dtype=event_dtype)
+        # does the same thing as the two following lines, but is 10x slower.
+        # The 'mapper' could actually be optimized further by calculating it globally, once per events table,
+        # but this is fast enough and saves having to pass through another variable.
+        # mapper = dict(zip(evs['flash_id'],evs))
+        # events = np.fromiter( (mapper[fl['flash_id']] for fl in fls), dtype=event_dtype)
+
+        target.send((events, fls))
+
+# @coroutine
+# def extract_events(target, flashID_key='flash_id'):
+#     """ Takes a large table of events and grabs only the events belonging to the flash.
+#         This is useful if you filter out a bunch of flashes before going to the trouble of
+#         reading the flashes in.
+#     """
+#     while True:
+#         evs, flash = (yield)
+#         flash_id = flash[flashID_key]
+#         event_dtype = evs[0].dtype
+#         # events = [ev[:] for ev in evs if ev[flashID_key] == flash_id]
+#         # events = np.asarray(events, dtype=event_dtype)
+#         # events = evs[:]
+#         events = evs[evs[flashID_key]==flash_id]
+#         # events = np.fromiter((ev[:] for ev in evs if ev[flashID_key] == flash_id), dtype=event_dtype)
+#         target.send((events, flash))
+
+@coroutine
+def no_projection(x_coord, y_coord, z_coord, target, use_flashes=False):
+    while True:
+        events, flashes = (yield)
+        if use_flashes==True:
+            points = flashes
+        else:
+            points = events
+        x,y,z = points[x_coord], points[y_coord], points[z_coord]
+        target.send((events, flashes, x,y,z))
+    
+
+@coroutine
+def project(x_coord, y_coord, z_coord, mapProj, geoProj, target, use_flashes=False):
+    """ Adds projected coordinates to the flash and events stream"""
+    while True:
+        events, flashes = (yield)
+        if use_flashes==True:
+            points = flashes
+        else:
+            points = events
+        x,y,z = mapProj.fromECEF( 
+                *geoProj.toECEF(points[x_coord], points[y_coord], points[z_coord])
+                )
+        target.send((events, flashes, x,y,z))
+
+@coroutine
+def footprint_mean(flash_id_key='flash_id', area_key='area'):
+    """ Takes x, y, z flash locations and gets 
+        Extent density unique pixels, average all flashes 
+    """
+    while True:
+        events, flash, x,y,z = (yield)
+        # print 'Doing extent density',
+        x_i = np.floor( (x-x0)/dx ).astype('int32')
+        y_i = np.floor( (y-y0)/dy ).astype('int32')
+
+        if len(x_i) > 0:
+            footprints = dict(zip(flash[flash_id_key], flash[area_key]))
+            # print 'with points numbering', len(x_i)
+            unq_idx = unique_vectors(x_i, y_i, events['flash_id'])
+            # if x[unq_idx].shape[0] > 1:
+            fl_id = events['flash_id'][unq_idx]
+            areas = [footprints[fi] for fi in fl_id] #puts areas in same order as x[unq_idx], y[unq_idx]
+            # counts normalized by areas 
+            target.send((x[unq_idx],y[unq_idx],areas))
+        # else:
+            # print ''
+
+@coroutine
+def point_density(target):
+    """ Sends event x, y, z location directly
+    """
+    while True:
+        events, flash, x, y, z = (yield)
+        # print 'Doing point density',
+        if len(x) > 0:
+            print 'with points numbering', len(x)
+            target.send((x, y, None))
+        # else:
+            # print ''
+
+@coroutine
+def extent_density(x0, y0, dx, dy, target, flash_id_key='flash_id', weight_key=None):
+    """ This function assumes a regular grid in x and y with spacing dx, dy
+        
+        x0, y0 is the x coordinate of the lower left corner of the lower-left grid cell, 
+        i.e., the lower left node of the grid mesh in cartesian space
+        
+        Eliminates duplicate points in gridded space and sends the reduced
+        set of points to the target.
+    """
+    while True:
+        # assumes x,y,z are in same order as events
+        events, flash, x,y,z = (yield)
+        # print 'Doing extent density',
+        x_i = np.floor( (x-x0)/dx ).astype('int32')
+        y_i = np.floor( (y-y0)/dy ).astype('int32')
+
+        if len(x_i) > 0:
+            print 'extent with points numbering', len(x_i), ' with weights', weight_key
+            unq_idx = unique_vectors(x_i, y_i, events[flash_id_key])
+            # if x[unq_idx].shape[0] > 1:
+            if weight_key <> None:
+                weight_lookup = dict(zip(flash[flash_id_key], flash[weight_key]))
+                weights = [weight_lookup[fi] for fi in events[unq_idx]['flash_id']] #puts weights in same order as x[unq_idx], y[unq_idx]
+            else:
+                weights = None
+            target.send((x[unq_idx], y[unq_idx], weights))
+        # else:
+            # print ''
+
+
+@coroutine
+def accumulate_points_on_grid(grid, xedge, yedge, out=None, label=''):
+    assert xedge.shape[0] == grid.shape[0]+1
+    assert yedge.shape[0] == grid.shape[1]+1
+    if out == None:
+        out = {}
+    # grid = None
+    try:
+        while True:
+            x, y, weights = (yield)
+            if len(x) > 0:
+                x = np.atleast_1d(x)
+                y = np.atleast_1d(y)
+                print 'accumulating ', len(x), 'points for ', label
+                count, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=None, normed=False)    
+                
+                if weights <> None:
+                    # histogramdd sums up weights in each bin for normed=False
+                    total, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=weights, normed=False)
+                    # return the mean of the weights in each bin
+                    bad = (count <= 0)
+                    count = np.asarray(total, dtype='float32')/count
+                    count[bad] = 0.0 
+                    
+                # try:
+                    # count, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=weights)
+                # except AttributeError:
+                #     # if x,y are each scalars, need to make 1D arrays
+                #     x = np.asarray((x,))
+                #     y = np.asarray((y,))
+                #     count, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=weights)
+                # using += (as opposed to grid = grid + count) is essential
+                # so that the user can retain a reference to the grid object
+                # outside this routine.
+                if grid == None:
+                    grid = count
+                    out['out'] = grid
+                else:
+                    grid += count
+    except GeneratorExit:
+        out['out'] = grid
+        
+        
+    
+        
+# if __name__ == '__main__':
+#     do_profile=False
+#     if do_profile:
+#         import hotshot
+#         from hotshot import stats
+#         prof = hotshot.Profile("density_test_profile")
+#         prof.runcall(example)
+#         prof.close()
+#         s=stats.load("density_test_profile")
+#         s.sort_stats("time").print_stats()
+#     else:
+#         x_coord, y_coord, lons, lats, test_grid = example()
+
diff --git a/density_tools.py b/density_tools.py
new file mode 100644
index 0000000..58f29c0
--- /dev/null
+++ b/density_tools.py
@@ -0,0 +1,118 @@
+import numpy as np
+
+def unique_vectors(x_i, y_i, g_id):
+    """ x_i, y_i: Discretized (bin index) values of point locations.
+        id:       Entity ids that group points
+        
+        Returns a (3, n_points) array of unique points and ids.
+        """
+        
+    # Because we do casting below, ascontiguousarray is important
+    locs = np.ascontiguousarray(np.vstack( (
+                           x_i.astype('int32'), 
+                           y_i.astype('int32'), 
+                           g_id.astype('int32')
+                       ) ).T )
+
+	# Find unique rows (unique grid cells occupied per flash) by casting to a set
+	# of strings comprised of the bytes that make up the rows. Strings are not serialized, just ints cast to byte strings.
+	# Based on the example at
+	# http://www.mail-archive.com/numpy-discussion@scipy.org/msg04176.html
+    # which doesn't quite get it right: it returns unique elements.
+    unq, unq_idx = np.unique(locs.view(
+                    ('S%d'%(locs.itemsize*locs.shape[1]))*locs.shape[0]), 
+                    return_index=True)
+    # these are unnecessary since we only need the row indices
+    # unq = unq.view('int32')
+    # unq = unq.reshape((len(unq)/locs.shape[1], locs.shape[1]))
+    
+    return unq_idx
+
+def extent_density(x, y, ids, x0, y0, dx, dy, xedge, yedge):
+    x_i = np.floor( (x-x0)/dx ).astype('int32')
+    y_i = np.floor( (y-y0)/dy ).astype('int32')
+    
+    unq_idx = unique_vectors(x_i, y_i, ids)
+    # if x[unq_idx].shape[0] > 1:
+    density,edges = np.histogramdd((x[unq_idx],y[unq_idx]), bins=(xedge,yedge))
+    return density, edges
+    
+
+
+def test_extent_density():
+    # create a set of x,y,id, that have one point in each grid location defined by x0,x1,dx
+    x = np.asarray((0.5, 1.5, 2.5, 3.5, 4.5))
+    x0, x1 = 0.0, 5.0
+    dx = 1.0
+    y = x - 2
+    y0, y1 = x0 - 2, x1 - 2
+    dy = dx
+    x_cover, y_cover = np.meshgrid(x,y)
+    xedge = np.arange(x0, x1+dx, dx)
+    yedge = np.arange(y0, y1+dy, dy)
+    ids = np.ones(y_cover.flatten().shape[0], dtype=int)
+    
+    # ------------------------------
+    # replicate the x and y points to have two points in each grid location
+    x_doubled = np.hstack((x_cover.flatten(), x_cover.flatten()))
+    y_doubled = np.hstack((y_cover.flatten(), y_cover.flatten()))
+
+    # replicate the ids such that the doubled points in each grid belong to the same entity
+    ids_doubled = np.hstack((ids, ids))
+    density, edges = extent_density(x_doubled, y_doubled, ids_doubled, 
+                                    x0, y0, dx, dy, xedge, yedge)
+    assert (density == 1).all()
+
+    # replicate the ids such that the doubled points in each grid belong to different entities
+    ids_doubled = np.hstack((ids, ids+1))
+    density, edges = extent_density(x_doubled, y_doubled, ids_doubled, 
+                                    x0, y0, dx, dy, xedge, yedge)
+    assert (density == 2).all()
+    
+    # ------------------------------
+    # replicate the x and y points to have two points in each grid location, but leave out
+    # one of the points (0.5, -1.5); lower-left in space, upper left in the density array printout
+    x_doubled = np.hstack((x_cover.flatten()[1:], x_cover.flatten()))
+    y_doubled = np.hstack((y_cover.flatten()[1:], y_cover.flatten()))
+
+    # replicate the ids such that the doubled points in each grid belong to the different entities
+    ids_doubled = np.hstack((ids[1:], ids+1))
+    density, edges = extent_density(x_doubled, y_doubled, ids_doubled, 
+                                    x0, y0, dx, dy, xedge, yedge)
+    assert (density == np.array([[ 1.,  2.,  2.,  2.,  2.],
+                                [ 2.,  2.,  2.,  2.,  2.],
+                                [ 2.,  2.,  2.,  2.,  2.],
+                                [ 2.,  2.,  2.,  2.,  2.],
+                                [ 2.,  2.,  2.,  2.,  2.]] ) ).all()
+    
+    # replicate the ids such that the doubled points in each grid belong to the same entity
+    ids_doubled = np.hstack((ids[1:], ids))
+    density, edges = extent_density(x_doubled, y_doubled, ids_doubled, 
+                                    x0, y0, dx, dy, xedge, yedge)
+    assert (density == 1).all()
+
+
+def test_unq():
+     locs = np.array([[ 0,  1,  2],
+                      [ 3,  4,  5],
+                      [ 8, 10, 11],
+                      [ 3,  5,  5],
+                      [ 3,  4,  5],
+                      [ 3,  4,  6],
+                      [ 9, 10, 11]])
+     
+     # this returns the unique elements
+     # unq, unq_idx = np.unique(locs.view('S%d'%locs.itemsize*locs.shape[0]), return_index=True)
+     
+     unq, unq_idx = np.unique(locs.view(('S%d'%(locs.itemsize*locs.shape[1]))*locs.shape[0]), return_index=True)
+     unq = unq.view('int32')
+     unq = unq.reshape((len(unq)/locs.shape[1], locs.shape[1]))
+     # print unq
+     # print locs[unq_idx] == unq
+     assert (locs[unq_idx] == unq).all()
+
+        
+if __name__ == '__main__':
+    test_unq()
+    test_extent_density()
+    print "Tests complete."
diff --git a/examples/sample_make_grid.py b/examples/sample_make_grid.py
new file mode 100644
index 0000000..d094624
--- /dev/null
+++ b/examples/sample_make_grid.py
@@ -0,0 +1,27 @@
+import glob
+from datetime import datetime
+
+from LMAtools.make_grids import grid_h5flashfiles
+
+
+h5_filenames = glob.glob('LYL*.flash.h5')
+start_time = datetime(2009,7,24, 1,00,00)
+end_time   = datetime(2009,7,24, 1,04,00)
+
+frame_interval=60.0
+dx=8.0e3
+dy=8.0e3
+x_bnd = (-100e3, 100e3)
+y_bnd = (-100e3, 100e3)
+
+# # KOUN
+# ctr_lat = 35.23833
+# ctr_lon = -97.46028
+
+# DC
+ctr_lat =  38.889444
+ctr_lon =  -77.035278
+
+
+grid_h5flashfiles(h5_filenames, start_time, end_time, frame_interval=frame_interval, 
+            dx=dx, dy=dy, x_bnd=x_bnd, y_bnd=y_bnd, ctr_lon=ctr_lon, ctr_lat=ctr_lat)
\ No newline at end of file
diff --git a/examples/sample_plot.py b/examples/sample_plot.py
new file mode 100644
index 0000000..cbecc21
--- /dev/null
+++ b/examples/sample_plot.py
@@ -0,0 +1,8 @@
+from LMAtools.multiples_nc import make_plot
+
+n_cols=2
+
+make_plot('LMA_20090724_010000_240_flash_extent.nc', 'flash_extent', n_cols=n_cols)
+make_plot('LMA_20090724_010000_240_flash_init.nc', 'flash_initiation', n_cols=n_cols)
+make_plot('LMA_20090724_010000_240_source.nc', 'lma_source', n_cols=n_cols)
+make_plot('LMA_20090724_010000_240_footprint.nc', 'flash_footprint', n_cols=n_cols)
\ No newline at end of file
diff --git a/make_grids.py b/make_grids.py
new file mode 100644
index 0000000..8e7ec54
--- /dev/null
+++ b/make_grids.py
@@ -0,0 +1,282 @@
+import glob 
+from datetime import datetime, timedelta
+
+import numpy as np
+import tables
+
+import density_to_files
+
+from acuity.coordinateSystems import MapProjection, GeographicSystem
+
+
+def to_seconds(dt):
+    'convert datetime.timedelta object to float seconds'
+    return dt.days * 86400 + dt.seconds + dt.microseconds/1.0e6
+
+
+def write_cf_netcdf(outfile, t_start, t, xloc, yloc, lon_for_x, lat_for_y, ctr_lat, ctr_lon, grid, grid_var_name, grid_description, format='i'):
+    """ Write a Climate and Forecast Metadata-compliant NetCDF file. 
+    
+        Should display natively in conformant packages like McIDAS-V.
+        
+    """
+
+    import pupynere as nc
+
+    missing_value = -9999
+    
+    nc_out = nc.NetCDFFile(outfile, 'w')
+    nc_out.createDimension('nx', xloc.shape[0])
+    nc_out.createDimension('ny', yloc.shape[0])
+    nc_out.createDimension('ntimes', t.shape[0])  #unlimited==None
+
+    proj = nc_out.createVariable('Lambert_Azimuthal_Equal_Area', 'i', ())
+    proj.grid_mapping_name = 'lambert_azimuthal_equal_area'
+    proj.longitude_of_projection_origin = ctr_lon
+    proj.latitude_of_projection_origin = ctr_lat
+    proj.false_easting = 0.0
+    proj.false_northing = 0.0
+    
+    # x_coord = nc_out.createVariable('longitude', 'f', ('nx',))
+    # x_coord.long_name="longitude"
+    # x_coord.standard_name="longitude"
+    # x_coord.units = "degrees_east"
+    x_coord = nc_out.createVariable('x', 'f', ('nx',))
+    x_coord.units = "km"
+    x_coord.long_name = "x coordinate of projection"
+    x_coord.standard_name = 'projection_x_coordinate'
+
+    # y_coord = nc_out.createVariable('latitude', 'f', ('nx',))
+    # y_coord.long_name="latitude"
+    # y_coord.standard_name="latitude"
+    # y_coord.units = "degrees_north"
+    y_coord = nc_out.createVariable('y', 'f', ('ny',))
+    y_coord.units = "km"
+    y_coord.long_name = "y coordinate of projection"
+    y_coord.standard_name = 'projection_y_coordinate'
+
+    times = nc_out.createVariable('time', 'f', ('ntimes',) )#, filters=no_compress)
+    times.long_name="time"
+    times.units = "seconds since %s" % t_start.strftime('%Y-%m-%d %H:%M:%S')
+    
+    lons = nc_out.createVariable('lons', 'd', ('nx','ny') )#, filters=no_compress)
+    lons.long_name="longitude"
+    lons.standard_name="longitude"
+    lons.units = "degrees_east"
+    
+    lats = nc_out.createVariable('lats', 'd', ('nx','ny') )#, filters=no_compress)
+    lats.long_name="latitude"
+    lats.standard_name="latitude"
+    lats.units = "degrees_north"
+
+    lightning2d = nc_out.createVariable(grid_var_name, format, ('ntimes','nx','ny') )#, filters=no_compress)
+    lightning2d.long_name=grid_description #'LMA VHF event counts (vertically integrated)'
+    lightning2d.units='dimensionless'
+    lightning2d.coordinates='time lons lats'
+    lightning2d.grid_mapping = "Lambert_Azimuthal_Equal_Area"
+    lightning2d.missing_value = missing_value
+
+    x_coord[:] = xloc[:]
+    y_coord[:] = yloc[:]
+    times[:] = t[:]
+    lons[:] = lon_for_x[:]
+    lats[:] = lat_for_y[:]
+
+    for i in range(grid.shape[2]):
+        lightning2d[i,:,:] = grid[:,:,i]
+    nc_out.close()
+
+
+
+def read_flashes(h5LMAfiles, target, base_date=None, min_points=10):
+
+    for filename in h5LMAfiles:
+        h5 = tables.openFile(filename)
+        table_name = h5.root.events._v_children.keys()[0]
+
+        print filename
+
+        # event_dtype = getattr(h5.root.events, table_name)[0].dtype
+        # events_all = getattr(h5.root.events, table_name)[:]
+        flashes = getattr(h5.root.flashes, table_name)
+        events = getattr(h5.root.events, table_name)[:]
+
+        fl_dtype = flashes[0].dtype
+        flashes = np.fromiter((fl[:] for fl in flashes if fl['n_points'] >= min_points), dtype=fl_dtype)
+
+        # get the start date of the file, which is the seconds-of-day reference
+        h5_start =  datetime(*getattr(h5.root.events, table_name).attrs['start_time'][0:3])
+        if base_date==None:
+            base_date = h5_start
+        extra_dt = to_seconds(h5_start - base_date)
+        # adjust time in seconds of day to be in correct reference
+        events['time'] += extra_dt
+        flashes['start'] += extra_dt
+
+        n_flashes = len(flashes)
+        print '    ', n_flashes, ' total flashes, with extra dt of', extra_dt
+
+        target.send((events,flashes))
+
+        h5.close()
+        
+
+def grid_h5flashfiles(h5_filenames, start_time, end_time, 
+                        frame_interval=120.0, dx=4.0e3, dy=4.0e3,
+                        x_bnd = (-100e3, 100e3),
+                        y_bnd = (-100e3, 100e3),
+                        z_bnd = (-20e3, 20e3),
+                        ctr_lat = 35.23833, ctr_lon = -97.46028,
+                        ):
+    from math import ceil
+    """
+    
+    Create 2D plan-view density grids for events, flash origins, flash extents, and mean flash footprint
+    
+    frame_interval: Frame time-step in seconds
+    dx, dy: horizontal grid size in m
+    {x,y,z}_bnd: horizontal grid edges in m
+    ctr_lat, ctr_lon: coordinate center
+    
+    Uses an azimuthal equidistant map projection on the WGS84 ellipsoid.
+    
+    
+    read_flashes
+    
+    filter_flash
+    extract_events
+    flash_to_frame
+    
+    frame0_broadcast, frame1_broadcast, ...
+    
+    each broadcaster above sends events and flashes to:
+    projection( event_location), projection(flash_init_location), projection(event_location)
+    which map respectively to:
+    point_density->accum_on_grid(event density), point_density->accum_on_grid(flash init density), extent_density->accum_on_grid(flash_extent_density)
+
+    grids are in an HDF5 file. how to handle flushing?
+    """
+    
+    # reference time is the date part of the start_time
+    ref_date = start_time.date()
+    t_ref = datetime(ref_date.year, ref_date.month, ref_date.day)
+    
+    frame_dt = timedelta(0, frame_interval, 0)
+    
+    duration = end_time - start_time
+    n_frames = int(ceil(to_seconds(duration) / to_seconds(frame_dt)))
+
+    t_edges = [start_time + i*frame_dt for i in range(n_frames+1)]
+    
+    t_edges_seconds = [to_seconds(edge - t_ref) for edge in t_edges]
+    
+    xedge=np.arange(x_bnd[0], x_bnd[1]+dx, dx)
+    yedge=np.arange(y_bnd[0], y_bnd[1]+dy, dy)
+    
+    x0 = xedge[0]
+    y0 = yedge[0]
+    
+    mapProj = MapProjection(projection='aeqd', ctrLat=ctr_lat, ctrLon=ctr_lon, lat_ts=ctr_lat, lon_0=ctr_lon, lat_0=ctr_lat)
+    geoProj = GeographicSystem()
+    
+    
+    event_density_grid  = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='int32')
+    init_density_grid   = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='int32')
+    extent_density_grid = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='int32')
+    footprint_grid      = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='float32')
+        
+    all_frames = []
+    extent_frames = []
+    init_frames = []
+    event_frames = []
+    for i in range(n_frames):
+        extent_out = {'name':'extent'}
+        init_out   = {'name':'init'}
+        event_out  = {'name':'event'}
+        accum_event_density  = density_to_files.accumulate_points_on_grid(event_density_grid[:,:,i], xedge, yedge,  out=event_out, label='event')
+        accum_init_density   = density_to_files.accumulate_points_on_grid(init_density_grid[:,:,i], xedge, yedge,   out=init_out,  label='init')
+        accum_extent_density = density_to_files.accumulate_points_on_grid(extent_density_grid[:,:,i], xedge, yedge, out=extent_out,label='extent')
+        accum_footprint      = density_to_files.accumulate_points_on_grid(footprint_grid[:,:,i], xedge, yedge, label='footprint')
+        extent_out['func'] = accum_extent_density
+        init_out['func'] = accum_init_density
+        event_out['func'] = accum_event_density
+        extent_frames.append(extent_out)
+        init_frames.append(init_out)
+        event_frames.append(event_out)
+        
+        event_density_target  = density_to_files.point_density(accum_event_density)
+        init_density_target   = density_to_files.point_density(accum_init_density)
+        extent_density_target = density_to_files.extent_density(x0, y0, dx, dy, accum_extent_density)
+        mean_footprint_target = density_to_files.extent_density(x0, y0, dx, dy, accum_footprint, weight_key='area')
+
+        spew_to_density_types = density_to_files.broadcast( ( 
+                    density_to_files.project('lon', 'lat', 'alt', mapProj, geoProj, event_density_target, use_flashes=False),
+                    density_to_files.project('init_lon', 'init_lat', 'init_alt', mapProj, geoProj, init_density_target, use_flashes=True),
+                    density_to_files.project('lon', 'lat', 'alt', mapProj, geoProj, extent_density_target, use_flashes=False),
+                    density_to_files.project('lon', 'lat', 'alt', mapProj, geoProj, mean_footprint_target, use_flashes=False),
+                    ) )
+
+        all_frames.append( density_to_files.extract_events_for_flashes( spew_to_density_types ) )
+        
+    framer = density_to_files.flashes_to_frames(t_edges_seconds, all_frames, time_key='start')
+    
+    read_flashes( h5_filenames, framer, base_date=t_ref)
+    
+    # print 'event_density_grid ', id(event_density_grid[:,:,-1])
+    # print 'extent_density_grid', id(extent_density_grid[:,:,-1])
+    # print 'init_density_grid  ', id(init_density_grid[:,:,-1])
+    
+    
+    x_coord = (xedge[:-1] + xedge[1:])/2.0
+    y_coord = (yedge[:-1] + yedge[1:])/2.0
+    nx = x_coord.shape[0]
+    ny = y_coord.shape[0]
+    
+    x_all, y_all = (a.T for a in np.meshgrid(x_coord, y_coord))
+    assert x_all.shape == y_all.shape
+    z_all = np.zeros_like(x_all)
+            
+    lons, lats, alts = x,y,z = geoProj.fromECEF( *mapProj.toECEF(x_all, y_all, z_all) )
+    lons.shape=x_all.shape
+    lats.shape=y_all.shape
+    
+    outflile_basename = 'LMA_%s_%d_' % (start_time.strftime('%Y%m%d_%H%M%S'), to_seconds(duration))
+    
+    write_cf_netcdf(outflile_basename+'flash_extent.nc', t_ref, np.asarray(t_edges_seconds[:-1]),
+                    x_coord/1.0e3, y_coord/1.0e3, lons, lats, ctr_lat, ctr_lon, 
+                    extent_density_grid, 'flash_extent', 'LMA flash extent density')
+    write_cf_netcdf(outflile_basename+'flash_init.nc', t_ref, np.asarray(t_edges_seconds[:-1]),
+                    x_coord/1.0e3, y_coord/1.0e3, lons, lats, ctr_lat, ctr_lon, 
+                    init_density_grid, 'flash_initiation', 'LMA flash initiation density')
+    write_cf_netcdf(outflile_basename+'source.nc', t_ref, np.asarray(t_edges_seconds[:-1]),
+                    x_coord/1.0e3, y_coord/1.0e3, lons, lats, ctr_lat, ctr_lon, 
+                    event_density_grid, 'lma_source', 'LMA source density')
+    write_cf_netcdf(outflile_basename+'footprint.nc', t_ref, np.asarray(t_edges_seconds[:-1]),
+                    x_coord/1.0e3, y_coord/1.0e3, lons, lats, ctr_lat, ctr_lon, 
+                    footprint_grid, 'flash_footprint', 'LMA per-flash footprint', format='f')
+                    
+    print 'max extent is', extent_density_grid.max()
+
+    return x_coord, y_coord, lons, lats, extent_density_grid
+
+
+if __name__ == '__main__':
+    h5_filenames = glob.glob('data/LYL*090610_20*.h5')
+    start_time = datetime(2009,6,10, 20,0,0)
+    end_time   = datetime(2009,6,10, 21,0,0)
+    
+    frame_interval=120.0
+    dx=4.0e3
+    dy=4.0e3
+    x_bnd = (-100e3, 100e3)
+    y_bnd = (-100e3, 100e3)
+    # # KOUN
+    ctr_lat = 35.23833
+    ctr_lon = -97.46028
+
+    # DC
+    # ctr_lat =  38.889444
+    # ctr_lon =  -77.035278
+    
+    grid_h5flashfiles(h5_filenames, start_time, end_time, frame_interval=frame_interval, 
+                dx=dx, dy=dy, x_bnd=x_bnd, y_bnd=y_bnd, ctr_lon=ctr_lon, ctr_lat=ctr_lat)
\ No newline at end of file
diff --git a/multiples.py b/multiples.py
new file mode 100644
index 0000000..16c258b
--- /dev/null
+++ b/multiples.py
@@ -0,0 +1,293 @@
+from small_multiples import small_multiples_plot
+from acuity.data import Bounds, Data, DataTransform, DataSelection
+from acuity.coordinateSystems import MapProjection, GeographicSystem
+from acuity.LMA.LMAdataHDF import LMAdataManagerHDF
+# from acuity.LMA.LMAarrayFile import LMAdataFile
+from acuity.LMA.LMAdata import LMAdataManager
+from acuity.views import AcuityView
+
+from density_tools import extent_density
+
+from mx.DateTime import DateTime, DateTimeDelta
+import numpy as np
+from pylab import figure, get_cmap, colorbar
+from matplotlib.figure import figaspect
+from matplotlib.colorbar import ColorbarBase
+from matplotlib.ticker import FuncFormatter
+from matplotlib.dates import mx2num, date2num, DateFormatter
+
+from math import ceil
+
+import pytz
+tz=pytz.timezone('US/Eastern') # Why, oh, why, is it using my local time zone?
+time_series_x_fmt = DateFormatter('%H%M:%S', tz=tz)
+
+
+def kilo(x, pos):
+    'The two args are the value and tick position'
+    return '%s' % (x/1000.0)
+
+kilo_formatter = FuncFormatter(kilo)
+
+
+loadLMA = True
+
+# /* DC LMA */
+DClat = 38.8888500 # falls church / western tip of arlington, rough centroid of stations
+DClon = -77.1685800
+kounLat = 35.23833
+kounLon = -97.46028
+kounAlt = 377.0
+radarLat=kounLat
+radarLon=kounLon
+radarAlt=0.0
+
+# # ARMOR
+# radarLat = 34.6461
+# radarLon = -86.7714
+# radarAlt = 180.0
+
+mapProj = MapProjection(projection='eqc', ctrLat=radarLat, ctrLon=radarLon, lat_ts=radarLat, lon_0=radarLon)
+geoProj = GeographicSystem()
+
+dx, dy = (8.0e3,)*2
+# dx, dy = 500.0, 500.0
+minute_intervals = [2.0]
+count_scale_factor = dx / 1000.0
+
+b = Bounds()
+#Crude limits to the domain of the storm
+b.x = (-60e3, 140e3)
+b.y = (-150e3, 50e3)
+b.z = (-20e3, 20e3)
+b.chi2 = (0,1.0)
+b.stations = (7,99)
+start_time = DateTime(2009,6,10,20,50,0)
+end_time   = DateTime(2009,6,10,21,00,0)
+max_count_baseline = 450 * count_scale_factor #/ 10.0
+
+pad = 0 #-25e3
+
+# approximate velocity of reference frame. used to adjust the viewport.
+# is average speed of LMA density center between 830 and 945 UTC
+u = 0 #17.8    # m/s
+v = 0 #15.6 
+view_dx = b.x[1]-b.x[0] #200.0e3
+view_dy = b.y[1]-b.y[0] #200.0e3
+# Position at some initial time
+x0 = b.x[1]-view_dx/2.0#-150.0e3
+y0 = b.y[1]-view_dy/2.0#-150.0e3
+t0 = DateTime(2009,6,10,22,40,0)
+
+source_density=False
+
+import glob
+
+if source_density==True:
+    LMAfiles = glob.glob("data/LYL*090610_20*.dat.gz") #+ glob.glob("data/LYL*090610_21*.dat.gz")
+    lmaManager = LMAdataManager(LMAfiles)
+    lma_view = AcuityView(DataSelection(lmaManager.data, b), mapProj, bounds=b)    
+else:
+    import tables
+    LMAfilesHDF = glob.glob('data/LYL*090610_20*.h5')
+    LMAtables = []
+    for hdffile in LMAfilesHDF:
+        h5 = tables.openFile(hdffile)
+        table_name = h5.root.events._v_children.keys()[0]
+        LMAtables.append('/events/'+table_name)
+    
+        # events = getattr(h5.root.events, table_name)[:]
+        # flashes = getattr(h5.root.flashes, table_name)[:]
+        # mapping = dict( ( fl, events[events['flash_id'] == fl['flash_id']] ) 
+        #                 for fl in flashes if (fl['n_points']>9)
+        #               )
+        h5.close()
+    HDFmanagers = [LMAdataManagerHDF(*args) for args in zip(LMAfilesHDF, LMAtables)]
+    # LMAtables   = [ '/events/LMA_080206_080000_3600', '/events/LMA_080206_090000_3600', '/events/LMA_080206_100000_3600' ]
+
+    # h5 = tables.openFile('data/LYLOUT_090610_180000_0600.dat.gz.flash.h5')
+    # table_name = h5.root.events._v_children.keys()[0]
+    # events = getattr(h5.root.events, table_name)[:]
+    # flashes = getattr(h5.root.flashes, table_name)[:]
+    # mapping = dict((fl, events[events['flash_id'] == fl['flash_id']]) for fl in flashes)
+    # h5.close()
+
+
+
+def runtest(lmaManager=None, lma_view=None, HDFmanagers=None):
+    # colormap = get_cmap('gist_yarg_r')
+    colormap = get_cmap('gist_earth')
+    
+    density_maxes = []
+    total_counts = []
+    all_t = []
+    
+    for delta_minutes in minute_intervals:
+        time_delta = DateTimeDelta(0, 0, delta_minutes, 0)
+        
+        n_frames   = int(ceil((end_time - start_time) / time_delta))
+        n_cols = 6
+        n_rows = int(ceil( float(n_frames) / n_cols ))
+        w, h = figaspect(float(n_rows)/n_cols)
+
+        xedge=np.arange(b.x[0], b.x[1]+dx, dx)
+        yedge=np.arange(b.y[0], b.y[1]+dy, dy)
+        x_range = b.x[1] - b.x[0]
+        y_range = b.y[1] - b.y[0]
+
+        min_count, max_count = 1, max_count_baseline*delta_minutes
+
+        f = figure(figsize=(w,h))
+        p = small_multiples_plot(fig=f, rows=n_rows, columns=n_cols)
+        p.label_edges(True)
+        
+        for ax in p.multiples.flat:
+            ax.yaxis.set_major_formatter(kilo_formatter)
+            ax.xaxis.set_major_formatter(kilo_formatter)
+
+        for i in range(n_frames):
+            frame_start = start_time + i*time_delta
+            frame_end   = frame_start + time_delta
+            b.sec_of_day = (frame_start.abstime, frame_end.abstime)
+            b.t = (frame_start, frame_end)
+            
+            do_plot = False
+            flash_extent_density = True
+            density = None
+            
+            if source_density==True:
+                lmaManager.refresh(b)
+                lma_view.transformed.cache_is_old()
+                x,y,t=lma_view.transformed['x','y','t']
+                density,edges = np.histogramdd((x,y), bins=(xedge,yedge))
+                do_plot=True
+            else:
+                for lmaManager in HDFmanagers:
+                    # yes, loop through every file every time and reselect data.
+                    # so wrong, yet so convenient.
+                    h5 = lmaManager.h5file
+                    if flash_extent_density == False:
+                        lmaManager.refresh(b)
+                        lma_view = AcuityView(DataSelection(lmaManager.data, b), mapProj, bounds=b)
+                        # lma_view.transformed.cache_is_old()
+                        x,y,t=lma_view.transformed['x','y','t']
+                        if x.shape[0] > 1: do_plot = True
+                        break
+                    else:
+                        # assume here that the bounds sec_of_day day is the same as
+                        # the dataset day
+                        t0, t1 = b.sec_of_day
+                        # events = getattr(h5.root.events, lmaManager.table.name)[:]
+                        # flashes = getattr(h5.root.flashes, lmaManager.table.name)[:]
+                        
+                        event_dtype = getattr(h5.root.events, lmaManager.table.name)[0].dtype
+                        events_all = getattr(h5.root.events, lmaManager.table.name)[:]
+                        flashes = getattr(h5.root.flashes, lmaManager.table.name)
+                        
+                        def event_yielder(evs, fls):
+                            these_events = []
+                            for fl in fls:
+                                if (    (fl['n_points']>9) & 
+                                        (t0 < fl['start']) & 
+                                        (fl['start'] <= t1) 
+                                    ):
+                                    these_events = evs[evs['flash_id'] == fl['flash_id']]
+                                    if len(these_events) <> fl['n_points']:
+                                        print 'not giving all ', fl['n_points'], ' events? ', these_events.shape
+                                    for an_ev in these_events:
+                                        yield an_ev
+
+                        
+                        # events = np.fromiter((an_ev for an_ev in ( events_all[events_all['flash_id'] == fl['flash_id']] 
+                        #                 for fl in flashes if (
+                        #                   (fl['n_points']>9) & (t0 < fl['start']) & (fl['start'] <= t1)
+                        #                 )
+                        #               ) ), dtype=event_dtype)
+                        events = np.fromiter(event_yielder(events_all, flashes), dtype=event_dtype)
+                        
+                        # print events['flash_id'].shape
+
+                        ### Flash extent density ###                        
+                        x,y,z = mapProj.fromECEF( 
+                                *geoProj.toECEF(events['lon'], events['lat'], events['alt'])
+                                )
+                                
+                        # Convert to integer grid coordinate bins
+                        #      0    1    2    3
+                        #   |    |    |    |    |
+                        # -1.5  0.0  1.5  3.0  4.5
+                    
+                        if x.shape[0] > 1:
+                            density, edges = extent_density(x,y,events['flash_id'].astype('int32'),
+                                                            b.x[0], b.y[0], dx, dy, xedge, yedge)
+                            do_plot = True                        
+                            break
+                # print 'density values: ', density.min(), density.max()
+                    
+            
+            if do_plot == True:  # need some data
+                # density,edges = np.histogramdd((x,y), bins=(xedge,yedge))
+                density_plot  = p.multiples.flat[i].pcolormesh(xedge,yedge,
+                                           np.log10(density.transpose()), 
+                                           vmin=-0.2,
+                                           vmax=np.log10(max_count),
+                                           cmap=colormap)
+                label_string = frame_start.strftime('%H%M:%S')
+                text_label = p.multiples.flat[i].text(b.x[0]-pad+x_range*.01, b.y[0]-pad+y_range*.01, label_string, color=(0.5,)*3, size=6)
+                density_plot.set_rasterized(True)
+                density_maxes.append(density.max())
+                total_counts.append(density.sum())
+                all_t.append(frame_start)
+                print label_string, x.shape, density.max(), density.sum()
+
+        color_scale = ColorbarBase(p.colorbar_ax, cmap=density_plot.cmap,
+                                           norm=density_plot.norm,
+                                           orientation='horizontal')
+        # color_scale.set_label('count per pixel')
+        color_scale.set_label('log10(count per pixel)')
+        
+        # moving reference frame correction. all panels will have same limits, based on time of last frame
+        view_dt = 0.0 # (frame_start - t0).seconds
+        x_ctr = x0 + view_dt*u
+        y_ctr = y0 + view_dt*v
+        view_x = (x_ctr - view_dx/2.0 - pad, x_ctr + view_dx/2.0 + pad)
+        view_y = (y_ctr - view_dy/2.0 - pad, y_ctr + view_dy/2.0 + pad)
+        # view_x  = (b.x[0]+view_dt*u, b.x[1]+view_dt*u)
+        # view_y  = (b.y[0]+view_dt*v, b.y[1]+view_dt*v)
+        
+        # print 'making timeseries',
+        # time_series = figure(figsize=(16,9))
+        # ts_ax = time_series.add_subplot(111)
+        # ts_ax.plot_date(mx2num(all_t),total_counts,'-', label='total sources', tz=tz)
+        # ts_ax.plot_date(mx2num(all_t),density_maxes,'-', label='max pixel', tz=tz)
+        # ts_ax.xaxis.set_major_formatter(time_series_x_fmt)
+        # ts_ax.legend()
+        # time_filename = 'out/LMA-timeseries_%s_%5.2fkm_%5.1fs.pdf' % (start_time.strftime('%Y%m%d_%H%M%S'), dx/1000.0, time_delta.seconds)
+        # time_series.savefig(time_filename)
+        # print ' ... done'
+        
+        print 'making multiples',
+        p.multiples.flat[0].axis(view_x+view_y)
+        filename = 'out/LMA-density_%s_%5.2fkm_%5.1fs.pdf' % (start_time.strftime('%Y%m%d_%H%M%S'), dx/1000.0, time_delta.seconds)
+        f.savefig(filename, dpi=150)
+        print ' ... done'
+        f.clf()
+        return events
+        
+
+        
+if __name__ == '__main__':
+    do_profile=True
+    if do_profile:
+        import hotshot
+        from hotshot import stats
+        prof = hotshot.Profile("multiples_test_profile")
+        prof.runcall(runtest, HDFmanagers=HDFmanagers)
+        prof.close()
+        s=stats.load("multiples_test_profile")
+        s.sort_stats("time").print_stats()
+    else:
+        if source_density:
+            res = runtest(lmaManager=lmaManager, lma_view=lma_view)
+        else:
+            res = runtest(HDFmanagers=HDFmanagers)
\ No newline at end of file
diff --git a/multiples_nc.py b/multiples_nc.py
new file mode 100644
index 0000000..78c0a32
--- /dev/null
+++ b/multiples_nc.py
@@ -0,0 +1,171 @@
+from small_multiples import small_multiples_plot
+# from acuity.data import Bounds, Data, DataTransform, DataSelection
+# from acuity.coordinateSystems import MapProjection, GeographicSystem
+
+from datetime import datetime, timedelta
+import numpy as np
+import pupynere as nc
+
+from pylab import figure, get_cmap, colorbar
+from matplotlib.figure import figaspect
+from matplotlib.colorbar import ColorbarBase
+from matplotlib.ticker import FuncFormatter
+from matplotlib.dates import mx2num, date2num, DateFormatter
+
+from math import ceil
+
+import pytz
+tz=pytz.timezone('US/Eastern') # Why, oh, why, is it using my local time zone?
+time_series_x_fmt = DateFormatter('%H%M:%S', tz=tz)
+
+
+def kilo(x, pos):
+    'The two args are the value and tick position'
+    return '%s' % (x/1000.0)
+
+kilo_formatter = FuncFormatter(kilo)
+
+# # /* DC LMA */
+# DClat = 38.8888500 # falls church / western tip of arlington, rough centroid of stations
+# DClon = -77.1685800
+# kounLat = 35.23833
+# kounLon = -97.46028
+# kounAlt = 377.0
+# radarLat=kounLat
+# radarLon=kounLon
+# radarAlt=0.0
+
+# # ARMOR
+# radarLat = 34.6461
+# radarLon = -86.7714
+# radarAlt = 180.0
+
+# mapProj = MapProjection(projection='eqc', ctrLat=radarLat, ctrLon=radarLon, lat_ts=radarLat, lon_0=radarLon)
+# geoProj = GeographicSystem()
+
+def centers_to_edges(x):
+    xedge=np.zeros(x.shape[0]+1)
+    xedge[1:-1] = (x[:-1] + x[1:])/2.0
+    dx = np.mean(np.abs(xedge[2:-1] - xedge[1:-2]))
+    xedge[0] = xedge[1] - dx
+    xedge[-1] = xedge[-2] + dx
+    return xedge        
+
+
+
+def make_plot(filename, grid_name, x_name='x', y_name='y', t_name='time', n_cols=6):
+    
+    f = nc.NetCDFFile(filename)
+    data = f.variables  # dictionary of variable names to nc_var objects
+    dims = f.dimensions # dictionary of dimension names to sizes
+    x = data[y_name]
+    y = data[x_name]
+    t = data[t_name]
+    grid = data[grid_name]
+    
+    assert len(x.shape) == 1
+    assert len(y.shape) == 1
+    assert len(t.shape) == 1
+        
+    grid_dims = grid.dimensions # tuple of dimension names
+    name_to_idx = dict((k, i) for i, k in enumerate(grid_dims))
+    
+    grid_t_idx = name_to_idx[t.dimensions[0]]
+    grid_x_idx = name_to_idx[x.dimensions[0]]
+    grid_y_idx = name_to_idx[y.dimensions[0]]
+        
+    n_frames = t.shape[0]
+    # n_cols = 6
+    n_rows = int(ceil( float(n_frames) / n_cols ))
+    
+    w, h = figaspect(float(n_rows)/n_cols)
+
+    colormap = get_cmap('gist_earth')
+    grey_color = (0.5,)*3
+    frame_color = (0.2,)*3
+    
+    density_maxes = []
+    total_counts = []
+    all_t = []
+        
+    xedge = centers_to_edges(x)
+    x_range = xedge.max() - xedge.min()
+    yedge = centers_to_edges(y)
+    y_range = yedge.max() - yedge.min()
+    dx = (xedge[1]-xedge[0])
+
+    # count_scale_factor = dx # / 1000.0
+    # max_count_baseline = 450 * count_scale_factor #/ 10.0
+    min_count, max_count = 1, grid[:].max() #max_count_baseline*(t[1]-t[0])
+
+    f = figure(figsize=(w,h))
+    p = small_multiples_plot(fig=f, rows=n_rows, columns=n_cols)
+    p.label_edges(True)
+    pad = 0.0 # for time labels in each frame
+            
+    for ax in p.multiples.flat:
+        ax.set_axis_bgcolor('black')
+        ax.spines['top'].set_edgecolor(frame_color)
+        ax.spines['bottom'].set_edgecolor(frame_color)
+        ax.spines['left'].set_edgecolor(frame_color)
+        ax.spines['right'].set_edgecolor(frame_color)
+    #     ax.yaxis.set_major_formatter(kilo_formatter)
+    #     ax.xaxis.set_major_formatter(kilo_formatter)
+    base_date = datetime.strptime(t.units, "seconds since %Y-%m-%d %H:%M:%S")
+    time_delta = timedelta(0,float(t[1]-t[0]),0)
+    start_time = base_date + time_delta
+        
+    indexer = [None,]*len(grid.shape)
+
+    for i in range(n_frames):
+        frame_start = base_date + timedelta(0,float(t[i]),0)
+        indexer[grid_t_idx] = i
+        
+        density = grid[indexer]
+        
+        # density,edges = np.histogramdd((x,y), bins=(xedge,yedge))
+        density_plot  = p.multiples.flat[i].pcolormesh(xedge,yedge,
+                                   np.log10(density.transpose()), 
+                                   vmin=-0.2,
+                                   vmax=np.log10(max_count),
+                                   cmap=colormap)
+        label_string = frame_start.strftime('%H%M:%S')
+        text_label = p.multiples.flat[i].text(xedge[0]-pad+x_range*.015, yedge[0]-pad+y_range*.015, label_string, color=grey_color, size=6)
+        density_plot.set_rasterized(True)
+        density_maxes.append(density.max())
+        total_counts.append(density.sum())
+        all_t.append(frame_start)
+        print label_string, x.shape, density.max(), density.sum()
+
+    color_scale = ColorbarBase(p.colorbar_ax, cmap=density_plot.cmap,
+                                       norm=density_plot.norm,
+                                       orientation='horizontal')
+
+    # color_scale.set_label('count per pixel')
+    color_scale.set_label('log10(count per pixel)')
+
+    view_x = (xedge.min(), xedge.max())
+    view_y = (yedge.min(), yedge.max())
+    
+    print 'making multiples',
+    p.multiples.flat[0].axis(view_x+view_y)
+    filename = 'LMA-%s_%s_%5.2fkm_%5.1fs.pdf' % (grid_name, start_time.strftime('%Y%m%d_%H%M%S'), dx, time_delta.seconds)
+    f.savefig(filename, dpi=150)
+    print ' ... done'
+    f.clf()
+        
+
+        
+if __name__ == '__main__':
+    do_profile=False
+    if do_profile:
+        import hotshot
+        from hotshot import stats
+        prof = hotshot.Profile("multiples_test_profile")
+        prof.runcall(runtest)
+        prof.close()
+        s=stats.load("multiples_test_profile")
+        s.sort_stats("time").print_stats()
+    else:
+        import sys
+        res = runtest(sys.argv[1], sys.argv[2])
\ No newline at end of file
diff --git a/small_multiples.py b/small_multiples.py
new file mode 100644
index 0000000..fbf829c
--- /dev/null
+++ b/small_multiples.py
@@ -0,0 +1,82 @@
+import matplotlib 
+matplotlib.rc('xtick', labelsize=6) 
+matplotlib.rc('ytick', labelsize=6) 
+
+from pylab import figure, show #, subplot, show
+from numpy import arange
+
+class small_multiples_plot(object):
+    
+    def __init__(self, fig=None, *args, **kwargs):
+        if fig is None:
+            fig = figure()
+        self.fig = fig
+        self.fig.subplots_adjust(bottom=0.20, left = 0.1, right=0.9, top=0.9)
+        self.colorbar_ax = fig.add_axes((0.1, 0.1, 0.8, 0.05))
+        self.multiples = small_multiples(self.fig, **kwargs)
+        
+    def label_edges(self, bool_val):
+        m = self.multiples
+        leftside = m[:,0]
+        for ax in leftside:
+            ax.yaxis.tick_left()
+            ax.yaxis.set_visible(bool_val)
+
+        #last row
+        bottomedge = m[-1,:]
+        for ax in bottomedge:
+            ax.xaxis.tick_bottom()
+            ax.xaxis.set_visible(bool_val)
+
+        
+def small_multiples(f, rows=4, columns=5, margin=(0.0,0.0), zoom_together=True):
+    """ Given a figure f, create linked subplots with given number of rows and columns.
+        Returns an object array of axes instances [rows, columns], with top left being [0,0].
+    """
+    # rows = 4    #number in y direction
+    # columns = 5 #number in x direction
+    
+    f.subplots_adjust(wspace=margin[0], hspace=margin[1])
+
+    # should use N.empty((rows,columns),dtype=object)
+    # and attribute name should perhaps be changed
+    multiples = arange(rows*columns, dtype=object)
+    multiples.shape=(rows, columns)
+    
+    # No axis defined to start with
+    commonaxis=None
+
+    for row in range(rows):
+        for column in range(columns):
+            nth_plot = row*columns + column 
+            ax = f.add_subplot(rows, columns, nth_plot + 1, sharex=commonaxis, sharey=commonaxis)
+            if not commonaxis and zoom_together:
+                commonaxis = ax
+
+            # leaves axes frame, but turns off axis labels and ticks
+            ax.xaxis.set_visible(False)
+            ax.yaxis.set_visible(False)
+            multiples[row, column] = ax
+            # ax.plot(range(10), range(10))
+            # ax.text(1,1,'%i, %i, %i' % (row, column, nth_plot))
+            # print row, column, nth_plot
+            
+    return multiples
+
+
+
+if __name__ == '__main__':
+    f=figure()
+    m=small_multiples(f)
+
+    #first column
+    leftside = m[:,0]
+    for ax in leftside:
+        ax.yaxis.set_visible(True)
+    
+    #last row
+    bottomedge = m[-1,:]
+    for ax in bottomedge:
+        ax.xaxis.set_visible(True)
+    
+    show()