MPLtools.py

 # -*- coding: utf-8 -*-
"""
mpltools.py
A bag of tools to be used in processing and interpreting MPL class data
collected by miniMPL
Created on Wed Apr 24 12:08:57 2013

@author: Paul Cottle

"""
import os,sys,site

from Tkinter import Tk
import tkFileDialog
import re
import numpy as np
import array, struct
import pandas as pan
import datetime
from scipy import constants as const    
from copy import deepcopy
from scipy.interpolate import interp1d
from scipy.ndimage.filters import generic_filter as genfilt
from collections import OrderedDict 
import h5py 

if sys.platform == 'win32': 
    from matplotlib import pyplot as plt
else:
    import matplotlib
    matplotlib.use('Agg')
    from matplotlib import pyplot as plt

class MPL:
    """
    This is a class type generated by unpacking a binary file generated by
    the mini-MPL lidar
    
    It includes two subclasses: header and data
    The metadata in the header is described in the MPL manual pp 37-38
    The data consists of a 2-D array of floats separated into channels
    
    copol = data measured in laser polarization
    crosspol = data measured in cross polarization
    
    """
        
    def __init__(self,filename=[]):
        
        self.header = None #slot for lidar header data
        self.data = None #slot for lidar raw data array
        self.rsq = None  #slot for range corrected, background subtracted data
        self.NRB = None  #slot for Normalized Relative Backscatter array 
        self.depolrat = None  #slot fo depol ratio array
        self.sigma = None   #slot for standard deviation data
        self.SNR = None     #slot for SNR data
        self.backscatter = None #slot for corrected backscatter array
        self.extinction = None  #slot for extinction array
        self.scenepanel = None  #slot for panel containing scene analysis features
        

    def copy(self):
        #currnetly not working
        
        return MPL(deepcopy(self))
    
    def append(self,MPLnew):
        
        if self.header is None:
            self.header = MPLnew.header
        else:
            maxprofnum = self.header['profnum'][-1]
            MPLnew.header['profnum'] = [p+maxprofnum for p in MPLnew.header['profnum']]
            
            self.header = self.header.append(MPLnew.header)
        
        numchans=self.header['numchans'][0]
        
        if self.data is None:
            self.data = MPLnew.data
        else:            
            for n in range(numchans):
                self.data[n] = self.data[n].append(MPLnew.data[n])
        
        if MPLnew.rsq is not None:
            if self.rsq is None:
                self.rsq = MPLnew.rsq
            else:
                for n in range(numchans):
                    self.rsq[n] = self.rsq[n].append(MPLnew.rsq[n])

        if MPLnew.NRB is not None:
            if self.NRB is None:
                self.NRB = MPLnew.NRB
            else:
                for n in range(numchans):
                    self.NRB[n] = self.NRB[n].append(MPLnew.NRB[n])
        
        if MPLnew.depolrat is not None:
            if self.depolrat is None:
                self.depolrat = MPLnew.depolrat
            else:
                self.depolrat[0] = self.depolrat[0].append(MPLnew.depolrat[0])
        
        if MPLnew.backscatter is not None:
            if self.backscatter is None:
                self.backscatter = MPLnew.backscatter
            else:
                for n in range(numchans):
                    self.backscatter[n] = self.backscatter[n].append(MPLnew.backscatter[n])
        
        if MPLnew.extinction is not None:
            if self.extinction is None:
                self.extinction = MPLnew.extinction
            else:
                for n in range(numchans):
                    self.extinction[n] = self.extinction[n].append(MPLnew.extinction[n])
        
        if MPLnew.scenepanel is not None:
            if self.scenepanel is None:
                self.scenepanel = MPLnew.scenepanel
            else:
                for n in range(len(self.scenepanel)):
                    for i in self.scenepanel[n].items:
                        self.scenepanel[n][i] = self.scenepanel[n][i].append(MPLnew.scenepanel[n][i])
        
        if MPLnew.sigma is not None:
            if self.sigma is None:
                self.sigma = MPLnew.sigma
            else:
                for key in self.sigma:
                    self.sigma[key]=self.sigma[key].append(MPLnew.sigma[key])
                    
        if MPLnew.SNR is not None:
            if self.SNR is None:
                self.SNR = MPLnew.SNR
            else:
                for key in self.SNR:
                    self.SNR[key]=self.SNR[key].append(MPLnew.SNR[key])            
        
        return self
    
    def fromMPL(self, filename):
        with open(filename,'rb') as binfile:
            profdat_copol = OrderedDict()
            profdat_crosspol = OrderedDict()
            header = OrderedDict()
            
            profnum = 0
            
            while True:
                try:
                    intarray16 = array.array('H')
                    intarray32 = array.array('I') # L is 8 byte on Xenon
                    floatarray = array.array('f')
                    byte_array = array.array('B')
                    copolvals = array.array('f')  
                    crosspolvals = array.array('f')  
                    
                    intarray16.fromfile(binfile, 8)
                    intarray32.fromfile(binfile, 8)
                    floatarray.fromfile(binfile, 2)
                    intarray16.fromfile(binfile, 1)
                    intarray32.fromfile(binfile, 1)
                    floatarray.fromfile(binfile, 2)
                    intarray16.fromfile(binfile, 3)
                    floatarray.fromfile(binfile, 8)
                    byte_array.fromfile(binfile, 2)
                    floatarray.fromfile(binfile, 2)
                    byte_array.fromfile(binfile, 1)
                    intarray16.fromfile(binfile, 1)
                    byte_array.fromfile(binfile, 1)
                    intarray16.fromfile(binfile, 3)
                    
                    headerdat = {}
                    headerdat['unitnum'] = intarray16[0]
                    headerdat['version'] = intarray16[1]
                    year = intarray16[2]
                    month = intarray16[3]
                    day = intarray16[4]
                    hour = intarray16[5]
                    minute = intarray16[6]
                    second = intarray16[7]
                                       
                    try:
                        dt = datetime.datetime(year,month,day,hour,minute,second)
                    except ValueError:
                        print "Error extracting data from {0}".format(filename)                        
                        break
            
                    headerdat['shotsum'] = intarray32[0]  #total number of shots collected per profile
                    headerdat['trigfreq'] = intarray32[1] #laser trigger frequency (usually 2500 Hz)
                    headerdat['energy'] = intarray32[2]/1000.0  #mean of laser energy monitor in uJ                      
                    headerdat['temp_0'] = intarray32[3]/100.0  #mean of A/D#0 readings*100
                    headerdat['temp_1'] = intarray32[4]/100.0  #mean of A/D#1 readings*100
                    headerdat['temp_2'] = intarray32[5]/100.0  #mean of A/D#2 readings*100
                    headerdat['temp_3'] = intarray32[6]/100.0  #mean of A/D#3 readings*100
                    headerdat['temp_4'] = intarray32[7]/100.0  #mean of A/D#4 readings*100
                    
                    headerdat['bg_avg1'] = floatarray[0] #mean background signal value for channel 1
                    headerdat['bg_std1'] = floatarray[1] #standard deviation of backgruond signal for channel 1
            
                    headerdat['numchans'] = intarray16[8] #number of channels
                    headerdat['numbins'] = intarray32[8] #total number of bins per channel
            
                    headerdat['bintime'] = floatarray[2]  #bin width in seconds
                    
                    headerdat['rangecal'] = floatarray[3]/1000.0 #range offset in km, default is 0
            
                    headerdat['databins'] = intarray16[9]  #number of bins not including those used for background
                    headerdat['scanflag'] = intarray16[10]  #0: no scanner, 1: scanner
                    headerdat['backbins'] = intarray16[11]  #number of background bins
            
                    headerdat['az'] = floatarray[4]  #scanner azimuth angle
                    headerdat['el'] = floatarray[5]  #scanner elevation angle
                    headerdat['deg'] = floatarray[6] #compass degrees (currently unused)
                    headerdat['pvolt0'] = floatarray[7] #currently unused
                    headerdat['pvolt1'] = floatarray[8] #currently unused
                    headerdat['gpslat'] = floatarray[9] #GPS latitude in decimal degreees (-999.0 if no GPS)
                    headerdat['gpslon'] = floatarray[10] #GPS longitude in decimal degrees (-999.0 if no GPS)
                    headerdat['cloudbase'] = floatarray[11] #cloud base height in [m]
            
                    headerdat['baddat'] = byte_array[0]  #0: good data, 1: bad data
                    headerdat['version'] = byte_array[1] #version of file format.  current version is 5
            
                    headerdat['bg_avg2'] = floatarray[12] #mean background signal for channel 2
                    headerdat['bg_std2'] = floatarray[13] #mean background standard deviation for channel 2
            
                    headerdat['mcs'] = byte_array[2]  #MCS mode register  Bit#7: 0-normal, 1-polarization
                                                 #Bit#6-5: polarization toggling: 00-linear polarizer control
                                                 #01-toggling pol control, 10-toggling pol control 11-circular pol control
            
                    headerdat['firstbin'] = intarray16[12]  #bin # of first return data
                    headerdat['systype'] = byte_array[3]   #0: standard MPL, 1: mini MPL
                    headerdat['syncrate'] = intarray16[13]  #mini-MPL only, sync pulses seen per second
                    headerdat['firstback'] = intarray16[14] #mini-MPL only, first bin used for background calcs
                    headerdat['headersize2'] = intarray16[15] #size of additional header data (currently unused)
                    headerdat['profnum'] = profnum
                    profnum += 1
                    
                    if headerdat['headersize2'] > 128:
                        byte_array.fromfile(binfile, 1)
                        floatarray.fromfile(binfile, 6)
                        intarray16.fromfile(binfile, 1)
                        floatarray.fromfile(binfile, 2)
                    
                        headerdat['Weatherstat'] = byte_array[4]   #0:weatehr station not used, 1:weather station used
                        headerdat['Int_temp'] = floatarray[14]  #Temperature inside in deg. Celsius
                        headerdat['Ext_temp'] = floatarray[15]  #Temp. outside
                        headerdat['Int_humid'] = floatarray[16] #Humidity inside in %
                        headerdat['Ext_humid'] = floatarray[17] #Hum. outside
                        headerdat['Dewpoint'] = floatarray[18]  #dewpoint in deg. Celsius
                        headerdat['Wnd_spd'] = floatarray[19]  #wind speed in km/h
                        headerdat['Wnd_dir'] = intarray16[16]  #wind direction in deg.
                        headerdat['Press'] = floatarray[20]  #Barometric pressure in hPa
                        headerdat['Rain'] = floatarray[21]  #Rain rate in mm/hr

                    numbins = headerdat['numbins']
                    numchans = headerdat['numchans'] 
                    altstep = headerdat['bintime']*const.c/2000.0 #altitude step in km
                    maxalt = numbins*altstep
                    firstbin=headerdat['firstbin']
                    minalt = headerdat['rangecal']+firstbin*altstep
                    altrange = np.arange(minalt,maxalt,altstep,dtype='float')
                    if len(altrange) != numbins:
                        altrange=altrange[:numbins]
                    
                    if numchans == 2:
                        crosspolvals.fromfile(binfile, numbins) 
                        temp = np.array(crosspolvals)                       
                        profdat_crosspol[dt] = temp[firstbin:]
                        copolvals.fromfile(binfile, numbins) 
                        temp = np.array(copolvals)
                        profdat_copol[dt] = temp[firstbin:]
                    else:
                        raise ValueError('Wrong number of channels')
                    
                    header[dt] = headerdat
                
                except EOFError:
                    break
            
        df_copol = pan.DataFrame.from_dict(profdat_copol,orient = 'index')
        df_copol.columns = altrange
        #minimum usable altitude is 150m for overlap reasons
        if altrange[0]<=0.150:
            df_copol=df_copol.loc[:,0.150:]
        df_crosspol = pan.DataFrame.from_dict(profdat_crosspol,orient = 'index')
        df_crosspol.columns = altrange
        if altrange[0]<=0.150:
            df_crosspol=df_crosspol.loc[:,0.150:]
        df_header = pan.DataFrame.from_dict(header, orient = 'index')
                
        self.data = [df_copol, df_crosspol]
        
        
        self.header = df_header

        return self        
    
    def fromHDF(self, filename, verbose = False):
                
        copoldat = pan.read_hdf(filename,'copol_raw')
        crosspoldat = pan.read_hdf(filename,'crosspol_raw')
        header = pan.read_hdf(filename,'header')
            
        self.data = [copoldat,crosspoldat]
        self.header = header
        
        try:
            copoldat_rsq = pan.read_hdf(filename,'copol_rsq')
            crosspoldat_rsq = pan.read_hdf(filename,'crosspol_rsq')
            self.rsq = [copoldat_rsq,crosspoldat_rsq]
        except KeyError:
            if verbose:
                print "Warning:  No Range-squared file"
        
        try:
            copoldat_NRB = pan.read_hdf(filename,'copol_NRB')
            crosspoldat_NRB = pan.read_hdf(filename,'crosspol_NRB')
            self.NRB = [copoldat_NRB,crosspoldat_NRB]
        except KeyError:
            if verbose:
                print "Warning: No NRB file"
        
        try:
            self.depolrat = [pan.read_hdf(filename,'depolrat')]
        except KeyError:
            if verbose:
                print "Warning: No Depol Ratio file"
        
        try:
            copoldat_back=pan.read_hdf(filename,'copol_backscatter')
            self.backscatter=[copoldat_back]
        except KeyError:
            if verbose:
                print "Warning: No Copol Backscatter file"
        try:
            crosspoldat_back=pan.read_hdf(filename,'crosspol_backscatter')
            self.backscatter.append(crosspoldat_back)
        except KeyError:
            if verbose:
                print "Warning: No Crosspol Backscatter file"
        
        try:
            copoldat_ext=pan.read_hdf(filename,'copol_extinction')
            self.extinction=[copoldat_ext]
        except KeyError:
            if verbose:
                print "Warning: No Copol Extinction file"
        try:
            crosspoldat_ext=pan.read_hdf(filename,'crosspol_extinction')
            self.extinction.append(crosspoldat_ext)
        except KeyError:
            if verbose:
                print "Warning: No Crosspol Extinction file"
        
        try:
            scenedat=pan.read_hdf(filename,'scenepanel')
            self.scenepanel=[scenedat]
        except KeyError:
            if verbose:
                print "Warning: No Scene Analysis file"
                
        sigmadict={}
        try:            
            tempsigma_copol=pan.read_hdf(filename,'sigma_copol_data')
            tempsigma_depol=pan.read_hdf(filename,'sigma_crosspol_data')
            sigmadict['data']=[tempsigma_copol,tempsigma_depol]
        except KeyError:
            if verbose:
                print "Warning: No sigma-data file"
        try:            
            tempsigma_copol=pan.read_hdf(filename,'sigma_copol_rsq')
            tempsigma_depol=pan.read_hdf(filename,'sigma_crosspol_rsq')
            sigmadict['rsq']=[tempsigma_copol,tempsigma_depol]
        except KeyError:
            if verbose:
                print "Warning: No sigma-rsq  file"
        try:            
            tempsigma_copol=pan.read_hdf(filename,'sigma_copol_NRB')
            tempsigma_depol=pan.read_hdf(filename,'sigma_crosspol_NRB')
            sigmadict['NRB']=[tempsigma_copol,tempsigma_depol]
        except KeyError:
            if verbose:
                print "Warning: No sigma-NRB file"                
        try:            
            tempsigma=pan.read_hdf(filename,'sigma_depolrat')
            sigmadict['depolrat']=[tempsigma]
        except KeyError:
            if verbose:
                print "Warning: No sigma-depolrat file"
        try:            
            tempsigma_copol=pan.read_hdf(filename,'sigma_copol_backscatter')
            sigmadict['backscatter']=[tempsigma_copol]
        except KeyError:
            if verbose:
                print "Warning: No sigma-copol-backscatter file"
        try:        
            tempsigma_crosspol=pan.read_hdf(filename,'sigma_crosspol_backscatter')
            sigmadict['backscatter'].append(tempsigma_crosspol)
        except KeyError:
            if verbose:
                print "Warning: No sigma-crosspol-backscatter file"
        try:            
            tempsigma_copol=pan.read_hdf(filename,'sigma_copol_extinction')
            sigmadict['extinction']=[tempsigma_copol]
        except KeyError:
            if verbose:
                print "Warning: No sigma-copol-extinction file"
        try:
            temp_sigma_crosspol=pan.read_hdf(filename,'sigma_crosspol_extinction')
            sigmadict['extinction'].append(temp_sigma_crosspol)
        except KeyError:
            if verbose:
                print "Warning: No sigma-crosspol-extinction file"
        
        if sigmadict:
            self.sigma=sigmadict
            
        SNRdict={}
        try:            
            tempSNR_copol=pan.read_hdf(filename,'SNR_copol_data')
            tempSNR_depol=pan.read_hdf(filename,'SNR_crosspol_data')
            SNRdict['data']=[tempSNR_copol,tempSNR_depol]
        except KeyError:
            if verbose:
                print "Warning: No SNR-data file"
        try:            
            tempSNR_copol=pan.read_hdf(filename,'SNR_copol_rsq')
            tempSNR_depol=pan.read_hdf(filename,'SNR_crosspol_rsq')
            SNRdict['rsq']=[tempSNR_copol,tempSNR_depol]
        except KeyError:
            if verbose:
                print "Warning: No SNR-rsq  file"
        try:            
            tempSNR_copol=pan.read_hdf(filename,'SNR_copol_NRB')
            tempSNR_depol=pan.read_hdf(filename,'SNR_crosspol_NRB')
            SNRdict['NRB']=[tempSNR_copol,tempSNR_depol]
        except KeyError:
            if verbose:
                print "Warning: No SNR-NRB file"                
        try:            
            tempSNR=pan.read_hdf(filename,'SNR_depolrat')
            SNRdict['depolrat']=[tempSNR]
        except KeyError:
            if verbose:
                print "Warning: No SNR-depolrat file"
        try:            
            tempSNR_copol=pan.read_hdf(filename,'SNR_copol_backscatter')
            SNRdict['backscatter']=[tempSNR_copol]
        except KeyError:
            if verbose:
                print "Warning: No SNR-copol-backscatter file"
        try:
            tempSNR_crosspol=pan.read_hdf(filename,'SNR_crosspol_backscatter')
            SNRdict['backscatter'].append(tempSNR_crosspol)
        except KeyError:
            if verbose:
                print "Warning: No SNR-crosspol-backscatter file"
        try:            
            tempSNR_copol=pan.read_hdf(filename,'SNR_copol_extinction')
            SNRdict['extinction']=[tempSNR_copol]
        except KeyError:
            if verbose:
                print "Warning: No SNR-copol-extinction file"
        try:
            temp_SNR_crosspol=pan.read_hdf(filename,'SNR_crosspol_extinction')
            SNRdict['extinction'].append(temp_SNR_crosspol)
        except KeyError:
            if verbose:
                print "Warning: No SNR-extinction file"
                
        if SNRdict:
            self.SNR=SNRdict
            
        return self
        
    def save_to_HDF(self, filename):
        
        store = pan.HDFStore(filename)
        
        store['header'] = self.header
        
        df_copol = self.data[0]
        df_crosspol = self.data[1]
        store['copol_raw'] = df_copol
        store['crosspol_raw'] = df_crosspol
        
        if self.rsq is not None:
            df_copol = self.rsq[0]
            df_crosspol = self.rsq[1]
            store['copol_rsq'] = df_copol
            store['crosspol_rsq'] = df_crosspol
        
        if self.NRB is not None:
            df_copol = self.NRB[0]
            df_crosspol = self.NRB[1]
            store['copol_NRB'] = df_copol
            store['crosspol_NRB'] = df_crosspol
        
        if self.depolrat is not None:
            df_depolrat = self.depolrat[0]
            store['depolrat'] = df_depolrat
        
        if self.backscatter is not None:
            df_backscatter_copol = self.backscatter[0]
            store['copol_backscatter'] = df_backscatter_copol
            if len(self.backscatter)>1:
                df_backscatter_crosspol = self.backscatter[1]
                store['crosspol_backscatter'] = df_backscatter_crosspol
        
        if self.extinction is not None:
            df_extinction_copol = self.extinction[0]
            store['copol_extinction'] = df_extinction_copol
            if len(self.extinction)>1:
                df_extinction_crosspol = self.extinction[1]
                store['crosspol_extinction'] = df_extinction_crosspol
        
        if self.scenepanel is not None:
            scenepanel = self.scenepanel[0]
            store['scenepanel'] = scenepanel
            
        if self.SNR is not None:
            for k,v in self.SNR.iteritems():
                if len(v)==2:
                    savename1='SNR_copol_{0}'.format(k)
                    savename2='SNR_crosspol_{0}'.format(k)
                    tempdf_copol=v[0]
                    tempdf_crosspol=v[1]
                    store[savename1]=tempdf_copol
                    store[savename2]=tempdf_crosspol
                else:
                    savename='SNR_{0}'.format(k)
                    tempdf_copol=v[0]
                    store[savename]=tempdf_copol
        
        if self.sigma is not None:
            for k,v in self.sigma.iteritems():
                if len(v)==2:
                    savename1='sigma_copol_{0}'.format(k)
                    savename2='sigma_crosspol_{0}'.format(k)
                    tempdf_copol=v[0]
                    tempdf_crosspol=v[1]
                    store[savename1]=tempdf_copol
                    store[savename2]=tempdf_crosspol
                else:
                    savename='sigma_{0}'.format(k)
                    tempdf_copol=v[0]
                    store[savename]=tempdf_copol
        
        store.close()
    
    def save_to_IDL(self,filename):
        #needs updating to include new attributes SNR and sigma
        def datetime_to_epoch(d):
            epoch=np.datetime64(datetime.datetime(1970,1,1))
            t=(d-epoch)/np.timedelta64(1,'s')
            return t
            
        with h5py.File(filename,'w') as f:        
            headergrp=f.create_group('header')
            headergrp.create_dataset('keys',data=np.array(self.header.columns.values,dtype='str'))
            datetime_index=[datetime_to_epoch(val) for val in self.header.index.values]
            headergrp.create_dataset('timetag',data=datetime_index)
            headergrp.create_dataset('values',data=self.header.values)
            
            if self.data:
                df_copol = self.data[0]
                df_crosspol = self.data[1]
                tempgrp=f.create_group('raw_data')
                tempcopol=tempgrp.create_group('copol')
                datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                tempcopol.create_dataset('timetag',data=datetime_index)
                tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                tempcopol.create_dataset('values',data=df_copol.values)
                tempcrosspol=tempgrp.create_group('crosspol')
                tempcrosspol.create_dataset('timetag',data=datetime_index)
                tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                tempcrosspol.create_dataset('values',data=df_crosspol.values)
            if self.rsq:
                df_copol = self.rsq[0]
                df_crosspol = self.rsq[1]
                tempgrp=f.create_group('range_squared')
                tempcopol=tempgrp.create_group('copol')
                datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                tempcopol.create_dataset('timetag',data=datetime_index)
                tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                tempcopol.create_dataset('values',data=df_copol.values)
                tempcrosspol=tempgrp.create_group('crosspol')
                tempcrosspol.create_dataset('timetag',data=datetime_index)
                tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                tempcrosspol.create_dataset('values',data=df_crosspol.values)
            if self.NRB:
                df_copol = self.NRB[0]
                df_crosspol = self.NRB[1]
                tempgrp=f.create_group('NRB')
                tempcopol=tempgrp.create_group('copol')
                datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                tempcopol.create_dataset('timetag',data=datetime_index)
                tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                tempcopol.create_dataset('values',data=df_copol.values)
                tempcrosspol=tempgrp.create_group('crosspol')
                tempcrosspol.create_dataset('timetag',data=datetime_index)
                tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                tempcrosspol.create_dataset('values',data=df_crosspol.values)
            if self.depolrat:
                df = self.depolrat[0]
                tempgrp=f.create_group('depol_ratio')
                datetime_index=[datetime_to_epoch(d) for d in df.index.values]
                tempgrp.create_dataset('timetag',data=datetime_index)
                tempgrp.create_dataset('altitude',data=df.columns.values)
                tempgrp.create_dataset('values',data=df.values)
            if self.SNR:
                tempgrp=f.create_group('SNR')
                for k,v in self.SNR.iteritems():
                    tempsubgrp=tempgrp.create_group(k)
                    if len(v)==2:
                        df_copol=v[0]
                        df_crosspol=v[1]
                        tempcopol=tempsubgrp.create_group('copol')
                        datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                        tempcopol.create_dataset('timetag',data=datetime_index)
                        tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                        tempcopol.create_dataset('values',data=df_copol.values)
                        tempcrosspol=tempsubgrp.create_group('crosspol')
                        tempcrosspol.create_dataset('timetag',data=datetime_index)
                        tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                        tempcrosspol.create_dataset('values',data=df_crosspol.values)
                    else:
                        df=v[0]
                        datetime_index=[datetime_to_epoch(d) for d in df.index.values]
                        tempsubgrp.create_dataset('timetag',data=datetime_index)
                        tempsubgrp.create_dataset('altitude',data=df.columns.values)
                        tempsubgrp.create_dataset('values',data=df.values)
            if self.backscatter:
                df_copol = self.backscatter[0]
                df_crosspol = self.backscatter[1]
                tempgrp=f.create_group('Backscatter')
                tempcopol=tempgrp.create_group('copol')
                datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                tempcopol.create_dataset('timetag',data=datetime_index)
                tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                tempcopol.create_dataset('values',data=df_copol.values)
                tempcrosspol=tempgrp.create_group('crosspol')
                tempcrosspol.create_dataset('timetag',data=datetime_index)
                tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                tempcrosspol.create_dataset('values',data=df_crosspol.values)   
            if self.extinction:
                df_copol = self.extinction[0]
                df_crosspol = self.extinction[1]
                tempgrp=f.create_group('Extinction')
                tempcopol=tempgrp.create_group('copol')
                datetime_index=[datetime_to_epoch(d) for d in df_copol.index.values]
                tempcopol.create_dataset('timetag',data=datetime_index)
                tempcopol.create_dataset('altitude',data=df_copol.columns.values)
                tempcopol.create_dataset('values',data=df_copol.values)
                tempcrosspol=tempgrp.create_group('crosspol')
                tempcrosspol.create_dataset('timetag',data=datetime_index)
                tempcrosspol.create_dataset('altitude',data=df_crosspol.columns.values)
                tempcrosspol.create_dataset('values',data=df_crosspol.values)
    
    def alt_resample(self,altrange,sigma_winsize=10,SNR_winsize=10,verbose=False):
        #takes a pandas dataframe generated by mplreader and resamples on regular
        #intervals in altitude and resets the limits of the set
        #note: limits of altrange must be within original limits of altitude data
        
        
        if verbose:
            print 'Altitude step resampling in progress ...'
        
        #resample raw MPL data
        templist=[]
        for dftemp in self.data:            
            templist.append(resample_cols(dftemp,altrange,verbose))
        self.data=templist
        
        #resample range corrected data
        if self.rsq is not None:
            templist=[]
            for dftemp in self.rsq:            
                templist.append(resample_cols(dftemp,altrange,verbose))
            self.rsq=templist
        else:
            if verbose:
                print "No Range-Squared Profiles"
        
        #resample NRB data
        if self.NRB is not None:
            templist=[]
            for dftemp in self.NRB:            
                templist.append(resample_cols(dftemp,altrange,verbose))
            self.NRB=templist
        else:
            if verbose:
                print "No NRB Profiles"
        
        if self.depolrat is not None:       
            templist=[]
            for dftemp in self.depolrat:            
                templist.append(resample_cols(dftemp,altrange,verbose))
            self.depolrat=templist
        else:
            if verbose:
                print "No Depol Ratio Profiles" 

        if self.backscatter is not None:       
            templist=[]
            for dftemp in self.backscatter:
                templist.append(resample_cols(dftemp,altrange,verbose))
            self.backscatter=templist
        else:
            if verbose:
                print "No Backscatter Profiles"

        if self.extinction is not None:       
            templist=[]
            for dftemp in self.extinction:
                templist.append(resample_cols(dftemp,altrange,verbose))                    
            self.extinction=templist
        else:
            if verbose:
                print "No Extinction Profiles"

        if self.scenepanel is not None:
            paneldict={}
            for i in self.scenepanel[0].items:
                dftemp = self.scenepanel[0].loc[i]
                paneldict[i]=resample_cols(dftemp,altrange,verbose,method='ffill')                
            self.scenepanel=[pan.Panel.from_dict(paneldict)]
        else:
            if verbose:
                print "No Scene Analysis"
        
        if self.sigma is not None:
            self.calculate_sigma(winsize=10)
        
        if self.SNR is not None:
            self.calculate_SNR(winsize=10)
            
        if verbose:
            print '... Done!'
                    
        self.header['numbins'] = [len(altrange) for db in self.header['numbins']]
        self.header['databins'] = [db - bb for (db,bb) in zip(self.header['numbins'],self.header['backbins'])]
        self.header['bintime'] = [(altrange[1]-altrange[0])/const.c for bt in self.header['bintime']]
        self.header['firstback'] = [len(altrange)+1 for fb in self.header['firstback']]
        
        return self
    
    def time_resample(self, timestep=None, starttime=None,endtime=None, datamethod = 'mean',
                      sigma_winsize=10,SNR_winsize=10,verbose=False):
        #resamples a pandas dataframe generated by mplreader on a regular timestep
        #and optionally limits it to a preset time range
        #timestep must be in timeSeries period format: numF where num=step size and
        #F = offset alias.  Ex: H = hours, M = minutes, S = seconds, L = millieconds
        
        temphead = {}
        self.header.sort_index(inplace=True)
        if starttime is not None:
            self.header = self.header.loc[self.header.index>=starttime]
        if endtime is not None:
            self.header = self.header.loc[self.header.index<=endtime]
        if timestep is not None:
            for col in self.header:           
                sumcols = ['shotsum']
                firstcols = ['unitnum','version','numchans','scanflag','version','mcs','systype','numchans']
                intmeancols = ['numbins','databins','backbins','firstbin','firstback']
                maxcols = ['baddat']            
                if col in sumcols: headermethod = 'sum'            
                elif col in firstcols: headermethod = 'first'            
                elif col in maxcols:  headermethod = 'max'
                else: headermethod = 'mean'
                            
                temphead[col] = self.header[col].resample(timestep, how = headermethod)            
                
                if (col in intmeancols) or (col in firstcols) or (col in maxcols):
                    try:
                        temphead[col] = temphead[col].astype('int32')
                    except ValueError:
                        whereisna = np.isnan(temphead[col])
                        temphead[col][whereisna] = -999
                        temphead[col] = temphead[col].astype('int32')

        
            self.header = pan.DataFrame(temphead)
        if verbose:
                print 'Time step regularization in progress ...'
        
        templist=[]        
        for dftemp in self.data:                
            if starttime is not None:
                dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
            if endtime is not None:
                dftemp = dftemp.loc[dftemp.index<=endtime]  
            if timestep is not None:
                dftemp = dftemp.resample(timestep, how = datamethod)
            templist.append(dftemp)
        self.data=templist
    
        if self.rsq is not None:
            templist=[]        
            for dftemp in self.rsq:                
                if starttime is not None:
                    dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                if endtime is not None:
                    dftemp = dftemp.loc[dftemp.index<=endtime]  
                if timestep is not None:
                    dftemp = dftemp.resample(timestep, how = datamethod)
                templist.append(dftemp)
            self.rsq=templist

        if self.NRB is not None:
            templist=[]        
            for dftemp in self.NRB:                
                if starttime is not None:
                    dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                if endtime is not None:
                    dftemp = dftemp.loc[dftemp.index<=endtime]  
                if timestep is not None:
                    dftemp = dftemp.resample(timestep, how = datamethod)
                templist.append(dftemp)
            self.NRB=templist

        if self.depolrat is not None:
            templist=[]        
            for dftemp in self.depolrat:                
                if starttime is not None:
                    dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                if endtime is not None:
                    dftemp = dftemp.loc[dftemp.index<=endtime]  
                if timestep is not None:
                    dftemp = dftemp.resample(timestep, how = datamethod)
                templist.append(dftemp)
            self.depolrat=templist
            
        if self.backscatter is not None:
            templist=[]        
            for dftemp in self.backscatter: 
                if starttime is not None:
                    dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                if endtime is not None:
                    dftemp = dftemp.loc[dftemp.index<=endtime]  
                if timestep is not None:
                    dftemp = dftemp.resample(timestep, how = datamethod)
                templist.append(dftemp)
            self.backscatter=templist
        
        if self.extinction is not None:
            templist=[]        
            for dftemp in self.extinction: 
                if starttime is not None:
                    dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                if endtime is not None:
                    dftemp = dftemp.loc[dftemp.index<=endtime]  
                if timestep is not None:
                    dftemp = dftemp.resample(timestep, how = datamethod)
                templist.append(dftemp)
            self.extinction=templist
        
        if self.scenepanel is not None:
            templist=[]
            for paneltemp in self.scenepanel:
                panelout=pan.Panel(items=paneltemp.items,major_axis=paneltemp.major_axis,
                                   minor_axis=paneltemp.minor_axis)
                for i in paneltemp.items:
                    dftemp=paneltemp[i]
                    if starttime is not None:
                        dftemp = dftemp.loc[dftemp.index>=starttime]                                                     
                    if endtime is not None:
                        dftemp = dftemp.loc[dftemp.index<=endtime]  
                    if timestep is not None:
                        dftemp = dftemp.resample(timestep, how ='ffill')
                    panelout[i]=dftemp
                templist.append(panelout)
            self.scenepanel=templist
        if verbose:
            print '... Done!'
        
        if self.sigma:
            self.calculate_sigma(winsize=sigma_winsize)
        
        if self.SNR:
            self.calculate_SNR(winsize=SNR_winsize)
                
        return self
        
    def range_cor(self):
        
        dataout = deepcopy(self.data)        
        bg = [self.header['bg_avg2'],self.header['bg_avg1']]
      
        for n in range(self.header['numchans'][0]): 
            rsq = (np.array(dataout[n].columns, dtype=float))**2
            for i in self.data[n].index:
                dataout[n].ix[i] = (self.data[n].ix[i] - bg[n].ix[i])*rsq
        
        self.rsq = dataout        
        return self
        
    def calculate_NRB(self, showplots = False,verbose=False):
        
        """
        Extracts data from MPL calibration files and applies it
        to a range-corrected set of mini-MPL data to convert from counts to attenuated backscatter
        """
        if sys.platform == 'win32':
            topdir = 'C:\Users\dashamstyr\Dropbox\Lidar Files\MPL Data\Calibration File Archive'
        else:
            topdir = '/data/lv1/pcottle/MPLCalibration'
        
        #if data were collected before June 2013, they were collected with MPL5008 and require
        #the associated calibration files
        version=np.int(self.header['version'][0])
        unitnum=np.int(self.header['unitnum'][0])
        tempdat = deepcopy(self.data)
        altvals = np.array(tempdat[0].columns, dtype='float')
        numchans=self.header['numchans'][0]
        if verbose:
            print "Calculating NRB"
        
        #first obtain filenames for deadtime, afterpulse, and overlap corrections
        if unitnum==5004:
            deadtimefile = os.path.join(topdir,'deadtimepoly_5004.bin')
            overlapfile = os.path.join(topdir,'MiniMPL5004_Horizontal_201301141700.bin')
            afterpulsefile = os.path.join(topdir,'5004_afterpusle_201402130600.bin')
        elif unitnum==5008:
            deadtimefile = os.path.join(topdir,'MMPL5008_deadtime.bin')
            overlapfile = os.path.join(topdir,'MMPL5008_overlap.bin')
            afterpulsefile = os.path.join(topdir,'MMPL5008_afterpulse.bin')
        elif unitnum==5012:
            deadtimefile = os.path.join(topdir,'MMPL5012_SPCM22625_deadtime7.bin')
            overlapfile = os.path.join(topdir,'MMPL5012_Overlap_201307310000.bin')
            afterpulsefile = os.path.join(topdir,'MMPL5012_Afterpulse_201308051500.mpl.bin')
        else:
            print "{0} is not a recognized Unit Number!".format(unitnum)
            return unitnum
        
        
        #depending on version, extract values for correction factors and interpolate to match altitudes in self.data
        if unitnum==5008:
            #extract deadtime correction factors
            with open(deadtimefile,'rb') as binfile:
                deadtimedat = array.array('f')
                while True:        
                    try:
                        deadtimedat.fromfile(binfile, 1)
                    except EOFError:
                        break
                
            coeffs = np.array(deadtimedat[::-1])
            
            deadtimecor = np.empty([numchans,len(tempdat[0].index),len(tempdat[0].columns)])
            for n in range(numchans):
                for i in range(len(tempdat[n].index)):
                    deadtimecor[n,i,:] = np.polynomial.polynomial.polyval(tempdat[n].iloc[i],coeffs)
            
            #extract afterpusle correction factors   
            with open(afterpulsefile, 'rb') as binfile:
                afterpulsedat = array.array('d')
                
                filedat = os.stat(afterpulsefile)
                numvals = filedat.st_size/8
                afterpulsedat.fromfile(binfile,numvals)
                
            numpairs = (numvals-1)/2
            mean_energy = np.array(afterpulsedat[0])
            aprange = np.array(afterpulsedat[1:numpairs+1])
            apvals_copol = np.array(afterpulsedat[numpairs+1:])/mean_energy
            apvals_crosspol=apvals_copol
            apvals=[apvals_copol,apvals_crosspol]
            
            #extract overlap correction factors
            with open(overlapfile, 'rb') as binfile:
                overlapdat = array.array('d')
                 
                filedat = os.stat(overlapfile)
                numvals = filedat.st_size/8
                overlapdat.fromfile(binfile,numvals)
            
            numpairs = numvals/2
            overrange = np.array(overlapdat[:numpairs])
            overvals = np.array(overlapdat[numpairs:])    
            
            altvals = np.array(tempdat[0].columns, dtype='float')
            interp_overlap = np.interp(altvals,overrange,overvals)
                
            for v in range(len(altvals)):
                if altvals[v] > max(overrange):
                    interp_overlap[v] = 1.0
    
        elif unitnum==5004 or unitnum==5012: 
            
            #extract deadtime correction factors
            with open(deadtimefile,'rb') as binfile:
                deadtimedat = array.array('f')
                while True:        
                    try:
                        deadtimedat.fromfile(binfile, 1)
                    except EOFError:
                        break
                
            coeffs = np.array(deadtimedat[::-1])
        
            deadtimecor = np.empty([numchans,len(tempdat[0].index),len(tempdat[0].columns)])
            for n in range(numchans):
                for i in range(len(tempdat[n].index)):
                    deadtimecor[n,i,:] = np.polynomial.polynomial.polyval(tempdat[n].iloc[i],coeffs)
                    
            #extract overlap correction factors
            with open(afterpulsefile, 'rb') as binfile:
                aprange = array.array('d')  #array of range values
                apvals_copol = array.array('d')  #array of correction values for copol
                apvals_crosspol = array.array('d')  #array of correction values for crosspol
        
                           
                temp = binfile.read(4)
                temp = binfile.read(2)
                apversion = struct.unpack('H',temp)[0] #version number for correction file
                temp = binfile.read(1)
                apnumchan = struct.unpack('B',temp)[0] #number of channels
                temp = binfile.read(4)
                apnumbins = struct.unpack('I',temp)[0] #number of correction bins
                temp = binfile.read(8)
                apenergy = struct.unpack('d',temp)[0] #laser energy during calibration in uJ
                
                if apversion != 3:
                    if verbose:
                        print 'WARNING!  This version of calibration file might be incompatible!'
                
                   
                if apnumchan == 2:
                    temp = binfile.read(8)
                    apbg_copol = struct.unpack('d',temp)[0]
                    temp = binfile.read(8)
                    apbg_crosspol = struct.unpack('d',temp)[0]
                    
                    aprange.fromfile(binfile,apnumbins)
                    apvals_copol.fromfile(binfile,apnumbins)
                    copol_norm = [n/apenergy for n in apvals_copol] #convert to counts/us
                    apvals_crosspol.fromfile(binfile,apnumbins)
                    crosspol_norm = [n/apenergy for n in apvals_crosspol] #convert to counts/us
                    apbg = [apbg_copol,apbg_crosspol]
                    apvals = [copol_norm,crosspol_norm]
                else:
                    if verbose:
                        print "Warning - wrong number of channels detected!"
            
            #extract overlap correction factors
            with open(overlapfile, 'rb') as binfile:
                overlapdat = array.array('d')
                 
                filedat = os.stat(overlapfile)
                numvals = filedat.st_size/8
                overlapdat.fromfile(binfile,numvals)
                
            numpairs = numvals/2
            overrange = np.array(overlapdat[:numpairs])
            overvals = np.array(overlapdat[numpairs:])    
        
            interp_overlap = np.interp(altvals,overrange,overvals)
                
            for v in range(len(altvals)):
                if altvals[v] > max(overrange):
                    interp_overlap[v] = 1.0  
        
        #now combine correction factors to calculate NRB from tempdat
        
        bg=(self.header['bg_avg2'],self.header['bg_avg1'])
        energy=self.header['energy']
        NRBout=[]
        
        for n in range(numchans):
            interp_afterpulse = np.interp(altvals,aprange,apvals[n])
            NRBtemp=pan.DataFrame(index=tempdat[n].index,columns=tempdat[n].columns)
            rsq = (np.array(tempdat[n].columns, dtype=float))**2
            for i in range(len(tempdat[n].index)):
                tempval = ((tempdat[n].iloc[i] - bg[n].iloc[i])*deadtimecor[n,i] - interp_afterpulse)/energy[i]                
                NRBtemp.iloc[i] = tempval*rsq/interp_overlap
            
            NRBout.append(NRBtemp.astype('float64'))
        
        self.NRB=NRBout
        
        if verbose:
            print "NRB calculation complete!"
        
        if showplots:
            temp = self.data[0].values
            maxval = temp.max()
            numsteps = maxval/100.
            x = np.arange(0,maxval,numsteps)
            y = np.polynomial.polynomial.polyval(x, coeffs)
            fig = plt.figure()
            fig.clf()
            ax = fig.add_subplot(111)
            ax.plot(x,y)
            ax.tick_params(axis='both', which='major', labelsize=20)
            ax.set_title('Deadtime Correction Curve',fontsize=30)
            ax.set_xlabel('Signal [mJ]',fontsize=25,labelpad=15)
            ax.set_ylabel('Correction Factor',fontsize=25,labelpad=15)

            fig = plt.figure()
            ax = fig.add_subplot(111)
            ax.plot(aprange[:10],apvals[1][:10],label='Copol Channel')
            ax.plot(aprange[:10],apvals[0][:10],label='Crosspol Channel')
            ax.legend()
            ax.tick_params(axis='both', which='major', labelsize=20)
            ax.set_title('Afterpulse Correction',fontsize=30)
            ax.set_xlabel('Altitude [m]', fontsize=25,labelpad=15)
            ax.set_ylabel('Correction Factor [mJ]', fontsize=25,labelpad=15)

            fig = plt.figure()
            ax = fig.add_subplot(111)
            ax.plot(overrange,overvals)
            ax.tick_params(axis='both', which='major', labelsize=20)
            ax.set_title('Overlap Correction',fontsize=30)
            ax.set_xlabel('Altitude [m]',fontsize=25,labelpad=15)
            ax.set_ylabel('Correction Factor',fontsize=25,labelpad=15)
            plt.show()                
        
        return self
    
    def calculate_depolrat(self,method='NRB',verbose=False):
        if verbose:
            print "Calculating Depolarization Ratio"
            
        if method=='NRB':
            copol = self.NRB[0]
            crosspol = self.NRB[1]
        elif method=='data':
            copol = self.data[0]
            crosspol = self.data[1]
        elif method=='rsq':
            copol = self.rsq[0]
            crosspol = self.rsq[1]
        
        depolMPL = (crosspol/copol).fillna(0.0)
        
        depolvals = (depolMPL/(depolMPL+1)).fillna(0.0)
        self.depolrat = [depolvals]
        
        if verbose:
            print "Depolarization ratio calculation complete!"
        
        return self
#    
    def calculate_sigma(self,winsize=10,verbose=False, datatypes=['all']):
        """
        Calculates stnadard deviations for mpl data
        
        inputs:
        num profs = number of vertical profiles to average together, defaults to 1
        datatypes = list of data types to callculate sigma for.  Could be 
                    'raw','rsq','NRB','depolrat', or 'all'
        
        output:
        self.sigma = a dict of pandas dataframes with datatype keys containing 
                    standard deviation values
        
        """
        
        if verbose:
            print "Calculating sigma"
        
        datasets=[]
        sigmadict={}
            
        for d in datatypes:
            if d=='data' or d=='all':
                datasets.append(('data',self.data))
            if d=='rsq' or d=='all':
                if self.rsq:
                    datasets.append(('rsq',self.rsq))
                elif verbose:
                    print "No RSQ data available for SNR calc"
            if d=='NRB' or d=='all':
                if self.NRB:
                    datasets.append(('NRB',self.NRB))
                elif verbose:
                    print "No NRB data available for SNR calc"
            if d=='depolrat' or d=='all':
                if self.depolrat:
                    datasets.append(('depolrat',self.depolrat))
                elif verbose:
                    print "No depolrat available for SNR calc"
#            if d=='backscatter' or d=='all':
#                if self.backscatter:
#                    datasets.append(('backscatter',self.backscatter))
#                elif verbose:
#                    print "No Backscatter available for SNR calc"
#            if d=='extinction' or d=='all':
#                if self.extinction:
#                    datasets.append(('extinction',self.extinction))
#                elif verbose:
#                    print "No extinction available for SNR calc"
        
        for dset_name,dset in datasets: 
            sigmadict[dset_name] = []
            if verbose:
                print "Calculating sigma values for {0}".format(dset_name)
            for n in range(len(dset)): 
                tempdat=dset[n]
                stdarray=pan.DataFrame(genfilt(tempdat,np.std,winsize),index=tempdat.index,
                                       columns=tempdat.columns)
                sigmadict[dset_name].append(stdarray)
                
        self.sigma=sigmadict

        if verbose:
            print "Sigma calculation done!"
            
        return self
    
    def calculate_SNR(self,bg_alt=None,winsize=10,verbose=False, datatypes=['all']):
        """
        Calculates signal to noise ratios for mpl data
        
        inputs:
        dfin = a pandas dataframe
        bg_alt = altitude above which signal is assumed to be purely background
                 if empty, topmost 100 data points are used
        num profs = number of vertical profiles to average together, defaults to 1
        datatypes = list of data types to callculate SNR for.  Could be 
                    'raw','rsq','NRB','depolrat', or 'all'
        
        output:
        self.SNR = a dict of pandas dataframes with datatype keys containing 
                    SNR values
        
        """
        bg = [self.header['bg_avg2'],self.header['bg_avg1']]
        
        if verbose:
            print "Calculating SNR"
        
        datasets=[]
        SNRdict={}
            
        for d in datatypes:
            if d=='data' or d=='all':
                datasets.append(('data',self.data))
            if d=='rsq' or d=='all':
                if self.rsq:
                    datasets.append(('rsq',self.rsq))
                elif verbose:
                    print "No RSQ data available for SNR calc"
            if d=='NRB' or d=='all':
                if self.NRB:
                    datasets.append(('NRB',self.NRB))
                elif verbose:
                    print "No NRB data available for SNR calc"
            if d=='depolrat' or d=='all':
                if self.depolrat:
                    datasets.append(('depolrat',self.depolrat))
                elif verbose:
                    print "No depolrat available for SNR calc"
#            if d=='backscatter' or d=='all':
#                if self.backscatter:
#                    datasets.append(('backscatter',self.backscatter))
#                elif verbose:
#                    print "No Backscatter available for SNR calc"
#            if d=='extinction' or d=='all':
#                if self.extinction:
#                    datasets.append(('extinction',self.extinction))
#                elif verbose:
#                    print "No extinction available for SNR calc"
        
        for dset_name,dset in datasets: 
            SNRdict[dset_name] = []
            for n in range(len(dset)): 
                tempdat=dset[n]
                if dset_name=='data':
                    if not bg_alt:
                        bg_alt=tempdat.columns[-10]
                    SNRtemp=pan.DataFrame(np.empty_like(tempdat.values),index=tempdat.index,columns=tempdat.columns)
                    for i in tempdat.index:
                        tempprof=tempdat.ix[i]                    
                        tempback=bg[n].ix[i]
                        tempvals=tempprof[bg_alt:].dropna()
                        sigmatemp=np.std(tempvals)
                        Ctemp=sigmatemp/np.mean(np.sqrt(np.abs(tempvals)))
                        SNR = lambda x: (x-tempback)/(Ctemp*np.sqrt(np.abs(x)))
                        SNRprof=np.array([SNR(v) for v in tempprof.values]).clip(0)
                        SNRtemp.ix[i]=SNRprof
                else:   
                    tempdat=dset[n].div((dset[n].columns.values)**2,axis=1)  #other datasets have been r-squared corrected, and this effect must be cancelled out
                    stdarray=pan.DataFrame(genfilt(tempdat,np.std,winsize),index=tempdat.index,
                                           columns=tempdat.columns)
                    meanarray=pan.DataFrame(genfilt(tempdat,np.mean,winsize),index=tempdat.index,
                                            columns=tempdat.columns)
                    SNRtemp=(meanarray/stdarray).fillna(0.0)
                SNRdict[dset_name].append(SNRtemp)
                
        self.SNR=SNRdict

        if verbose:
            print "SNR calculation done!"
            
        return self
    
    def calc_all(self,bg_alt=None,depol_method='NRB',winsize=10,showplots=False,verbose=False):
        """
        calculates all uncalculated fields for an mpl object
        """
        if self.rsq is None:
            self.range_cor()
        
        if self.NRB is None:
            self.calculate_NRB(showplots=showplots,verbose=verbose)
        
        if self.depolrat is None:
            self.calculate_depolrat(method=depol_method,verbose=verbose)
        
        if self.sigma is None:
            self.calculate_sigma(winsize=winsize,verbose=verbose)
            
        if self.SNR is None:
            self.calculate_SNR(bg_alt=bg_alt,winsize=winsize,verbose=verbose)
        
        return self
   
        
def set_dir(titlestring):
     
    # Make a top-level instance and hide since it is ugly and big.
    root = Tk()
    root.withdraw()
    
    # Make it almost invisible - no decorations, 0 size, top left corner.
    root.overrideredirect(True)
    root.geometry('0x0+0+0')
#    
    # Show window again and lift it to top so it can get focus,
    # otherwise dialogs will end up behind the terminal.
    root.deiconify()
    root.attributes("-topmost",1)
    root.focus_force()
    
    file_path = tkFileDialog.askdirectory(parent=root,title=titlestring)
     
    if file_path != "":
       return str(file_path)
     
    else:
       print "you didn't open anything!"
    
    # Get rid of the top-level instance once to make it actually invisible.
    root.destroy() 
          
def get_files(titlestring,filetype = ('.mpl','*.mpl')):
     
    # Make a top-level instance and hide since it is ugly and big.
    root = Tk()
    root.withdraw()
    
    # Make it almost invisible - no decorations, 0 size, top left corner.
    root.overrideredirect(True)
    root.geometry('0x0+0+0')
#    
    # Show window again and lift it to top so it can get focus,
    # otherwise dialogs will end up behind the terminal.
    root.deiconify()
    root.attributes("-topmost",1)
    root.focus_force()

    filenames = []
     
    filenames = tkFileDialog.askopenfilename(title=titlestring, filetypes=[filetype],multiple='True')
    
    #do nothing if already a python list
    if filenames == "": 
        print "You didn't open anything!"  
        return
    
    root.destroy()
    
    if isinstance(filenames,list):
        result = filenames   
    elif isinstance(filenames,tuple): 
        result = list(filenames)
    else:
        #http://docs.python.org/library/re.html
        #the re should match: {text and white space in brackets} AND anynonwhitespacetokens
        #*? is a non-greedy match for any character sequence
        #\S is non white space
    
        #split filenames string up into a proper python list
        result = re.findall("{.*?}|\S+",filenames)
    
        #remove any {} characters from the start and end of the file names
        result = [ re.sub("^{|}$","",i) for i in result ] 
    result.sort()
    return result

def resample_cols(dfin,newcols,verbose=False,method='interp'):
    oldcols=dfin.columns
        
    mincol=oldcols[0]
    maxcol=oldcols[-1]
    
    if mincol>newcols[0]:
        newcols=sorted([c for c in newcols if c>=mincol])
        if verbose:
            print "WARNING: Minimum column value reset to {0}".format(newcols[0])
        
    if maxcol<newcols[-1]:
        newcols=sorted([c for c in newcols if c<=maxcol])
        if verbose:
            print "WARNING: Maximum column value reset to {0}".format(newcols[-1])
    
    if len(newcols)==len(oldcols) and all(newcols==oldcols):
        return dfin
    else:
        newvalues=[]   
        for row in dfin.iterrows():
            if method=='interp':
                f=interp1d(oldcols,row[1].values)
                newvalues.append(f(newcols))
            elif method=='ffill':
                newrow=[]
                for col in newcols:
                    edgeval=row[1].groupby(row[1].index<=col).groups[True][-1]
                    newrow.append(row[1][edgeval])
                newvalues.append(newrow)    
            elif method=='bfill':
                newrow=[]
                for col in newcols:
                    edgeval=row[1].groupby(row[1].index<=col).groups[False][0]
                    newrow.append(row[1][edgeval])
                newvalues.append(newrow)            
        dfout=pan.DataFrame(data=newvalues,index=dfin.index,columns=newcols)       
        return dfout
        
    
def MPLtoHDF(filename, appendflag = 'False'):
    #not updated recently - do not use as is!
    #07/07/2015
    h5filename = filename.split('.')[0]+'_proc.h5'

    with open(filename,'rb') as binfile:
    
        profdat_copol = OrderedDict()
        profdat_crosspol = OrderedDict()
        header = OrderedDict()
        headerdat = OrderedDict()
        
        profnum = 0
        
        while True:
            try:
                intarray16 = array.array('H')
                intarray32 = array.array('I') # L is 8 byte on Xenon
                floatarray = array.array('f')
                byte_array = array.array('B')
                copolvals = array.array('f')  
                crosspolvals = array.array('f')  
                
                intarray16.fromfile(binfile, 8)
                intarray32.fromfile(binfile, 8)
                floatarray.fromfile(binfile, 2)
                intarray16.fromfile(binfile, 1)
                intarray32.fromfile(binfile, 1)
                floatarray.fromfile(binfile, 2)
                intarray16.fromfile(binfile, 3)
                floatarray.fromfile(binfile, 8)
                byte_array.fromfile(binfile, 2)
                floatarray.fromfile(binfile, 2)
                byte_array.fromfile(binfile, 1)
                intarray16.fromfile(binfile, 1)
                byte_array.fromfile(binfile, 1)
                intarray16.fromfile(binfile, 3)
                
                headerdat['unitnum'] = intarray16[0]
                headerdat['version'] = intarray16[1]
                year = intarray16[2]
                month = intarray16[3]
                day = intarray16[4]
                hour = intarray16[5]
                minute = intarray16[6]
                second = intarray16[7]
                dt = datetime.datetime(year,month,day,hour,minute,second)
            
                headerdat['shotsum'] = intarray32[0]  #total number of shots collected per profile
                headerdat['trigfreq'] = intarray32[1] #laser trigger frequency (usually 2500 Hz)
                headerdat['energy'] = intarray32[2]/1000.0  #mean of laser energy monitor in uJ
                headerdat['temp_0'] = intarray32[3]/1000.0  #mean of A/D#0 readings*100
                headerdat['temp_1'] = intarray32[4]/1000.0  #mean of A/D#1 readings*100
                headerdat['temp_2'] = intarray32[5]/1000.0  #mean of A/D#2 readings*100
                headerdat['temp_3'] = intarray32[6]/1000.0  #mean of A/D#3 readings*100
                headerdat['temp_4'] = intarray32[7]/1000.0  #mean of A/D#4 readings*100
                
                headerdat['bg_avg1'] = floatarray[0] #mean background signal value for channel 1
                headerdat['bg_std1'] = floatarray[1] #standard deviation of backgruond signal for channel 1
            
                headerdat['numchans'] = intarray16[8] #number of channels
                headerdat['numbins'] = intarray32[8] #total number of bins per channel
            
                headerdat['bintime'] = floatarray[2]  #bin width in seconds
                
                headerdat['rangecal'] = floatarray[3] #range offset in meters, default is 0
            
                headerdat['databins'] = intarray16[9]  #number of bins not including those used for background
                headerdat['scanflag'] = intarray16[10]  #0: no scanner, 1: scanner
                headerdat['backbins'] = intarray16[11]  #number of background bins
            
                headerdat['az'] = floatarray[4]  #scanner azimuth angle
                headerdat['el'] = floatarray[5]  #scanner elevation angle
                headerdat['deg'] = floatarray[6] #compass degrees (currently unused)
                headerdat['pvolt0'] = floatarray[7] #currently unused
                headerdat['pvolt1'] = floatarray[8] #currently unused
                headerdat['gpslat'] = floatarray[9] #GPS latitude in decimal degreees (-999.0 if no GPS)
                headerdat['gpslon'] = floatarray[10]#GPS longitude in decimal degrees (-999.0 if no GPS)
                headerdat['cloudbase'] = floatarray[11] #cloud base height in [m]
            
                headerdat['baddat'] = byte_array[0]  #0: good data, 1: bad data
                headerdat['version'] = byte_array[1] #version of file format.  current version is 1
            
                headerdat['bg_avg2'] = floatarray[12] #mean background signal for channel 2
                headerdat['bg_std2'] = floatarray[13] #mean background standard deviation for channel 2
            
                headerdat['mcs'] = byte_array[2]  #MCS mode register  Bit#7: 0-normal, 1-polarization
                                             #Bit#6-5: polarization toggling: 00-linear polarizer control
                                             #01-toggling pol control, 10-toggling pol control 11-circular pol control
            
                headerdat['firstbin'] = intarray16[12]  #bin # of first return data
                headerdat['systype'] = byte_array[3]   #0: standard MPL, 1: mini MPL
                headerdat['syncrate'] = intarray16[13]  #mini-MPL only, sync pulses seen per second
                headerdat['firstback'] = intarray16[14] #mini-MPL only, first bin used for background calcs
                headerdat['headersize2'] = intarray16[15] #size of additional header data (currently unused)
                
                if headerdat['headersize2'] > 128:
                    byte_array.fromfile(binfile, 1)
                    floatarray.fromfile(binfile, 6)
                    intarray16.fromfile(binfile, 1)
                    floatarray.fromfile(binfile, 2)
                
                    headerdat['Weatherstat'] = byte_array[4]   #0:weatehr station not used, 1:weather station used
                    headerdat['Int_temp'] = floatarray[14]  #Temperature inside in deg. Celsius
                    headerdat['Ext_temp'] = floatarray[15]  #Temp. outside
                    headerdat['Int_humid'] = floatarray[16] #Humidity inside in %
                    headerdat['Ext_humid'] = floatarray[17] #Hum. outside
                    headerdat['Dewpoint'] = floatarray[18]  #dewpoint in deg. Celsius
                    headerdat['Wnd_spd'] = floatarray[19]  #wind speed in km/h
                    headerdat['Wnd_dir'] = intarray16[16]  #wind direction in deg.
                    headerdat['Press'] = floatarray[20]  #Barometric pressure in hPa
                    headerdat['Rain'] = floatarray[21]  #Rain rate in mm/hr
                
                
                numbins = headerdat['numbins']
                numchans = headerdat['numchans'] 
                altstep = headerdat['bintime']*const.c #altitude step in meters
                maxalt = numbins*altstep
                minalt = headerdat['rangecal']
                altrange = np.arange(minalt,maxalt,altstep,dtype='float')
                
                if numchans == 2:
                    copolvals.fromfile(binfile, numbins) 
                    temp = np.array(copolvals)
                    profdat_crosspol[dt] = pan.Series(temp, index = altrange)
                    crosspolvals.fromfile(binfile, numbins) 
                    temp = np.array(crosspolvals)
                    profdat_copol[dt] = pan.Series(temp, index = altrange)
                else:
                    raise ValueError('Wrong number of channels')
                
                headerdat['profnum'] = profnum
                profnum += 1
                header[dt] = pan.Series(headerdat)
            except EOFError:
                break
        
        df_copol = pan.DataFrame.from_dict(profdat_copol, orient='index')
        df_copol.columns=altrange
        df_crosspol = pan.DataFrame.from_dict(profdat_crosspol, orient='index')
        df_crosspol.columns=altrange
        df_header = pan.DataFrame.from_dict(header, orient='index')
        
    store = pan.HDFStore(h5filename)
    store['copol_raw'] = df_copol
    store['crosspol_raw'] = df_crosspol
    store['header'] = df_header
    store.close()

def partdepolratcalc(depolin,beta_parallel,beta_mol,moldepolrat=0.0035):
    
    #default moldepolrat: narrow double filter allows only Cabannes line (see SPIE proc reference)
#    A = moldepolrat/(1+moldepolrat)
#    partdepolrat=(depolin*beta_parallel-A*beta_mol)/(beta_parallel-A*beta_mol)
    beta_p=beta_parallel-beta_mol    
    partdepolrat=(beta_parallel*depolin - beta_mol*moldepolrat)/beta_p
    
    return partdepolrat

def buffered_array(data,(x,y)):

    #create buffer around dataframe
    datashape = np.shape(data)
    
    b = int(np.ceil(y/2))
    t = y-b
    
    if len(datashape)==1:
        rows=datashape[0]
        newsize=(rows+y)
        newarray=np.empty(newsize)
        (newrows)=newarray.shape
        newarray[b:-t]=data
        newarray[:b]=data[:b]
        newarray[-t:]=data[-t:]
        
    else:
        rows=datashape[0]
        columns=datashape[1]    
        #simply copy first and last values to fill in buffers
        if x > 1:
            l = int(np.ceil(x/2))
            r = x-l
            newsize=(rows+x,columns+y)
            newarray=np.empty(newsize)
            (newrows,newcolums)=newarray.shape
            newarray[:l,b:-t]=data[0,:]
            newarray[l:-r,b:-t]=data
            newarray[-r:,b:-t]=data[-1,:]
        else:
            l=0
            r=0
            newsize=(rows,columns+y)
            newarray=np.empty(newsize)
            (newrows,newcolums)=newarray.shape
            newarray[:,b:-t]=data
    
        newarray[:,:b]=newarray[:,b:2*b]
        newarray[:,-t:]=newarray[:,-2*t:-t] 
    
    return newarray


    
def slopecalc(prof, winsize = 5):

    """
    Calculates slope of data for a single profile using a smoothing window of
    predetermined size
    
    inputs:
    prof:  a pandas series where index is altitude
    n:  number of consecutive values to use to calculate slope
    
    output:
    slope: output series,same size as input,with profile slopes
    """
    data = prof.values
    altrange = prof.index
    
    numvals = len(data)
    
    #calculate slopes of profile
      
    slopevals = np.empty_like(data)
    if winsize < 2:
        print "Sorry, need at least two values to calculate slope!"
        return
    elif winsize == 2:
        for i in range(numvals-1):
            slopevals[i] = (data[i+1]-data[i])/(altrange[i+1]-altrange[i])
        slopevals[-1] = slopevals[-2]
    else:
        #set up empty arrays for calculating slope with buffers on each end
        tempdat = np.empty(numvals+winsize)
        tempalt = np.empty_like(tempdat)
        #define buffer widths on left and right
        l = int(np.ceil(winsize/2))
        r = winsize-l
        #populate central segment with data and altrange
        tempdat[l:r] = data
        tempalt[l:r] = altrange
        #define window values for filling data and index buffers 
        lwinvals = data[:winsize]
        lwindex = altrange[:winsize]
        rwinvals = data[-winsize:]
        rwindex = altrange[-winsize:]
        #determine slope and intercept for left and right data windows      
        LA = np.array([lwindex,np.ones(winsize)])
        [lint,lslope]=np.linalg.lstsq(LA.T,lwinvals)
        RA = np.array([rwindex,np.ones(winsize)])
        [rint,rslope]=np.linalg.lstsq(RA.T,rwinvals)
        #calculate altitude index values for left and right buffers
        laltstep = altrange[1]-altrange[0]
        raltstep = altrange[-1]-altrange[-2]
        lold=altrange[0]
        rold=altrange[-1]
        for n in np.arange(l)[::-1]:
            tempalt[n]=lold-laltstep
            lold=tempalt[n]
        for m in np.arange(r):
            tempalt[m]=rold+raltstep
            rold=tempalt[m]
        #fill in data buffers
        tempdat[:l]=lint+lslope*tempalt[:l]
        tempdat[-r:]=rint+rslope*tempalt[-r:]
        
        #use linear regression to calculate slopes for sliding windows
        for v in range(numvals):
            windat = tempdat[v:v+winsize]
            winalt = tempalt[v:v+winsize]
            A = np.array([winalt,np.ones(winsize)])
            [inter,slopevals[v]]=np.linalg.lstsq(A.T,windat)
        
    slope = pan.Series(data=slopevals,index=altrange)
    return slope 

def SNR_mask_depol(mplin,**kwargs):
    
    SNRthreshold=kwargs.get('SNRthreshold',3)
    numprofs=kwargs.get('numprofs',1)
    bg_alt=kwargs.get('bg_alt',None)
    nopassval=kwargs.get('nopassval',float('nan'))
    inplace=kwargs.get('inplace',False)
    recalc=kwargs.get('recalc',False)
    datatype=kwargs.get('datatype','NRB')
    
    if recalc or not mplin.SNR:
        mplin = mplin.calculate_SNR(bg_alt,numprofs,datatypes=[datatype])
  
    #start by creating mask where areas that fall below SNRthreshold are zeroed out
    SNRmask = mplin.SNR[datatype][0]>=SNRthreshold
    
    if inplace:
        mplout=mplin
    else:
        mplout=deepcopy(mplin)
           
    if mplout.depolrat:
        mplout.depolrat[0]=mplin.depolrat[0]*SNRmask
        mplout.depolrat[0].replace(0,nopassval,inplace=True)
    return mplout 

def SNR_mask_scene(mplin,**kwargs):
    
    SNRthreshold=kwargs.get('SNRthreshold',3)
    numprofs=kwargs.get('numprofs',1)
    bg_alt=kwargs.get('bg_alt',None)
    inplace=kwargs.get('inplace',False)
    recalc=kwargs.get('recalc',False)
    datatype=kwargs.get('datatype','NRB')
    nopassval=kwargs.get('nopassval',None)
    
    if nopassval is None:
        nopassval={'Base':np.nan,
                   'Top':np.nan,
                   'Lidar_Ratio':0.0,
                   'Delta':0.0,
                   'Type':'Insufficient Signal',
                   'Sub-Type':'Insufficient Signal',
                   'Depol':0.0,
                   'colormask':9}
    
    if recalc or mplin.SNR is None:
        mplin = mplin.calculate_SNR(bg_alt,numprofs,datatypes=[datatype])
  
    #start by creating mask where areas that fall below SNRthreshold are zeroed out
    SNRmask = (mplin.SNR[datatype][0]>=SNRthreshold)#|(mplin.scenepanel[0]['Type']=='Clear Air')
#    SNRmask.replace(False,np.nan,inplace=True)
    
    if inplace:
        mplout=mplin
    else:
        mplout=deepcopy(mplin)
           
    if mplout.scenepanel:
        for sceneitem,maskval in nopassval.iteritems():
            mplout.scenepanel[0][sceneitem]=mplin.scenepanel[0][sceneitem]*SNRmask
            mplout.scenepanel[0][sceneitem].fillna(nopassval,inplace=True)
    return mplout 
    
def SNR_mask_all(mplin,**kwargs):
    
    SNRthreshold=kwargs.get('SNRthreshold',3)
    numprofs=kwargs.get('numprofs',1)
    bg_alt=kwargs.get('bg_alt',None)
    nopassval=kwargs.get('nopassval',float('nan'))
    inplace=kwargs.get('inplace',False)
    recalc=kwargs.get('recalc',False)
    datatype=kwargs.get('datatype','NRB')
    if recalc or not mplin.SNR:
        mplin = mplin.calculate_SNR(bg_alt,numprofs,datatype=['data'])
  
    #start by creating mask where areas that fall below SNRthreshold are zeroed out
    SNRmask = mplin.SNR[datatype][0]>=SNRthreshold
    
    if inplace:
        mplout=mplin
    else:
        mplout=deepcopy(mplin)
    
    for n in range(mplout.header['numchans'][0]):
        mplout.data[n]=mplin.data[n]*SNRmask
        mplout.data[n].replace(0,nopassval,inplace=True)
        if mplout.rsq:
            mplout.rsq[n]=mplin.rsq[n]*SNRmask
            mplout.rsq[n].replace(0,nopassval,inplace=True)        
        if mplout.NRB:
            mplout.NRB[n]=mplin.NRB[n]*SNRmask
            mplout.NRB[n].replace(0,nopassval,inplace=True)        
    if mplout.depolrat:
        mplout.depolrat[0]=mplin.depolrat[0]*SNRmask
        mplout.depolrat[0].replace(0,nopassval,inplace=True)
    return mplout  
    
def NRB_mask_create(dfin,**kwargs):
    
    """
    generates threshold altitudes to avoids spurious results by removing all 
    data beyond strong signal spikes
    
    """
    NRBthreshold=kwargs.get('NRBthreshold',3)
    NRBmin=kwargs.get('NRBmin',0.05)
    minalt=kwargs.get('minalt',0.150)
    numprofs=kwargs.get('numprofs',1)
    winsize=kwargs.get('winsize',5)
    
    #start by creating array of threshold altitudes and masking NRB copol
    #creates a new array with buffers to account for numprofs, winsize
    
    data = dfin.values
    altrange=dfin.columns.values
    (rows,columns) = data.shape
    minalt_index=np.where(altrange>=minalt)[0][0]
    newarray = buffered_array(data,(numprofs,winsize))
    (newrows,newcolums) = newarray.shape

    #set default values for cutoff to maximum altitude 
    threshalts=np.ones(len(dfin.index))*altrange[-1]
   
    for r in range(rows):
        tempprof=np.mean(newarray[r:r+numprofs],axis=0)
        for c in np.arange(minalt_index,columns):
            tempval = np.mean(tempprof[c:c+winsize])
            if tempval >= NRBthreshold:
                for c2 in np.arange(c,columns):
                    tempval = np.mean(tempprof[c2:c2+winsize])
                    if tempval <= NRBmin:
                        threshalts[r]=altrange[c2]
                        break 
    threshseries=pan.Series(data=threshalts,index=dfin.index)
    return threshseries

def NRB_mask_apply(dfin,threshseries,nopassval=np.nan,inplace=True):
    
    if inplace:
        dfout=dfin
    else:
        dfout=deepcopy(dfin)
    altvals = dfin.columns.values    
    for r in dfin.index:
        tempval=[x for x in altvals if x>=threshseries.ix[r]][0]
        dfout.ix[r,tempval:]=nopassval
    return dfout

def NRB_mask_all(MPLin,**kwargs):
    """
        uses a list of threshold altitudes, or generates one based on kwargs
        and applies it to all data sets within an MPL class object
    """
    
    threshseries=kwargs.get('threshseries',None)
    NRBmasktype=kwargs.get('NRBmasktype','profile')
    NRBthreshold=kwargs.get('NRBthreshold',3)
    NRBmin=kwargs.get('NRBmin',0.5)
    minalt=kwargs.get('minalt',0.150)
    numprofs=kwargs.get('numprofs',1)
    winsize=kwargs.get('winsize',3)
    nopassval=kwargs.get('nopassval',np.nan)
    inplace=kwargs.get('inplace',True)
    
    if inplace:
        MPLout=MPLin
    else:
        MPLout=deepcopy(MPLin)
    
    if threshseries is None:
        threshkwargs= {'NRBthreshold':NRBthreshold,'masktype':NRBmasktype,'NRBmin':NRBmin,'minalt':minalt,
                       'numprofs':numprofs,'winsize':winsize,'nopassval':nopassval}
        try:
            threshseries=NRB_mask_create(MPLout.NRB[0],**threshkwargs)
        except IndexError:
            print "NRB Mask can only work for MPL class object where NRB has been calculated!"
            return MPLout
    
    for n in range(MPLout.header['numchans'][0]):
        MPLout.data[n]=NRB_mask_apply(MPLout.data[n],threshseries)
        
        if MPLout.rsq is not None:
            MPLout.rsq[n]=NRB_mask_apply(MPLout.rsq[n],threshseries)
        if MPLout.NRB is not None:
            MPLout.NRB[n]=NRB_mask_apply(MPLout.NRB[n],threshseries)
        if MPLout.backscatter is not None:
            try:
                MPLout.backscatter[n]=NRB_mask_apply(MPLout.backscatter[n],threshseries)
            except IndexError:
                continue
        if MPLout.extinction is not None:
            try:
                MPLout.extinction[n]=NRB_mask_apply(MPLout.extinction[n],threshseries)
            except IndexError:
                continue
           
    if MPLout.depolrat is not None:
        MPLout.depolrat[0]=NRB_mask_apply(MPLout.depolrat[0],threshseries)
    
    if MPLout.scenepanel is not None:
        tempscene=MPLout.scenepanel[0]
        tempscene['Type']=NRB_mask_apply(tempscene['Type'],threshseries,nopassval='Insufficient Signal')
        tempscene['Sub-Type']=NRB_mask_apply(tempscene['Sub-Type'],threshseries,nopassval='Insufficient Signal')
        tempscene['colormask']=NRB_mask_apply(tempscene['colormask'],threshseries,nopassval=9)
        tempscene['Lidar_Ratio']=NRB_mask_apply(tempscene['Lidar_Ratio'],threshseries,nopassval=np.nan)
        tempscene['Depol']=NRB_mask_apply(tempscene['Depol'],threshseries,nopassval=np.nan)
        tempscene['Mean']=NRB_mask_apply(tempscene['Mean'],threshseries,nopassval=np.nan)
        tempscene['Base']=NRB_mask_apply(tempscene['Base'],threshseries,nopassval=np.nan)
        tempscene['Top']=NRB_mask_apply(tempscene['Top'],threshseries,nopassval=np.nan)
    return MPLout
    
#    def save_to_MPL(self,filename):
#        #currently not working!
#        import numpy as np
#        import array
#        
#        with open(filename, 'wb') as MPLout:
#        
#            intarray16 = array.array('H')
#            intarray32 = array.array('L')
#            floatarray = array.array('f')
#            byte_array = array.array('B')
#            datavals = array.array('f')
#                       
#            intarray16.append(self.header.unitnum)
#            intarray16.append(self.header.version)
#            d = self.header.datetime
#            
#            intarray16.append(d.year)
#            intarray16.append(d.month)
#            intarray16.append(d.day)
#            intarray16.append(d.hour)
#            intarray16.append(d.minute)
#            intarray16.append(d.second)
#    
#            intarray32.append(self.header.shotsum) #total number of shots collected per profile
#            intarray32.append(self.header.trigfreq) #laser trigger frequency (usually 2500 Hz)
#            intarray32.append(int(self.header.energy*1000))  #mean of laser energy monitor in uJ                      
#            intarray32.append(int(self.header.energy*1000))  #mean of A/D#0 readings*100
#            intarray32.append(int(self.header.temp_1*1000)) #mean of A/D#1 readings*100
#            intarray32.append(int(self.header.temp_2*1000))  #mean of A/D#2 readings*100
#            intarray32.append(int(self.header.temp_3*1000))  #mean of A/D#3 readings*100
#            intarray32.append(int(self.header.temp_4*1000))  #mean of A/D#4 readings*100
#            
#            floatarray.append(self.header.bg_avg1) #mean background signal value for channel 1
#            floatarray.append(self.header.bg_std1) #standard deviation of backgruond signal for channel 1
#    
#            intarray16.append(self.header.numchans) #number of channels
#            intarray32.append(self.header.numbins) #total number of bins per channel
#    
#            floatarray.append(self.header.bintime)  #bin width in seconds
#            
#            floatarray.append(self.header.rangecal) #range offset in meters, default is 0
#    
#            intarray16.append(self.header.databins) #number of bins not including those used for background
#            intarray16.append(self.header.scanflag)  #0: no scanner, 1: scanner
#            intarray16.append(self.header.backbins) #number of background bins
#    
#            floatarray.append(self.header.az)  #scanner azimuth angle
#            floatarray.append(self.header.el) #scanner elevation angle
#            floatarray.append(self.header.deg) #compass degrees (currently unused)
#            floatarray.append(self.header.pvolt0) #currently unused
#            floatarray.append(self.header.pvolt1) #currently unused
#            floatarray.append(self.header.gpslat) #GPS latitude in decimal degreees (-999.0 if no GPS)
#            floatarray.append(self.header.gpslon) #GPS longitude in decimal degrees (-999.0 if no GPS)
#            floatarray.append(self.header.cloudbase) #cloud base height in [m]
#    
#            byte_array.append(self.header.baddat)  #0: good data, 1: bad data
#            byte_array.append(self.header.version) #version of file format.  current version is 1
#    
#            floatarray.append(self.header.bg_avg2) #mean background signal for channel 2
#            floatarray.append(self.header.bg_std2) #mean background standard deviation for channel 2
#    
#            byte_array.append(self.header.mcs) #MCS mode register  Bit#7: 0-normal, 1-polarization
#                                         #Bit#6-5: polarization toggling: 00-linear polarizer control
#                                         #01-toggling pol control, 10-toggling pol control 11-circular pol control
#    
#            intarray16.append(self.header.firstbin)  #bin # of first return data
#            byte_array.append(self.header.systype)  #0: standard MPL, 1: mini MPL
#            intarray16.append(self.header.syncrate)  #mini-MPL only, sync pulses seen per second
#            intarray16.append(self.header.firstback) #mini-MPL only, first bin used for background calcs
#            intarray16.append(self.header.headersize2) #size of additional header data (currently unused)
#                     
#            copoldat = self.data[0].values
#            crosspoldat = self.data[1].values
#            
#            profile_buffer = np.zeros(32,dtype='float')
#            
#            for n in range(np.shape(copoldat)[0]):
#                datavals.fromlist(crosspoldat[n].tolist())
#                datavals.fromlist(copoldat[n].tolist())
#                datavals.fromlist(profile_buffer.tolist())
#        
#            temparray = intarray16[:8]                                    
#            temparray.tofile(MPLout)
#            temparray = intarray32[:8]                 
#            temparray.tofile(MPLout)
#            temparray = floatarray[:2]             
#            temparray.tofile(MPLout)
#            temparray = array.array('H',[intarray16[8]])             
#            temparray.tofile(MPLout)
#            temparray = array.array('L',[intarray32[8]])             
#            temparray.tofile(MPLout)
#            temparray = floatarray[2:4]            
#            temparray.tofile(MPLout)
#            temparray = intarray16[9:12]             
#            temparray.tofile(MPLout)
#            temparray = floatarray[4:12]
#            temparray.tofile(MPLout)
#            temparray = byte_array[:2]            
#            temparray.tofile(MPLout)
#            temparray = floatarray[12:14]            
#            temparray.tofile(MPLout)        
#            temparray = array.array('B',[byte_array[2]])             
#            temparray.tofile(MPLout)
#            temparray = array.array('H',[intarray16[12]])             
#            temparray.tofile(MPLout)
#            temparray = array.array('B',[byte_array[3]])             
#            temparray.tofile(MPLout)
#            temparray = intarray16[13:]             
#            temparray.tofile(MPLout)
#            datavals.tofile(MPLout)

if __name__ == '__main__':

    olddir = os.getcwd()
    
    os.chdir('C:\Users\dashamstyr\Dropbox\Lidar Files\MPL Data\DATA\Ucluelet Files')
    
    filename = get_files('Select raw file',filetype=('.mpl','*.mpl'))[0]
    
#    print 'Testing MPLtoHDF'
#    MPLtoHDF(filename)    
#    print 'done'
    
    print 'Testing MPL class functions'
    
    print 'Import MPL data from .mpl file'
    
    MPLtest = MPL()
    MPLtest.fromMPL(filename)
    
    print 'Done'
    
    print 'Calculate all corrections'
#    
    MPLtest.calc_all(verbose=True)
    
    print 'Calculate NRB mask'
    
    MPLmasked=NRB_mask_all(MPLtest,inplace=False)
#    os.chdir(olddir)