MPLread.py

def mplreader(filename):

    import numpy as np
    import array
    import datetime as dt

    #create a class instance to contain the data and header

    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,binfile):
            import numpy as np
            import array
            import datetime as dt
            import pandas as pan
            from scipy import constants as const
            
            class header:
                def __init__(self,intarray16,intarray32,floatarray,bytearray):
                    self.unitnum = intarray16[0]
                    self.version = intarray16[1]
                    year = intarray16[2]
                    month = intarray16[3]
                    day = intarray16[4]
                    hour = intarray16[5]
                    minute = intarray16[6]
                    second = intarray16[7]

                    self.datetime = dt.datetime(year,month,day,hour,minute,second)

                    self.shotsum = intarray32[0]  #total number of shots collected per profile
                    self.trigfreq = intarray32[1] #laser trigger frequency (usually 2500 Hz)
                    self.energy = intarray32[2]/1000.0  #mean of laser energy monitor in uJ
                    self.temp_0 = intarray32[3]/1000.0  #mean of A/D#0 readings*100
                    self.temp_1 = intarray32[4]/1000.0  #mean of A/D#1 readings*100
                    self.temp_2 = intarray32[5]/1000.0  #mean of A/D#2 readings*100
                    self.temp_3 = intarray32[6]/1000.0  #mean of A/D#3 readings*100
                    self.temp_4 = intarray32[7]/1000.0  #mean of A/D#4 readings*100
                    
                    self.bg_avg1 = floatarray[0] #mean background signal value for channel 1
                    self.bg_std1 = floatarray[1] #standard deviation of backgruond signal for channel 1

                    self.numchans = intarray16[8] #number of channels
                    self.numbins = intarray32[8] #total number of bins per channel

                    self.bintime = floatarray[2]  #bin width in seconds
                    
                    self.rangecal = floatarray[3] #range offset in meters, default is 0

                    self.databins = intarray16[9]  #number of bins not including those used for background
                    self.scanflag = intarray16[10]  #0: no scanner, 1: scanner
                    self.backbins = intarray16[11]  #number of background bins

                    self.az = floatarray[4]  #scanner azimuth angle
                    self.el = floatarray[5]  #scanner elevation angle
                    self.deg = floatarray[6] #compass degrees (currently unused)
                    self.pvolt0 = floatarray[7] #currently unused
                    self.pvolt1 = floatarray[8] #currently unused
                    self.gpslat = floatarray[9] #GPS latitude in decimal degreees (-999.0 if no GPS)
                    self.gpslon = floatarray[10]#GPS longitude in decimal degrees (-999.0 if no GPS)
                    self.cloudbase = floatarray[11] #cloud base height in [m]

                    self.baddat = bytearray[0]  #0: good data, 1: bad data
                    self.version = bytearray[1] #version of file format.  current version is 1

                    self.bg_avg2 = floatarray[12] #mean background signal for channel 2
                    self.bg_std2 = floatarray[13] #mean background standard deviation for channel 2

                    self.mcs = bytearray[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

                    self.firstbin = intarray16[12]  #bin # of first return data
                    self.systype = bytearray[3]   #0: standard MPL, 1: mini MPL
                    self.syncrate = intarray16[13]  #mini-MPL only, sync pulses seen per second
                    self.firstback = intarray16[14] #mini-MPL only, first bin used for background calcs
                    self.headersize2 = intarray16[15] #size of additional header data (currently unused)
            
            integer16 = array.array('H')
            integer32 = array.array('L')
            floats = array.array('f')
            bytes = array.array('B')
            datavals = array.array('f')
            
            
            integer16.fromfile(binfile, 8)
            integer32.fromfile(binfile, 8)
            floats.fromfile(binfile, 2)
            integer16.fromfile(binfile, 1)
            integer32.fromfile(binfile, 1)
            floats.fromfile(binfile, 2)
            integer16.fromfile(binfile, 3)
            floats.fromfile(binfile, 8)
            bytes.fromfile(binfile, 2)
            floats.fromfile(binfile, 2)
            bytes.fromfile(binfile, 1)
            integer16.fromfile(binfile, 1)
            bytes.fromfile(binfile, 1)
            integer16.fromfile(binfile, 3)
            
            print integer16
            print integer32
            print floats
            print bytes
            
            self.header = header(integer16,integer32,floats, bytes)
            
            
            numbins = self.header.numbins
            numchans = self.header.numchans          
            
            while True:
                try:
                    datavals.fromfile(binfile,(numbins*numchans))
                    
                except EOFError:
                    break
                
            altstep = self.header.bintime*const.c #altitude step in meters
            maxalt = numbins*altstep
            minalt = self.header.rangecal
            altrange = np.arange(minalt,maxalt,altstep,dtype='float')
            
            timestep = dt.timedelta(seconds = self.header.shotsum/self.header.trigfreq) #time between profiles in seconds
            dt = self.header.datetime
          
            if numchans == 2:
                
                profdat_copol = []
                profdat_crosspol = []
                
                n = 0
                itermax = 4*numbins
                while n < len(datavals):
                    profdat_crosspol.append(datavals[n:n+numbins]) 
                    n += numbins                               
                    profdat_copol.append(datavals[n:n+numbins])
                    
                    n += numbins + 32
                    
                    #note: this more adaptive algorithm should have worked
                    #I'm not sure why it doesn't so I'm going with the
                    #patch for now - 04-25-13
#                            m = n+numbins
#                            temp_prof_crosspol = datarray[n:m]
#                            temp_prof_copol = datarray[m:m+numbins]
#                            
#                            
#                            if all(temp_prof_crosspol > 1e-6):                             
#                                profdat_crosspol.append(temp_prof_crosspol)                                
#                                profdat_copol.append(temp_prof_copol)
#                                n += 2*numbins
#                            else:
#                                n += 1
                copoldat = np.array(profdat_copol)
                crosspoldat = np.array(profdat_crosspol)
                numsteps = np.shape(copoldat)[0]
    
                timerange = []
                for t in range(0,numsteps):
                    timerange.append(dt)
                    dt += timestep  
                    
                self.data = [pan.DataFrame(copoldat,index = timerange,columns = altrange),\
                pan.DataFrame(crosspoldat,index = timerange,columns = altrange)]
                
            else:
                raise ValueError('Wrong number of channels')

    
    binfile = open(filename,'rb')
    
    MPL_out = MPL(binfile)
    
    return MPL_out


if __name__=='__main__':

    import os, sys
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.colors import Normalize
    from matplotlib import cm

    olddir = os.getcwd()

    os.chdir('C:\Program Files (x86)\SigmaMPL\DATA')

    filename = '201304250900.mpl'

    MPLdat = mplreader(filename)
    
    copol = MPLdat.data[0]
    crosspol = MPLdat.data[1]
    
    print 'Copol Shape is'
    print np.shape(copol)
    print 'Crosspol Shape is'
    print np.shape(crosspol)
    
    cmap = cm.jet
    cmap.set_over('r')
    cmap.set_under('k')
    vmin= 0
    vmax= 1
    the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
    
    
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    im1 = ax1.pcolormesh(copol.values.T,cmap=cmap,norm=the_norm)
    cb = fig.colorbar(im1,extend='both')

    ax2 = fig.add_subplot(122)
    im2 = ax2.pcolormesh(crosspol.values.T,cmap=cmap,norm=the_norm)
    cb = fig.colorbar(im1,extend='both')
    
    plt.show()

    os.chdir(olddir)