@coroutine
def accumulate_points_on_grid(grid, xedge, yedge, out=None, label='',
        grid_frac_weights=False):
    assert xedge.shape[0] == grid.shape[0]+1
    assert yedge.shape[0] == grid.shape[1]+1
    if out == None:
        out = {}
    # grid = None

    # When we do a calculation like average flash area, we need to sum the
    # areas, and sum the flashes, and divide at the end. Otherwise, if we have
    # a frame that spans multiple data files (and therefore chunks of
    # x,y,weights) we will calculate the sum of the averages due to each chunk
    # instead of getting the true average. Therefore, create a set of grids for
    # tracking the accumulation, and update the final grid with the new
    # accumulation each time through the loop.
    count_hist = grid.copy()
    total_hist = grid.copy()

    have_weights = False

    try:
        while True:
            if grid_frac_weights:
                x, y, weights, grid_frac = (yield)
                # There is an issue with small weights being rounded to
                # zero in histogramdd, so multiply by some large value
                # and divide it back out later.
                # https://github.com/numpy/numpy/issues/9465
                # seems to not be necessary for the dynamic range we have
                # grid_frac = grid_frac.astype('f8')
                # grid_frac_scale = 1.0e5
                # grid_frac = grid_frac.astype('f8')*grid_frac_scale
            else:
                x, y, weights = (yield)
                grid_frac=None
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y), bins=(xedge, yedge),
                    weights=grid_frac, normed=False)

                count_hist += count.astype(count_hist.dtype)
                # if grid_frac_weights:
                #     count /= grid_frac_scale
                if weights is not None:
                    have_weights = True
                    # histogramdd sums up weights in each bin for normed=False
                    if grid_frac is not None:
                        weights = weights*grid_frac
                    total, edges = np.histogramdd((x,y), bins=(xedge, yedge),
                        weights=weights, normed=False)
                    total_hist += total
                    del total, edges
                    
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    if have_weights:
                        bad = (count_hist <= 0)
                        avg = np.asarray(total_hist, dtype='float32')/count_hist
                        avg[bad] = 0.0
                        del bad
                    else:
                        avg = count_hist
                    grid[:] = avg[:].astype(grid.dtype)
                    # grid += count.astype(grid.dtype)
                del count, avg
            del x, y, weights, grid_frac
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def accumulate_points_on_grid_sdev(grid, grid2, xedge, yedge, out=None, label='', grid_frac_weights=True):
    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:
            if grid_frac_weights:
                x, y, weights, grid_frac = (yield)
            else:
                x, y, weights = (yield)
                grid_frac=None
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=grid_frac, normed=False)    
                
                if weights is not None:
                    # histogramdd sums up weights in each bin for normed=False
                    total, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=np.asarray(weights)**2., 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 
                    
                    del total, edges, bad
                    
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    grid += count.astype(grid.dtype) 
                    grid = np.sqrt(grid - (grid2)**2.)
                del count
            del x, y, weights
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def accumulate_energy_on_grid(grid, xedge, yedge, out=None, label='', grid_frac_weights=True):
    """
    Like accumulate_points_on_grid, but doesn't normalize by the total count
    """
    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:
            if grid_frac_weights:
                x, y, weights, grid_frac = (yield)
            else:
                x, y, weights = (yield)
                grid_frac=None
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y), bins=(xedge, yedge), weights=grid_frac, normed=False)    
                
                if weights is not None:
                    have_weights = True
                    # histogramdd sums up weights in each bin for normed=False
                    if grid_frac is not None:
                        weights = weights*grid_frac
                    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
                    del total, edges, bad
    
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    grid += count.astype(grid.dtype)
                del count
            del x, y, weights
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def accumulate_points_on_grid_3d(grid, xedge, yedge, zedge, out=None, label=''):
    assert xedge.shape[0] == grid.shape[0]+1
    assert yedge.shape[0] == grid.shape[1]+1
    assert zedge.shape[0] == grid.shape[2]+1
    if out == None:
        out = {}
    # grid = None
    try:
        while True:
            x, y, z, weights = (yield)
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                z = np.atleast_1d(z)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), weights=None, normed=False)    
                
                if weights != None:
                    # histogramdd sums up weights in each bin for normed=False
                    total, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), 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 
                    
                    del total, edges, bad
        
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    grid += count.astype(grid.dtype)
                del count
            del x, y, z, weights
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def accumulate_points_on_grid_sdev_3d(grid, grid2, xedge, yedge, zedge, out=None, label=''):
    assert xedge.shape[0] == grid.shape[0]+1
    assert yedge.shape[0] == grid.shape[1]+1
    assert zedge.shape[0] == grid.shape[2]+1
    if out == None:
        out = {}
    # grid = None
    try:
        while True:
            x, y, z, weights = (yield)
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                z = np.atleast_1d(z)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), weights=None, normed=False)    
                
                if weights != None:
                    # histogramdd sums up weights in each bin for normed=False
                    total, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), weights=np.asarray(weights)**2., 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 
                    
                    del total, edges, bad
                    
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    grid += count.astype(grid.dtype)
                    grid = np.sqrt(grid - (grid2)**2.)
                del count
            del x, y, z, weights
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def accumulate_energy_on_grid_3d(grid, xedge, yedge, zedge, out=None, label=''):
    assert xedge.shape[0] == grid.shape[0]+1
    assert yedge.shape[0] == grid.shape[1]+1
    assert zedge.shape[0] == grid.shape[2]+1
    if out == None:
        out = {}
    # grid = None
    try:
        while True:
            x, y, z, weights = (yield)
            if len(x) > 0:
                x = np.atleast_1d(x)
                y = np.atleast_1d(y)
                z = np.atleast_1d(z)
                log.info(('accumulating ', len(x), 'points for ', label))
                count, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), weights=None, normed=False)    
                
                if weights != None:
                    # histogramdd sums up weights in each bin for normed=False
                    total, edges = np.histogramdd((x,y,z), bins=(xedge, yedge, zedge), weights=np.abs(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 
                    
                    del total, edges, bad
                    
                # 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 is None:
                    grid = count
                    out['out'] = grid
                else:
                    grid += count.astype(grid.dtype)
                del count
            del x, y, z, weights
            gc.collect()
    except GeneratorExit:
        out['out'] = grid @coroutine
def point_density(target, weight_key=None, weight_flashes=True,
        flash_id_key='flash_id', event_grid_area_fraction_key=None):
    """ Sends event x, y, z location directly. If weight_key is provided
    also extract the weights from the flash data with variable name matching
    weight_key. if weight_flashes=False, use the event data instead of the
    flash data.
    """
    while True:
        events, flash, x, y, z = (yield)
        # print 'Doing point density',
        x = np.atleast_1d(x)
        y = np.atleast_1d(y)
        if len(x) > 0:
            if weight_key is not None:
                if weight_flashes:
                    weight_lookup = dict(list(zip(flash[flash_id_key],
                        flash[weight_key])))
                    #puts weights in same order as x, y
                    weights = np.fromiter((weight_lookup[fi] for fi in
                        events['flash_id']), dtype='float64')
                else:
                    weights = events[weight_key]
            else:
                weights = None
            log.debug('with points numbering %s'.format(len(x)))
            target.send((x, y, weights))
        del events, flash ,x,y,z @coroutine
def extent_density(x0, y0, dx, dy, target, flash_id_key='flash_id', 
    weight_key=None, event_grid_area_fraction_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')
        test_flash_id = 53735
        if len(x_i) > 0:
            log.info(('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(list(zip(flash[flash_id_key], flash[weight_key])))
                #puts weights in same order as x[unq_idx], y[unq_idx]
                weights = np.fromiter((weight_lookup[fi] for fi in
                    events[unq_idx]['flash_id']), dtype='float64')
                # del weight_lookup
            else:
                weights = None
            if event_grid_area_fraction_key is not None:
                # Each event with a unique index is replicated above
                # with the representative value for the flash (weights = None
                # implies a weight of +1 for each flash). If there is knowledge
                # of how much of the underlying grid cell each event fills
                # (e.g. from pixel-based event detector), then we can modify
                # the weights by how much of the grid cell is filled.
                # The logic here presumes that any of the events in the grid
                # cell cover as much area as any other, i.e., that the pixels
                # doing the event detection don't move during the time of the
                # flash.
                grid_frac = events[unq_idx][event_grid_area_fraction_key]
            else:
                grid_frac = None
                
            # Diagnostics
            # test_flash_mask = (events['flash_id'] == test_flash_id)
            # test_events = events[test_flash_mask]
            # if (test_flash_mask.sum() > 0) & (weight_key == 'area'):
            #     print("Data for flash {0}".format(test_flash_id))
            #     mesh_xi = test_events['mesh_xi']
            #     mesh_yi = test_events['mesh_yi']
            #     mesh_frac = test_events['mesh_frac']
            #     mesh_t = test_events['time']
            #     for vals in zip(mesh_t, mesh_frac,
            #                     mesh_xi, x_i[test_flash_mask],
            #                     mesh_yi, y_i[test_flash_mask],
            #                     ):
            #         print(vals, weight_lookup[test_flash_id])
            #
            # test_flash_mask = (events[unq_idx]['flash_id'] == test_flash_id)
            # test_events = events[unq_idx][test_flash_mask]
            # if (test_flash_mask.sum() > 0) & (weight_key == 'area'):
            #     print("Unique data for flash {0}".format(test_flash_id))
            #     mesh_xi = test_events['mesh_xi']
            #     mesh_yi = test_events['mesh_yi']
            #     mesh_frac = test_events['mesh_frac']
            #     mesh_t = test_events['time']
            #     for vals in zip(mesh_t, mesh_frac,
            #                     mesh_xi, x_i[unq_idx][test_flash_mask],
            #                     mesh_yi, y_i[unq_idx][test_flash_mask],
            #                     weights[test_flash_mask]):
            #         print(vals, weight_lookup[test_flash_id])
            
            target.send((x[unq_idx], y[unq_idx], weights, grid_frac))
            del weights, grid_frac, unq_idx
        # else:
            # print ''
        del events, flash, x, y, z, x_i, y_i @coroutine
def point_density_3d(target, weight_key=None, weight_flashes=True,
        flash_id_key='flash_id'):
    """ Sends event x, y, z location directly. If weight_key is provided
    also extract the weights from the flash data with variable name matching 
    weight_key. if weight_flashes=False, use the event data instead of the
    flash data.
    """
    while True:
        events, flash, x, y, z = (yield)
        # print 'Doing point density',
        if len(x) > 0:
            if weight_key is not None:
                if weight_flashes:
                    weight_lookup = dict(list(zip(flash[flash_id_key], 
                        flash[weight_key])))
                    #puts weights in same order as x, y
                    weights = np.fromiter((weight_lookup[fi] for fi in
                        events['flash_id']), dtype='float64')
                else:
                    weights = events[weight_key]
            else:
                weights = None
            log.debug('with points numbering %s'.format(len(x)))
            target.send((x, y, z, weights))
        del events, flash ,x,y,z @coroutine
def extent_density_3d(x0, y0, z0, dx, dy, dz, 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')
        z_i = np.floor( (z-z0)/dz ).astype('int32')
        log.debug(len(x_i))
        if len(x_i) > 0:
            log.info(('extent with points numbering', len(x_i), ' with weights', weight_key))
            unq_idx = unique_vectors(x_i, y_i, z_i, events[flash_id_key])
            # if x[unq_idx].shape[0] > 1:
            if weight_key != None:
                weight_lookup = dict(list(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]
                del weight_lookup
            else:
                weights = None
            target.send((x[unq_idx], y[unq_idx], z[unq_idx], weights)) 
            del weights, unq_idx
        # else:
            # print ''
        del events, flash, x, y, z, x_i, y_i, z_i @coroutine
def project(x_coord, y_coord, z_coord, mapProj, geoProj, target, 
            use_flashes=False, transform=True):
    """ Adds projected coordinates to the flash and events stream"""
    while True:
        events, flashes = (yield)
        if use_flashes==True:
            points = flashes
        else:
            points = events
        if transform:
            x,y,z = mapProj.fromECEF(*geoProj.toECEF(
                        points[x_coord], points[y_coord], points[z_coord]))
        else:
            x,y,z = points[x_coord], points[y_coord], points[z_coord]
        target.send((events, flashes, np.atleast_1d(x),
                     np.atleast_1d(y), np.atleast_1d(z)))
        del events, flashes, x,y,z, points @coroutine
def flashes_to_frames(time_edges, targets, time_key='start', do_events=False,
    time_edges_datetime=None, flash_counter=None):
    """ time_edges_datetime is same len as time_edges but with datetime objects
        instead of floats.

        When paired with extract_events_for_flashes, and events=False, the
        flashes are placed in the correct time frame, and any events from that
        flash, including those that cross a time boundary, are included.

        if do_events='event_time_key', then also subset the events. This
        operation is naive, i.e., the events are selected by time with no
        attempt to keep events together with their parent flash. Therefore, it
        is important to ensure that events and flashes are sent together in
        chunks that do not cross time boundaries, which implies pre-aggregating
        and time-tagging the event data so that the events and flashes remain
        together when naively subset. If those conditions are met then this
        option allows one to set up a pipeline without an additional
        extract_events_for_flashes step.
    """
    if time_edges_datetime is None:
        # print "Datetime-style time edges not found, using time edges in seconds for flash count label"
        time_edges_datetime = time_edges
    
    flash_count_messages = []
    
    assert len(time_edges) == (len(time_edges_datetime))
    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:])]
        if do_events != False:
            ev_start_times = events[do_events]
            ev_sort_idx = np.argsort(ev_start_times)
            ev_idx = np.searchsorted(ev_start_times[ev_sort_idx], time_edges)
            ev_slices = [slice(*i) for i in zip(ev_idx[0:-1], ev_idx[1:])]
        else:
            ev_slices = range(len(time_edges))
        for target, s, ev_s, frame_start_time in zip(targets,
                    slices, ev_slices, time_edges_datetime[:-1]):
            these_flashes = flashes[sort_idx][s]
            if do_events != False:
                these_events = events[ev_sort_idx][ev_s]
            else:
                these_events = events
            if flash_counter is not None:
                flash_counter.send((these_flashes, frame_start_time))
            # flash_count_status = "Sending %s flashes to frame starting at %s" % (len(these_flashes), frame_start_time)
            # flash_count_messages += flash_count_status
            # print flash_count_status
            target.send((these_events, these_flashes))
        del events, flashes, start_times, sort_idx, idx, slices
    log.info(flash_count_messages) @coroutine
def flash_count_log(logfile, format_string="%s flashes in frame starting at %s"):
    """ Write flash count for some frame to a file-like object. File open/close should be handled
        by the calling routine."""
        
    # Track flash count for each frame
    frame_times = {}
    try:
        while True:
            # Receive list of flashes, frame start time
            flashes, frame_start_time = (yield)
            n_flashes = len(flashes)
        
            try:
                frame_times[frame_start_time] += n_flashes
            except KeyError:
                # Key doesn't exist, so can't increment flash count 
                frame_times[frame_start_time]  = n_flashes
    except GeneratorExit:
        all_times = list(frame_times.keys())
        all_times.sort()
        for frame_start_time in all_times:
            flash_count_status = format_string % (frame_times[frame_start_time], frame_start_time)
            if hasattr(logfile, 'write'):
                logfile.write(flash_count_status+'\n')
            else:
                logfile.info(flash_count_status) @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
        event_dtype = evs.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))
        del events, evs, fls @coroutine
def broadcast(targets):
    while True:
        stuff = (yield)
        for target in targets:
            target.send(stuff)
        del stuff def read_flashes(h5LMAfiles, target, base_date=None, min_points=10):
    """ This routine is the data pipeline source, responsible for pushing out
        events and flashes. Using a different flash data source format is a matter of
        replacing this routine to read in the necessary event and flash data."""

    # for filename in h5LMAfiles:
        # h5 = LMAh5File(filename, min_points=min_points)
        # for events, flashes in h5.gen_events_flashes():
        #
        #     # get the start date of the file, which is the seconds-of-day reference
        #     h5_start =  h5.base_date
        #     if base_date==None:
        #         base_date = h5_start
        #
        #     # adjust time in seconds of day to be in correct reference
        #     extra_dt = to_seconds(h5_start - base_date)
        #     events['time'] += extra_dt
        #     flashes['start'] += extra_dt
        #
        #     n_flashes = len(flashes)
        #     print('{0} -- {1} flashes > {2} pts; dt+={3}  '.format(filename, n_flashes, min_points, extra_dt))

    h5s = LMAh5Collection(h5LMAfiles, base_date=base_date, min_points=min_points)
    for push_out in h5s: # get a stream of events, flashes
        target.send(push_out)
        del push_out def to_seconds(dt):
    'convert datetime.timedelta object to float seconds'
    return dt.days * 86400 + dt.seconds + dt.microseconds/1.0e6 class MapProjection(CoordinateSystem):
    """Map projection coordinate system. Wraps proj4, and uses its projecion names. Defaults to
        equidistant cylindrical projection
    """

    def __init__(self, projection='eqc', ctrLat=None, ctrLon=None, ellipse='WGS84', datum='WGS84', **kwargs):
        self.ERSxyz = proj4.Proj(proj='geocent', ellps=ellipse, datum=datum)
        self.projection = proj4.Proj(proj=projection, ellps=ellipse, datum=datum, **kwargs)
        self.ctrLat=ctrLat
        self.ctrLon=ctrLon
        self.ctrAlt=0.0
        self.geoCS = GeographicSystem()
        self.cx, self.cy, self.cz = 0, 0, 0
        self.cx, self.cy, self.cz = self.ctrPosition()

    def ctrPosition(self):
        if (self.ctrLat != None) & (self.ctrLon != None):
            ex, ey, ez = self.geoCS.toECEF(self.ctrLon, self.ctrLat, self.ctrAlt)
            cx, cy, cz = self.fromECEF(ex, ey, ez)
        else:
            cx, cy, cz = 0, 0, 0
        return cx, cy, cz

    def toECEF(self, x, y, z):
        x += self.cx
        y += self.cy
        z += self.cz
        projectedData = array(proj4.transform(self.projection, self.ERSxyz, x, y, z ))
        if len(projectedData.shape) == 1:
            px, py, pz = projectedData[0], projectedData[1], projectedData[2]
        else:
            px, py, pz = projectedData[0,:], projectedData[1,:], projectedData[2,:]
        return px, py, pz

    def fromECEF(self, x, y, z):
        projectedData = array(proj4.transform(self.ERSxyz, self.projection, x, y, z ))
        if len(projectedData.shape) == 1:
            px, py, pz = projectedData[0], projectedData[1], projectedData[2]
        else:
            px, py, pz = projectedData[0,:], projectedData[1,:], projectedData[2,:]
        return px-self.cx, py-self.cy, pz-self.cz class GeographicSystem(CoordinateSystem):
    """
    Coordinate system defined on the surface of the earth using latitude,
    longitude, and altitude, referenced by default to the WGS84 ellipse.

    Alternately, specify the ellipse shape using an ellipse known
    to pyproj, or [NOT IMPLEMENTED] specify r_equator and r_pole directly.
    """
    def __init__(self, ellipse='WGS84', datum='WGS84',
                 r_equator=None, r_pole=None):
        if (r_equator is not None) | (r_pole is not None):
            if r_pole is None:
                r_pole=r_equator
            self.ERSlla = proj4.Proj(proj='latlong', a=r_equator, b=r_pole)
        else:
            # lat lon alt in some earth reference system
            self.ERSlla = proj4.Proj(proj='latlong', ellps=ellipse, datum=datum)
        self.ERSxyz = proj4.Proj(proj='geocent', ellps=ellipse, datum=datum)
    def toECEF(self, lon, lat, alt):
        lat = atleast_1d(lat) # proj doesn't like scalars
        lon = atleast_1d(lon)
        alt = atleast_1d(alt)
        if (lat.shape[0] == 0): return lon, lat, alt # proj doesn't like empties
        projectedData = array(proj4.transform(self.ERSlla, self.ERSxyz, lon, lat, alt ))
        if len(projectedData.shape) == 1:
            return projectedData[0], projectedData[1], projectedData[2]
        else:
            return projectedData[0,:], projectedData[1,:], projectedData[2,:]

    def fromECEF(self, x, y, z):
        x = atleast_1d(x) # proj doesn't like scalars
        y = atleast_1d(y)
        z = atleast_1d(z)
        if (x.shape[0] == 0): return x, y, z # proj doesn't like empties
        projectedData = array(proj4.transform(self.ERSxyz, self.ERSlla, x, y, z ))
        if len(projectedData.shape) == 1:
            return projectedData[0], projectedData[1], projectedData[2]
        else:
            return projectedData[0,:], projectedData[1,:], projectedData[2,:] def write_cf_netcdf_latlon(outfile, t_start, t, xloc, yloc, 
    lon_for_x=None, lat_for_y=None, ctr_lat=None, ctr_lon=None, 
    grid=None, grid_var_name='lightning',
    grid_description='some kind of lightning data', 
    format='i', missing_value=-9999, is_latlon=True, is_fixedgrid=False,
    **kwargs):
def write_cf_netcdf_3d(outfile, t_start, t, xloc, yloc, zloc, lon_for_x, lat_for_y, alt_for_z, ctr_lat, ctr_lon, ctr_alt, grid, grid_var_name, grid_description, format='i', **kwargs):
def write_cf_netcdf_3d_latlon(outfile, t_start, t, xloc, yloc, zloc, lon_for_x, lat_for_y, alt_for_z, ctr_lat, ctr_lon, ctr_alt, grid, grid_var_name, grid_description, format='i', **kwargs):
import glob
import os, sys
import logging
import numpy as np
from datetime import datetime, timedelta
from .density_to_files import (accumulate_points_on_grid,
    accumulate_points_on_grid_sdev, accumulate_energy_on_grid,
    accumulate_points_on_grid_3d, accumulate_points_on_grid_sdev_3d,
    accumulate_energy_on_grid_3d,
    point_density, extent_density, point_density_3d, extent_density_3d, project,
    flashes_to_frames, flash_count_log, extract_events_for_flashes)
from lmatools.stream.subset import broadcast
from lmatools.io.LMA_h5_file import read_flashes, to_seconds
from lmatools.coordinateSystems import MapProjection, GeographicSystem
from .cf_netcdf import (write_cf_netcdf, write_cf_netcdf_3d,
    write_cf_netcdf_latlon, write_cf_netcdf_3d_latlon,
    write_cf_netcdf_noproj, write_cf_netcdf_fixedgrid)
        zedge=np.arange(z_bnd[0], z_bnd[1]+dz, dz)

        if self.proj_name == 'latlong':
            dx_units = '{0:6.4f}deg'.format(float(dx))
            mapProj = GeographicSystem()
        elif self.proj_name == 'pixel_grid':
            dx_units = 'pixel'
            mapProj = pixel_coords
        elif self.proj_name == 'geos':
            dx_units = '{0:03d}urad'.format(int(dx*1e6))
            mapProj = pixel_coords
        else:
            dx_units = '{0:5.1f}m'.format(float(dx))
            mapProj = MapProjection(projection=self.proj_name, ctrLat=ctr_lat, ctrLon=ctr_lon, lat_ts=ctr_lat,
                                lon_0=ctr_lon, lat_0=ctr_lat, lat_1=ctr_lat, ellipse=proj_ellipse, datum=proj_datum)
        geoProj = GeographicSystem()

        self.t_ref = t_ref
        self.t_edges = t_edges
        self.t_edges_seconds = t_edges_seconds
        self.duration = duration
        self.xedge = xedge
        self.yedge = yedge
        self.zedge = zedge
        self.mapProj = mapProj
        self.geoProj = geoProj
        self.dx_units = dx_units

        self.pipeline_setup()
        self.output_setup()

    def pipeline_setup(self):
        """
        Create the grids for each variable of interest.

        Set up the data processing pipeline up to the point at which data are
        written to grids in memory.
        """
        event_grid_area_fraction_key=self.event_grid_area_fraction_key
        energy_grids=self.energy_grids
        n_frames = self.n_frames
        xedge, yedge, zedge = self.xedge, self.yedge, self.zedge
        dx, dy, dz = self.dx, self.dy, self.dz
        x0 = xedge[0]
        y0 = yedge[0]
        z0 = zedge[0]
        mapProj = self.mapProj
        geoProj = self.geoProj

        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')
        if self.event_grid_area_fraction_key is not None:
            extent_density_grid = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='float32')
        else:
            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')

        specific_energy_grid = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='float32')
        total_energy_grid = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='float32')
        flashsize_std_grid  = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, n_frames), dtype='float32')

        if self.do_3d == True:
            event_density_grid_3d  = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='int32')
            init_density_grid_3d   = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='int32')
            extent_density_grid_3d = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='int32')
            footprint_grid_3d      = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='float32')

            specific_energy_grid_3d = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='float32')
            total_energy_grid_3d    = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='float32')
            flashsize_std_grid_3d   = np.zeros((xedge.shape[0]-1, yedge.shape[0]-1, zedge.shape[0]-1, n_frames), dtype='float32')

        self.outgrids = (extent_density_grid,
                    init_density_grid,
                    event_density_grid,
                    footprint_grid,
                    specific_energy_grid,
                    flashsize_std_grid,
                    total_energy_grid,
                    )

        if self.do_3d == True:
            self.outgrids_3d = (extent_density_grid_3d,
                    init_density_grid_3d,
                    event_density_grid_3d,
                    footprint_grid_3d,
                    specific_energy_grid_3d,
                    flashsize_std_grid_3d,
                    total_energy_grid_3d
                    )
        else:
            self.outgrids_3d = None


        all_frames = []
        for i in range(n_frames):
            extent_out = {'name':'extent'}
            init_out   = {'name':'init'}
            event_out  = {'name':'event'}
            std_out    = {'name':'std'}

            extent_out_3d = {'name':'extent_3d'}
            init_out_3d   = {'name':'init_3d'}
            event_out_3d  = {'name':'event_3d'}
            std_out_3d    = {'name':'std_3d'}

            accum_event_density  = accumulate_points_on_grid(event_density_grid[:,:,i], xedge, yedge,  out=event_out, label='event')
            accum_init_density   = accumulate_points_on_grid(init_density_grid[:,:,i], xedge, yedge,   out=init_out,  label='init')
            accum_extent_density = accumulate_points_on_grid(extent_density_grid[:,:,i], xedge, yedge, out=extent_out,label='extent',  grid_frac_weights=True)
            accum_footprint      = accumulate_points_on_grid(footprint_grid[:,:,i], xedge, yedge, label='footprint',  grid_frac_weights=True)

            accum_specific_energy = accumulate_energy_on_grid(specific_energy_grid[:,:,i], xedge, yedge, out=extent_out, label='specific_energy',  grid_frac_weights=True)
            accum_flashstd       = accumulate_points_on_grid_sdev(flashsize_std_grid[:,:,i], footprint_grid[:,:,i], xedge, yedge, out=extent_out, label='flashsize_std',  grid_frac_weights=True)
            accum_total_energy   = accumulate_energy_on_grid(total_energy_grid[:,:,i], xedge, yedge, out=extent_out, label='total_energy',  grid_frac_weights=True)

            if self.do_3d == True:
                accum_event_density_3d  = accumulate_points_on_grid_3d(event_density_grid_3d[:,:,:,i], xedge, yedge, zedge,  out=event_out_3d, label='event_3d')
                accum_init_density_3d   = accumulate_points_on_grid_3d(init_density_grid_3d[:,:,:,i], xedge, yedge, zedge,   out=init_out_3d,  label='init_3d')
                accum_extent_density_3d = accumulate_points_on_grid_3d(extent_density_grid_3d[:,:,:,i], xedge, yedge, zedge, out=extent_out_3d,label='extent_3d')
                accum_footprint_3d      = accumulate_points_on_grid_3d(footprint_grid_3d[:,:,:,i], xedge, yedge, zedge, label='footprint_3d')


                accum_specific_energy_3d = accumulate_energy_on_grid_3d(specific_energy_grid_3d[:,:,:,i], xedge, yedge, zedge, out=extent_out_3d,label='specific_energy_3d')