class LMAgridFileCollection(object):
    def __init__(self, filenames, grid_name, 
                 x_name='x', y_name='y', t_name='time'):
        """ Load NetCDF files given by filenames and a known grid type.
            The NetCDF files must have been created by lmatools,
            though you might have some luck if the grids are cf-compliant
            more generally.
        
            NCs = LMAgridFileCollection(filenames, 'lma_source')
            
            The data are retrievable by iteration over the 
            collection of files:
            
            >>> for t, xedge, yedge, data in NCs:
            >>>     print(t)

            Or, if you know a time accurately, you can do:
            >>> from datetime import datetime
            >>> t = datetime(2013,6,6,3,0,0)
            >>> xedge, yedge, data = NCs.data_for_time(t)

            The data have been oriented such that a call to matplotlib's
            pcolormesh(xedge,yedge,data)
            will do the expected thing.
        
        """

        self.x_name = x_name
        self.y_name = y_name
        self.t_name = t_name
        self.grid_name = grid_name
        self._filenames = filenames
#         self._ncfiles = [NetCDFFile(f) for f in filenames]
        self._time_lookup = {} # populated by self._frame_times_for_file
        self.times = [t for t in self._all_times()]
        self.times.sort()
        
    def data_for_time(self, t0, return_nc=False):
        """ Read data from the file corresponding to datetime t.
            Returns xedge, yedge, and density for the file 
        
            If return_nc is True, also return the NetCDFFile corresponding
            to the data as a fourth return value
        """
        fname, i = self._time_lookup[t0]  #i is the frame id for this time in NetCDFFile f
        
        f = NetCDFFile(fname)
        
        data = f.variables  # dictionary of variable names to nc_var objects
        dims = f.dimensions # dictionary of dimension names to sizes
        x = data[self.x_name]
        y = data[self.y_name]
        t = data[self.t_name]
        grid = data[self.grid_name]
        indexer = [slice(None),]*len(grid.shape)
        
        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]]
                
        xedge = centers_to_edges(x)
        yedge = centers_to_edges(y)
        indexer[grid_t_idx] = i
        density = grid[indexer].transpose()
        out = xedge, yedge, density
        if return_nc:
            out += (f,)
        else:
            f.close()
        return out
        
    def _all_times(self):
        for f in self._filenames:
            for t in self._frame_times_for_file(f):
                yield t    
        
    def _frame_times_for_file(self, fname):
        """ Called once by init to set up frame lookup tables and yield 
            the frame start times. _frame_lookup goes from 
            datetime->(nc file, frame index)"""
        f = NetCDFFile(fname)

        data = f.variables  # dictionary of variable names to nc_var objects
        dims = f.dimensions # dictionary of dimension names to sizes
        t = data[self.t_name]
        
        try:
            base_date = datetime.strptime(t.units, "seconds since %Y-%m-%d %H:%M:%S")
        except ValueError:
            base_date = datetime.strptime(t.units, "seconds since %Y-%m-%d")
        for i in range(np.atleast_1d(t).shape[0]):
            frame_start = base_date + timedelta(0,float(t[i]),0)
            self._time_lookup[frame_start] = (fname, i)
            yield frame_start
        f.close()
        
    def get_projection(self):
        """ Returns GeographicSystem and MapProjection instances from
            lmatools.coordinateSystems corresponding
            to the coordinate system specified by the metadata of the
            first NetCDF file in self._filenames.
        """
        from lmatools.coordinateSystems import GeographicSystem, MapProjection
        geosys = GeographicSystem()
        
        f = NetCDFFile(self._filenames[0])
        
        # Surely someone has written an automated library to parse coordinate 
        # reference data from CF-compliant files.
        if 'Lambert_Azimuthal_Equal_Area' in list(f.variables.keys()):
            nc_proj = f.variables['Lambert_Azimuthal_Equal_Area']
            proj_name = 'laea'
            ctrlon, ctrlat  = (nc_proj.longitude_of_projection_origin,
                                      nc_proj.latitude_of_projection_origin,)
            try:
                ctralt = nc_proj.altitude_of_projection_origin
            except AttributeError:
                print("No altitude attribute in NetCDF projection data, setting to 0.0")
                ctralt = 0.0
            mapproj = MapProjection(proj_name, ctrLat = ctrlat, ctrLon=ctrlon, 
                                     lat_0=ctrlat, lon_0=ctrlon)
            # print geosys.fromECEF(*mapproj.toECEF((0,0), (0,0), (0,0)))
            return geosys, mapproj
        else:
            print("projection not found, assuming lat, lon grid")    
            return geosys, geosys
        f.close()

        
                    
    def __iter__(self):
        for t in self.times:
            xedge, yedge, density = self.data_for_time(t)
            yield t, xedge, yedge, density def plot_energy_from_area_histogram(histo, bin_edges, bin_unit='km', save=False, fig=None, color_cycle_length=1, color_map='gist_earth', duration=600.0):
    """ Histogram for flash width vs. count """
    
    fig, spectrum_ax, fivethirds_line_artist, spectrum_artist = energy_plot_setup()
    spectrum_ax.set_title(save.split('/')[-1].split('.')[0])
    
    flash_1d_extent, specific_energy = calculate_energy_from_area_histogram(histo, bin_edges, duration)
    spectrum_artist.set_data(flash_1d_extent, specific_energy)
    
    if save==False:
        plt.show()
    else:
        # ax.set_title(save)
        fig.savefig(save)
        fig.clf() def get_energy_spectrum_bins(min_pwr=-2, max_pwr = 4, delta_pwr=0.1):
    #defaults are 1e-2 to 1e4
    powers = np.arange(min_pwr, max_pwr+delta_pwr, delta_pwr)
    footprint_bin_edges = 10**powers
    return footprint_bin_edges def bin_center(bins):
    return (bins[:-1] + bins[1:]) / 2.0 def plot_energies(footprint_bin_edges,time_array,scalar_map,flashes_series,
                    flashes_in_poly_edges,spectrum_save_file_base_en,which_energy,title):
    '''
    Plots estimated total and specific energy fields in .h5 files. 
    Arguments:
            -) flash area bin edges: footprint_bin_edges
            -) times in polygon:     time_array
            -) color by time cmap:   scalar_map
            -) radar specific times: time_mask   
            -) flash table series:   flashes_series
            -) flashes in lasso:     flashes_in_poly_edges
            -) save path:            spectrum_save_file_base_en
            -) energy field:         which_energy
            -) title of plot:        title             
    '''
    import matplotlib.pyplot as plt
    s_m = scalar_map
    figure_energy = plt.figure(figsize=(14,9)) 
    ax_energy = figure_energy.add_subplot(111)
    # cmap = plt.cm.Reds_r
    
    tmin = min(flashes_in_poly_edges)
    tmax = max(flashes_in_poly_edges)
    
    flash_1d_extent = bin_center(np.sqrt(footprint_bin_edges))
    for f, (flashes, t0, t1) in enumerate(zip(flashes_series, flashes_in_poly_edges[:-1], flashes_in_poly_edges[1:])):
        if flashes.shape[0] > 1:
            # If there is no data in this time window, skip it
            histo_cd, edges_cd = np.histogram(np.sqrt(flashes['area']), 
                                              bins=np.sqrt(footprint_bin_edges), 
                                              weights=np.abs(flashes[which_energy]))        
            estimated, = ax_energy.loglog(flash_1d_extent[:],
                                np.abs(np.asarray(histo_cd))/np.sqrt(flash_1d_extent),
                                color=s_m.to_rgba(np.asarray(time_array)[f]),
                                alpha=0.7); 
                                                   
    if which_energy == 'Energy' or which_energy == 'total_energy':
        ax_energy.set_ylim(1e7,1e13)
        # plt.xlim(plt.xlim()[::-1])
        ax_energy.set_xlim(1e2,1e-1)
        wavenumber = (2.*np.pi)/flash_1d_extent
        inertialsubrange = 10**6 * (wavenumber*0.0002)**(-5.0/3.0)
        
    else:
        wavenumber = (2.*np.pi)/flash_1d_extent
        inertialsubrange = 10**1 * (wavenumber)**(5.0/3.0)
        ax_energy.set_xlim(1e2,1e-1)
        ax_energy.set_xlim(plt.xlim()[::-1])

    spectrum_save_file_en = spectrum_save_file_base_en.format(tmin.strftime('%y%m%d%H%M%S'),
                                                              tmax.strftime('%y%m%d%H%M%S'))

    cbar = plt.colorbar(s_m)
    ax_energy.loglog(flash_1d_extent, inertialsubrange,'k-',alpha=0.5);
    ax_energy.set_title('Estimated {0} Spectra for Cell Selection'.format(title),fontsize=15)
    ax_energy.set_xlabel(r'Flash width ($\sqrt{A_h}$, $km$)',fontsize=15)
    ax_energy.set_ylabel(r'Energy ($J$)',fontsize=15)
    cbar.ax.tick_params(labelsize=15)
    ax_energy.tick_params(labelsize=15)
    plt.savefig(spectrum_save_file_en)
    plt.close() def flash_size_stats(flashes):
    """ Using the area data in the flashes table, calculate the flash size moments
        return moments, energy, and energy_per_flash, where energy is the dimensional
        energy defined by Bruning and Thomas (2015, JGR).
    
        Returned is a single row array with column names as follows:
        result_dtype = [
            ('number','f4'), 
            ('mean', 'f4')
            ('variance', 'f4')
            ('skewness', 'f4')
            ('kurtosis', 'f4')
            ('energy', 'f4')
            ('energy_per_flash', 'f4')
           ]
    """
    area = flashes['area']
    length_scale = np.sqrt(area) # length scale per Bruning and Thomas (2015, JGR)
    
    result_dtype = [
        ('number','f4'), 
        ('mean', 'f4'),
        ('variance', 'f4'),
        ('skewness', 'f4'),
        ('kurtosis', 'f4'),
        ('energy', 'f4'),
        ('energy_per_flash', 'f4'),
        ('specific_energy', 'f4'),
        ('total_energy', 'f4'),
       ]
    result = np.zeros((1,), dtype=result_dtype)
     
    # Moment calculation reduces all flash lengths into a few representative moments.
    central, standard = central_moments_from_raw(raw_moments(length_scale, n_moments=5))
    number, mean, variance, skewness, kurtosis = standard
    result['number'] = number
    result['mean'] = mean
    result['variance'] = variance
    result['skewness'] = skewness
    result['kurtosis'] = kurtosis
    
    energy = number * (mean*mean + variance) # Bruning and Thomas (2015, JGR)
    result['energy'] = energy
    result['energy_per_flash'] = energy/number
    
    #Moments for energy:
    try:
        result['specific_energy'] = flashes['specific_energy'].sum()
        result['total_energy']  = flashes['total_energy'].sum()
    except ValueError:
        "Print skipping energy quantities"
    return result def plot_flash_stat_time_series(basedate, t_edges, stats, major_tick_every=1800):
    t_start, t_end = t_edges[:-1], t_edges[1:]
    fig = plt.figure(figsize=(11,8.5))
    starts = np.fromiter( ((s - basedate).total_seconds() for s in t_start), dtype=float )
    ends = np.fromiter( ((e - basedate).total_seconds() for e in t_end), dtype=float )
    t = (starts+ends) / 2.0
    window_size = (ends-starts)/60.0 # in minutes

    ax_energy = fig.add_subplot(2,2,2)
    ax_energy.plot(t, stats['energy']/window_size, label='Total energy')
    ax_energy.legend()

    ax_count = fig.add_subplot(2,2,1, sharex=ax_energy)
    # ax_count.plot(t, stats['energy_per_flash']/window_size, label='$E_T$ flash$^{-1}$ min$^{-1}$ ')
    ax_count.plot(t, stats['number']/window_size, label='Flash rate (min$^{-1}$)')
    ax_count.legend()#location='upper left')

    ax_mean = fig.add_subplot(2,2,3, sharex=ax_energy)
    mean, std = stats['mean'], stats['variance'], 
    skew, kurt = stats['skewness'], stats['kurtosis']
    sigma = np.sqrt(std)
    ax_mean.plot(t, mean, label = 'Mean flash size (km)')
    ax_mean.plot(t, sigma, label = 'Standard deviation (km)')
    ax_mean.set_ylim(0, 20)
    ax_mean.legend()
#         ax_mean.fill_between(t, mean+sigma, mean-sigma, facecolor='blue', alpha=0.5)

    ax_higher = fig.add_subplot(2,2,4, sharex=ax_energy)
    ax_higher.plot(t, skew, label='Skewness')
    ax_higher.plot(t, kurt, label='Kurtosis')
    ax_higher.set_ylim(-2,10)
    ax_higher.legend()

    for ax in fig.get_axes():
        ax.xaxis.set_major_formatter(SecDayFormatter(basedate, ax.xaxis))  
        ax.set_xlabel('Time (UTC)')
        ax.xaxis.set_major_locator(MultipleLocator(major_tick_every))
        ax.xaxis.set_minor_locator(MultipleLocator(major_tick_every/2.0))
    return fig def plot_energy_stats(size_stats, basedate, t_edges, outdir):
    t_start, t_end = t_edges[:-1], t_edges[1:]
    starts = np.fromiter( ((s - basedate).total_seconds() for s in t_start), dtype=float )
    ends = np.fromiter( ((e - basedate).total_seconds() for e in t_end), dtype=float )
    t = (starts+ends) / 2.0
    
    specific_energy = size_stats
    
    figure = plt.figure(figsize=(15,10))
    ax     = figure.add_subplot(111)
    ax.plot(t,specific_energy,'k-',label='Specific Energy',alpha=0.6)
    plt.legend()
    # ax.set_xlabel('Time UTC')
    ax.set_ylabel('Specific Energy (J/kg)')
    
    for axs in figure.get_axes():
        axs.xaxis.set_major_formatter(SecDayFormatter(basedate, axs.xaxis))  
        axs.set_xlabel('Time (UTC)')
        axs.xaxis.set_major_locator(MultipleLocator(1800))
        axs.xaxis.set_minor_locator(MultipleLocator(1800/2))
    
    return figure def plot_tot_energy_stats(size_stats, basedate, t_edges, outdir):
    t_start, t_end = t_edges[:-1], t_edges[1:]
    starts = np.fromiter( ((s - basedate).total_seconds() for s in t_start), dtype=float )
    ends = np.fromiter( ((e - basedate).total_seconds() for e in t_end), dtype=float )
    t = (starts+ends) / 2.0
    
    specific_energy = np.abs(size_stats)
    
    figure = plt.figure(figsize=(15,10))
    ax     = figure.add_subplot(111)
    ax.plot(t,specific_energy,'k-',label='Total Energy',alpha=0.6)
    plt.legend()
    # ax.set_xlabel('Time UTC')
    ax.set_ylabel('Total Energy (J)')
    
    for axs in figure.get_axes():
        axs.xaxis.set_major_formatter(SecDayFormatter(basedate, axs.xaxis))  
        axs.set_xlabel('Time (UTC)')
        axs.xaxis.set_major_locator(MultipleLocator(1800))
        axs.xaxis.set_minor_locator(MultipleLocator(1800/2))
    
    return figure class FractalLengthProfileCalculator(object):
    def __init__(self, D, b_s, alt_bins):
        self.alt_bins = alt_bins
        self.D = D
        self.b_s = b_s
        
    def process_events_flashes(self, events, flashes):
        """ Given a table of flashes and their corresponding events table,
            calculate the Bruning and Thomas (2015, JGR) fractal length
            for each flash and produce a vertical profile of flash lengths
            for all flashes.
        
            Returns:
            per_flash_data: List of dictionaries containing the flash length 
                and triangulation geometry data, including the flash volume.
            IC_totals: array of summary data for each profile; see dtype below.
            CG_totals: as above but for CGs only.
            These data are not normalized by time.
        
            totals_dtype = [
                ('fractal_length_2D_hull', 'f8'),
                ('fractal_length_3D_hull', 'f8'),
                ('source_density_profile', 'f8', alt_bin_shape), 
                ('length_density_profile', 'f8', alt_bin_shape),
                ]   
            
        """
        results = [result for result in gen_fractal_length_for_flashes(
                       events, flashes, self.D, self.b_s, self.alt_bins)
                  ]
                       
        totals_CG = self._sum_profiles(results, is_CG_flag=True)
        totals_IC = self._sum_profiles(results, is_CG_flag=False)

        # IC_totals = self._process_totals(totals_IC, duration)
        # CG_totals = self._process_totals(totals_CG, duration)
        
        return results, totals_IC, totals_CG

    def normalize_profiles(self, totals_seq, durations):
        """ Given a sequence of results from repeated calls to self.process_events_flashes,
            return profiles normalized by durations (an array of durations for each interval
            corresponding to the items in total_seq). also return a source density 
            
            returns:    
            lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin
            total 2D and 3D lengths, source and 3D length profiles, and the maximum 
            length per unit time per km altitude across all bins.
        """
        totals = np.asarray(list(totals_seq))
        L_profile_rate = np.squeeze(totals['length_density_profile'].T/durations)
        max_L_bin = L_profile_rate.max()
        src_profile_rate = totals['source_density_profile'].T/durations
        scaled_sources = np.squeeze(max_L_bin * src_profile_rate / src_profile_rate.max())
    
        lengths_2D = np.squeeze(totals['fractal_length_2D_hull'])/durations
        lengths_3D = np.squeeze(totals['fractal_length_3D_hull'])/durations
        return lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin

                           
    def _sum_profiles(self, results, is_CG_flag=False):
        
        alt_bins = self.alt_bins
        # one fewer density value (bin interval) than bin edges.
        alt_bin_shape = (alt_bins.shape[0]-1,)
        res_dtype = [ # this the dtype of the return value, i.e., the sums
            ('fractal_length_2D_hull', 'f8'),
            ('fractal_length_3D_hull', 'f8'),
            ('source_density_profile', 'f8', alt_bin_shape), 
            ('length_density_profile', 'f8', alt_bin_shape), #, (alt_bins.shape[0]-1,) ,)
            ]   
        
        if len(results) == 0:
            return np.zeros((1,), dtype=res_dtype)
        
        # Check to see if there are CG flags in the flash data table
        # just look at the first flash and assume the rest have it if one does.
        if 'CG' in results[0]['flash'].dtype.names:
            results_iter = (ri for ri in results if
                                    (ri['flash']['CG'] == is_CG_flag) )
        else:
            results_iter = results
            
        # def get_res_iter():
        res_iter = ( (
             r['2D']['length'], 
             r['3D']['length'],
             tuple(r['3D']['vertical_profile']['src_density']),
             tuple(r['3D']['vertical_profile']['length_density'])
            ) for r in results_iter
            )

        # http://stackoverflow.com/questions/19201868/how-to-set-dtype-for-nested-numpy-ndarray
        each_result = np.fromiter(res_iter, dtype=res_dtype)
        total = np.empty((1,), dtype=res_dtype)
        for colname in total.dtype.names:
            total[colname] = each_result[colname].sum(axis=0)
        return total
        
    def make_time_series_plot(self, basedate, t_edges, 
        lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin,
        label_t_every=3600., figsize=(7.5,10)):      
        """ Arguments: 
            t_edges: N+1 bin boundaries corresponding to the N time series values in the
                    other arguments. Units: seconds since start of day.
            lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin: 
                    The values returned by normalize_profiles
            label_t_every: interval in seconds at which to label the time axis
            figsize: (width, height) of figure in inches (passed to matplotlib.figure)
            
        """
        import matplotlib.pyplot as plt
        cmap = 'cubehelix_r'        
    
        fig = plt.figure(figsize=figsize)
        ax_L = fig.add_subplot(311)
        ax_prof_L = fig.add_subplot(312)
        ax_prof_src = fig.add_subplot(313)
        
        min_L_bin = 0*max_L_bin
        
        starts, ends = t_edges[:-1], t_edges[1:]
        t_centers = (starts+ends)/2.0
    
        src_pm = ax_prof_src.pcolormesh(t_edges, self.alt_bins, scaled_sources, 
            cmap=cmap, vmin=min_L_bin, vmax=max_L_bin)
        src_pm.set_rasterized(True)
        cbar_src = fig.colorbar(src_pm, ax=ax_prof_src, orientation='horizontal')
        cbar_src.set_label('Source count per height interval per time\n(scaled to max length, km/km/min)')
    
        L_pm = ax_prof_L.pcolormesh(t_edges, self.alt_bins, L_profile_rate, 
            cmap=cmap, vmin=min_L_bin, vmax=max_L_bin)
        L_pm.set_rasterized(True)
        cbar_L = fig.colorbar(L_pm, ax=ax_prof_L, orientation='horizontal')
        cbar_L.set_label('Length per height interval per time\n(km/km/min)')

        for ax in (ax_prof_src, ax_prof_L):
            ax.set_ylabel('Altitude (km)')
                    
        ax_L.plot(t_centers, lengths_2D, 
            label='2D Hull Area')
        ax_L.plot(t_centers, lengths_3D,
            label='3D Hull Volume')
        ax_L.legend()
        ax_L.set_ylabel('Fractal Length (km/min)')
    
        for ax in (ax_L, ax_prof_L, ax_prof_src):
            ax.xaxis.set_major_formatter(SecDayFormatter(basedate, ax.xaxis))  
            ax.set_xlabel('Time (UTC)')
            ax.xaxis.set_major_locator(MultipleLocator(label_t_every))
        return fig
        
    def write_profile_data(self, basedate, t_edge, outfilebase,
            lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin,
            partition_kind='total'):
        """ Arguments:
            outfilebase: file base name, including directory.
            basedate: date against which t_edge is referenced.
            t_edge: N+1 bin boundaries corresponding to the N time series values in the
                    other arguments. Units: seconds since start of day given by basedate.
            lengths_2D, lengths_3D, scaled_sources, L_profile_rate, max_L_bin: 
                    The values returned by normalize_profiles
            partition_kind: one of 'total', 'IC', 'CG', or another unique tag
                    classifying this kind of profile. Appended to filename just before.

            The filename is calculated as [outfile_base]_[parition_kind].txt

        """
        starts = t_edge[:-1]
        ends = t_edge[1:]
        header = ""
        header += "# LMA channel length distribution\n"
        header += "# Base date = " + basedate.isoformat() +"\n"
        header += "# Altitude_bins = " + str(self.alt_bins.tolist()) +"\n"
        header += "# starts, ends, lengths_2D, lengths_3D, [scaled_sources x Nbins], [L_profile_rate x Nbins]\n"
        text_dump = open(outfilebase+'_{0}.txt'.format(partition_kind), 'w')
        text_dump.write(header)
        # shape of scaled_sources and L_profile_rate are (N_alt_bins, N_times), and loop is over the first dimension, so take transpose so loop is over N_times
        for s0, e1, l2, l3, Sprof, Lprof in zip(starts, ends, lengths_2D, lengths_3D, scaled_sources.T, L_profile_rate.T):
            text_dump.write('{0}, {1}, {2}, {3}, {4}, {5}\n'.format(
                                                    s0, e1, l2, l3,
                                                    str(Sprof.tolist()), str(Lprof.tolist())
                                                    )) def read_poly_log_file(filename, lat_lon_filter=None):
    """ Import polygons from logger tool.
        
        Returns (polys, t_edges)
        polys: length N sequence of (lon,lat) tuples giving the polygon vertices 
        t_edges: length N+1 sequence of datetime instances giving the 
            start and end times to which each of the N polygons applies. Each
            polygon is valid beginning at the frame_time given in the log file
            for each polygon, and ends when next polygon begins. If the last 
            polygon has frame_end, then that time is used as the end time of the
            last polygon. Otherwise, the last time in t_edges is found by taking
            the frame_time of the last lasso and extending it by the frame_time 
            difference between the second to last and last lassos.
    """

    time_keys = {'frame_time':'%Y-%m-%dT%H:%M:%S',
                 'frame_end':'%Y-%m-%dT%H:%M:%S'}
    polys_log = read_polys(filename, 
                       sort_key = 'frame_time', time_keys=time_keys)

    if 'lon_verts' in polys_log[0]:
        lon_name, lat_name = 'lon_verts', 'lat_verts'
    else:
        "Assuming lon, lat data are in the x_verts, y_verts polygon log entries"
        lon_name, lat_name = 'x_verts', 'y_verts'
    
    flash_stat_polys = tuple(list(zip(p[lon_name], p[lat_name])) for p in polys_log)
    flash_stat_left_time_edges = tuple(p['frame_time'] for p in polys_log)
    if 'frame_end' not in polys_log[-1]:
        dt = flash_stat_left_time_edges[-1] - flash_stat_left_time_edges[-2]
        t_edges = flash_stat_left_time_edges + (flash_stat_left_time_edges[-1] + dt, )
    else:
        flash_stat_right_time_edges = tuple(p['frame_end'] for p in polys_log)
        t_edges = flash_stat_left_time_edges + (flash_stat_right_time_edges[-1],)
        
    if lat_lon_filter is not None:
        flash_stat_polys = apply_filter_poly(lat_lon_filter, polys_log, flash_stat_polys)
    
    return (flash_stat_polys, t_edges) def polys_to_bounding_box(polys):
    """ Given a sequence of (x,y) tuples in polys, return xspan, yspan
        where xspan is a tuple of left, right bounds
        and   yspan is a tuple of bottom, top bounds 
        that encompass all polygons
    """
    all_verts = np.asarray([p for p in itertools.chain(*polys)])
    sw_corner = all_verts.min(axis=0)
    ne_corner = all_verts.max(axis=0)
    lons = (sw_corner[0], ne_corner[0])
    lats = (sw_corner[1], ne_corner[1])
    return lons, lats def h5_files_from_standard_path(path_to_sort_results, date_min, date_max):
    """ Path contains a set of directories with h5 files as follows
        path/h5_files/yyyy/Mon/dd/LYLOUT*.h5
        where Mon is the three-letter month name (e.g., Dec).
        This is the directory structure produced by many of the lmatools examples.
        date_min and date_max are datetime objects given the range of dates to examine.
        Return a list of unique HDF5 filenames (sorted by filename) on that path.
    
        Right now, date_min and date_max can only be used to span a single day 
        boundary.
    """
    globber_start = path_to_sort_results+'/h5_files/{0}/*.h5'.format(date_min.strftime('%Y/%b/%d'))
    globber_end = path_to_sort_results+'/h5_files/{0}/*.h5'.format(date_max.strftime('%Y/%b/%d'))
    # print globber_start
    # print globber_end
    h5_filenames = glob.glob(globber_start)
    h5_filenames += glob.glob(globber_end)
    # Keep only the unique filenames
    h5_filenames = set(h5_filenames)

    reduced_h5 = []
    for h5i in h5_filenames:
        h5time, h5endtime = parse_lma_h5_filename(h5i)
        if ( (h5time - date_min).total_seconds() > -60.0 ) & ( (h5time - date_max).total_seconds() < 60.0 ):
            reduced_h5.append(h5i)
    h5_filenames = reduced_h5
    h5_filenames.sort()
    assert len(h5_filenames) > 0
    return h5_filenames def nc_files_from_standard_path(path_to_sort_results, field_file_code, date_min, date_max):
    filenames = glob.glob(path_to_sort_results+'/grid_files/{1}/*{2}.nc'.format('unused',date_min.strftime('%Y/%b/%d'), field_file_code))
    filenames += glob.glob(path_to_sort_results+'/grid_files/{1}/*{2}.nc'.format('unused',date_max.strftime('%Y/%b/%d'), field_file_code))
    # Keep only the unique filenames
    filenames = set(f for f in filenames if not (f.find('_3d.nc') >=0))

    # nc_names = glob.glob(grid_dir+'/20%s/*.nc' %(date.strftime('%y/%b/%d')))
    # nc_names_3d = glob.glob(grid_dir+'/20%s/*_3d.nc' %(date.strftime('%y/%b/%d')))
    # nc_names_2d = list(set(nc_names) - set(nc_names_3d))

    reduced_grids = []
    for gi in filenames:
        g_parts = os.path.split(gi)[-1].split('_')[1:3]
        g_duration_part = float(os.path.split(gi)[-1].split('_')[3])
        gtime = datetime.strptime('_'.join(g_parts), '%Y%m%d_%H%M%S')
        if ( ( (gtime - date_min).total_seconds() > -g_duration_part ) &
             ( (gtime - date_max).total_seconds() <  g_duration_part ) ):
            reduced_grids.append(gi)
    filenames = reduced_grids
    filenames.sort()
    # print filenames
    assert len(filenames) > 0
    return filenames class TimeSeriesPolygonFlashSubset(TimeSeriesFlashSubset):
    def __init__(self, *args, **kwargs):
        # could also accept coord_names and time_key kwargs here, but those
        # should be standardized in the lmatools h5 format, so we hard code
        # them below
        
        # This is here because the PolygonLassoFilter imports matplotlib
        # and we don't want other imports from this module that don't
        # need matplotlib to trigger a matplotlib import
        from lmatools.lasso.energy_stats import TimeSeriesPolygonLassoFilter 
        
        self.polys = kwargs.pop('polys', [])
        self.t_edges_polys = kwargs.pop('t_edges_polys', [])

        super(TimeSeriesPolygonFlashSubset, self).__init__(*args, **kwargs)
        
        # strictly speaking, we don't even need the time series part; that's been done
        # and we're just lassoin' the points. 
        # But, this code is known to work, so we just reuse it here.
        self.fl_lassos = TimeSeriesPolygonLassoFilter(coord_names=('init_lon', 'init_lat'), time_key='start',
                                         time_edges=self.t_edges_polys, polys=self.polys, basedate=self.base_date )
        self.ev_lassos = TimeSeriesPolygonLassoFilter(coord_names=('lon', 'lat'), time_key='time',
                                         time_edges=self.t_edges_polys, polys=self.polys, basedate=self.base_date)
        self.grid_lassos = TimeSeriesPolygonLassoFilter(coord_names=('lon', 'lat'), time_key='t',
                                         time_edges=self.t_edges_polys, polys=self.polys, basedate=self.base_date)
    
    def gen_chopped_events_flashes(self, *args, **kwargs):
        parent = super(TimeSeriesPolygonFlashSubset, self)
        for ch_ev_series, ch_fl_series in parent.gen_chopped_events_flashes(*args, **kwargs):
            # This gives a time series for each HDF5 LMA file. Next loop
            # over each chopped time series window.
            # Apply polygon filter to time series created by superclass, yield chopped events and flashes
            lassoed_ev = [ch_ev[self.ev_lassos.filter_mask(ch_ev)] for ch_ev in ch_ev_series]
            lassoed_fl = [ch_fl[self.fl_lassos.filter_mask(ch_fl)] for ch_fl in ch_fl_series]
            yield lassoed_ev, lassoed_fl class rho_retrieve(object):
    def __init__(self, area, d, zinit, separation, constant=False, arbitrary_rho=0.4e-9):
    def rho_br(self, z, separation):
    def energy_estimation(self, rho, d, area):
    def calculate(self):
import argparse
import os, sys, errno
import numpy as np
import matplotlib.pyplot as plt
    import xray as xr
    import pandas as pd
    import matplotlib.pyplot as plt
    import matplotlib
    import matplotlib.colors
from datetime import datetime, timedelta
from numpy.lib.recfunctions import append_fields, stack_arrays
from lmatools.grid.grid_collection import LMAgridFileCollection
from lmatools.flash_stats import plot_energy_from_area_histogram, get_energy_spectrum_bins,  bin_center, plot_energies
from lmatools.lasso.energy_stats import flash_size_stats, plot_flash_stat_time_series, plot_energy_stats, plot_tot_energy_stats
from lmatools.lasso.length_stats import FractalLengthProfileCalculator
from lmatools.lasso.cell_lasso_util import read_poly_log_file, polys_to_bounding_box, h5_files_from_standard_path, nc_files_from_standard_path
from lmatools.lasso.cell_lasso_timeseries import TimeSeriesPolygonFlashSubset
from lmatools.lasso import empirical_charge_density as cd ###NEW
    from scipy.stats import scoreatpercentile
    from matplotlib.cm import get_cmap
        fl_good = np.ones(flashes.shape, dtype=bool)
        
        for k, (v_min, v_max) in bounds.items():
            if k in events.dtype.names:
                ev_good &= (events[k] >= v_min) & (events[k] <= v_max)
            if k in flashes.dtype.names:
                fl_good &= (flashes[k] >= v_min) & (flashes[k] <= v_max)
        events_filt, flashes_filt = (events[ev_good], flashes[fl_good])
        
        yield (events_filt, flashes_filt)

# used to remove flashes not meeting certain critera
bounds={'area':(min_area, max_area),
        'ctr_alt':(min_alt,max_alt),
        'total_energy':(min_energy, max_energy)
}
filtered_time_series = gen_filtered_time_series(events_series, flashes_series, bounds)
events_series, flashes_series = zip(*filtered_time_series)

# =====
# plot the raw flash locations
# =====
class FlashCentroidPlotter(object):
    def __init__(self, t_sec_min, t_sec_max):
        self.t_sec_min, self.t_sec_max = t_sec_min, t_sec_max
        self.fig = plt.figure(figsize=(11,11))
        self.evax = self.fig.add_subplot(111)
        self.evax.set_ylabel('Latitude (deg)')
        self.evax.set_xlabel('Longitude (deg)')
        self.evax.set_title('Flash event-weighted centroid, \n{0} to {1} seconds'.format(
            t_sec_min, t_sec_max
            # date_min.strftime('%Y-%m-%d %H%M:%S'),
            # date_max.strftime('%Y-%m-%d %H%M:%S')
        ))
        self.colorbar = None
        self.scatters = []

    def plot(self, flashes):
        scart = self.evax.scatter(flashes['init_lon'], flashes['init_lat'],
                      c=flashes['start'], s=4, cmap='viridis',
                      vmin=self.t_sec_min, vmax=self.t_sec_max, edgecolor='none')
        if self.colorbar is None:
            self.colorbar = plt.colorbar(scart, ax=self.evax)
        self.scatters.append(scart)

flash_location_plotter = FlashCentroidPlotter(min(flashes_in_poly.t_edges_seconds),
                                              max(flashes_in_poly.t_edges_seconds))
for flashes in flashes_series:
    flash_location_plotter.plot(flashes)
flash_ctr_filename = 'flash_ctr_{0}_{1}.png'.format(t_start.strftime('%y%m%d%H%M%S'),
                                                    t_end.strftime('%y%m%d%H%M%S'))
flash_location_plotter.fig.savefig(os.path.join(outdir, flash_ctr_filename))


# =======================================
# Loop over each window in the time series and calculate some
# aggregate flash statistical properties
# =======================================

# Set up fractal length calcuations and channel height profiles
if do_3d:
    D = 1.5
    b_s = 200.0
    max_alt, d_alt = 20.0, 0.5
    alt_bins = np.arange(0.0,max_alt+d_alt, d_alt)
    length_profiler = FractalLengthProfileCalculator(D, b_s, alt_bins)
else:
    length_profiler = None

def write_flash_volume(outfile, field):
    dat = pd.DataFrame(field).T
    dat.to_csv(outfile)

def gen_flash_summary_time_series(events_series, flashes_series, length_profiler):
    for events, flashes in zip(events_series, flashes_series):
        # reduce all flashes in this time interval to representative moments
        size_stats = flash_size_stats(flashes)
        # for each flash, get a dictionary with 2D and 3D fractal lengths.
        # also includes the volume and point triangulation data.
        if length_profiler is not None:
            per_flash_data, IC_profile, CG_profile = length_profiler.process_events_flashes(events, flashes)
        else:
            IC_profile, CG_profile = None, None
        yield size_stats, IC_profile, CG_profile


time_series = gen_flash_summary_time_series(events_series, flashes_series,
                                            length_profiler)
size_stats, IC_profiles, CG_profiles = zip(*time_series)

size_stats = stack_arrays(size_stats)
iso_start, iso_end = flashes_in_poly.t_edges_to_isoformat(as_start_end=True)
# assume here that iso_start and iso_end are identical-length strings,
# as they should be if the iso format is worth anything.
size_stats = append_fields(size_stats, ('start_isoformat','end_isoformat'),
                                         data=(iso_start, iso_end), usemask=False)                                      
# =====
# Write flash size stats data (see harvest_flash_timeseries) and plot moment time series
# =====
def write_size_stat_data(outfile, size_stats):
    stat_keys = ('start_isoformat','end_isoformat', 'number', 'mean', 'variance',
                 'skewness', 'kurtosis', 'energy', 'energy_per_flash','total_energy','specific_energy')
    header = "# start_isoformat, end_isoformat, number, mean, variance, skewness, kurtosis, energy, energy_per_flash, total_energy, specific_energy\n"
    line_template = "{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}\n"
    f = open(outfile, 'w')
    #for start_t, end_t, stats in zip(starts, ends, size_stats):
    f.write(header)
    for stats in size_stats:
        stat_vals = [stats[k] for k in stat_keys]
        line = line_template.format(*stat_vals)
        f.write(line)
    f.close()

stats_filename = os.path.join(outdir,'flash_stats.csv')
stats_figure = os.path.join(outdir,'flash_stats_{0}_{1}.pdf'.format(t_start.strftime('%y%m%d%H%M%S'),
                                                    t_end.strftime('%y%m%d%H%M%S')))
write_size_stat_data(stats_filename, size_stats)
fig = plot_flash_stat_time_series(flashes_in_poly.base_date, flashes_in_poly.t_edges, size_stats)
fig.savefig(stats_figure)


energy_fig = plot_energy_stats(size_stats['specific_energy'], flashes_in_poly.base_date, flashes_in_poly.t_edges, outdir)