Initial files

Tue, 15 May 2012 16:43:26 +0200

author
ulalume3 <binietoglou@imaa.cnr.it>
date
Tue, 15 May 2012 16:43:26 +0200
changeset 0
9d2b98ecf23d
child 1
82b144ee09b2

Initial files

.hgignore file | annotate | diff | comparison | revisions
__init__.py file | annotate | diff | comparison | revisions
ciao.py file | annotate | diff | comparison | revisions
generic.py file | annotate | diff | comparison | revisions
licel.py file | annotate | diff | comparison | revisions
musa.py file | annotate | diff | comparison | revisions
musa_2009_netcdf_parameters.py file | annotate | diff | comparison | revisions
musa_netcdf_parameters.py file | annotate | diff | comparison | revisions
pearl.py file | annotate | diff | comparison | revisions
pearl_netcdf_parameters.py file | annotate | diff | comparison | revisions
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,3 @@
+syntax: glob
+*.pyc
+*.py~
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ciao.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,20 @@
+import milos
+
+class CiaoMixin:
+
+    def get_PT(self):
+        ''' Gets the pressure and temperature at station level from the Milos station.
+        The results are stored in the info dictionary.        
+        '''
+        
+        start_time = self.info['start_time']
+        stop_time = self.info['stop_time']
+        dt = stop_time - start_time
+        mean_time = start_time + dt/2
+        
+        # this guarantees that more that half the measurement period is taken into account
+        atm = milos.Atmospheric_condition(mean_time)
+        temperature = atm.get_mean('Air_Temperature', start_time, stop_time)
+        pressure = atm.get_mean('Air_Pressure', start_time, stop_time)
+        self.info['Temperature'] = temperature
+        self.info['Pressure'] = pressure
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/generic.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,543 @@
+# General imports
+import datetime
+from operator import itemgetter
+
+# Science imports
+import numpy as np
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+import netCDF4 as netcdf
+
+# CNR-IMAA specific imports
+import milos
+
+
+netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4'
+
+class BaseLidarMeasurement():
+    """ This is the general measurement object.
+    It is meant to become a general measurement object 
+    independent of the input files.
+    
+    Each subclass should implement the following:
+    * the import_file method.
+    * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters.
+    
+    You can override the get_PT method to define a custom procedure to get ground temperature and pressure.
+    The one implemented by default is by using the MILOS meteorological station data. 
+    
+    """
+    
+    def __init__(self, filelist= None):
+        self.info = {}
+        self.dimensions = {}
+        self.variables =  {}
+        self.channels = {}
+        self.attributes = {}
+        self.files = []
+        self.dark_measurement = None
+        if filelist:
+            self.import_files(filelist)
+        
+    def import_files(self,filelist):
+        for f in filelist:
+            self.import_file(f)
+        self.update()
+
+    def import_file(self,filename):
+        raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.')
+
+    def update(self):
+        '''
+        Update the the info, variables and dimensions of the lidar measurement based 
+        on the information found in the channels.
+        
+        Reading of the scan_angles parameter is not implemented.
+        '''
+        
+        # Initialize
+        start_time =[]
+        stop_time = []
+        points = []
+        all_time_scales = []
+        channel_name_list = []
+        
+        # Get the information from all the channels
+        for channel_name, channel in self.channels.items():
+            channel.update()
+            start_time.append(channel.start_time)
+            stop_time.append(channel.stop_time)
+            points.append(channel.points)
+            all_time_scales.append(channel.time)
+            channel_name_list.append(channel_name)
+            
+        # Find the unique time scales used in the channels
+        time_scales = set(all_time_scales)
+        
+        # Update the info dictionary
+        self.info['start_time'] = min(start_time)
+        self.info['stop_time'] = max(stop_time)
+        self.info['duration'] = self.info['stop_time'] - self.info['start_time']
+        
+        # Update the dimensions dictionary
+        self.dimensions['points'] = max(points)
+        self.dimensions['channels'] = len(self.channels)
+        # self.dimensions['scan angles'] = 1
+        self.dimensions['nb_of_time_scales'] = len(time_scales)
+        
+        # Update the variables dictionary
+        # Write time scales in seconds
+        raw_Data_Start_Time = []
+        raw_Data_Stop_Time = []
+        
+        for current_time_scale in list(time_scales): 
+            raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time
+            raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds
+            raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list
+            # Check if this time scale has measurements every 30 or 60 seconds.
+            
+            duration = self._get_duration(raw_start_in_seconds)
+            
+            raw_stop_in_seconds = raw_start_in_seconds + duration
+            raw_Data_Stop_Time.append(raw_stop_in_seconds)
+            
+        self.variables['Raw_Data_Start_Time']= raw_Data_Start_Time
+        self.variables['Raw_Data_Stop_Time']= raw_Data_Stop_Time
+    
+        # Make a dictionary to match time scales and channels
+        channel_timescales = []
+        for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales):
+            # The following lines are PEARL specific. The reason they are here is not clear.
+            # if channel_name =='1064BLR':
+            #     channel_name = '1064'
+            for (ts,n) in zip(time_scales, range(len(time_scales))):
+                if current_time_scale == ts:
+                    channel_timescales.append([channel_name,n])
+        self.variables['id_timescale'] = dict(channel_timescales)
+    
+    def _get_duration(self, raw_start_in_seconds):
+        ''' Return the duration for a given time scale. In some files (ex. Licel) this
+        can be specified from the files themselves. In others this must be guessed.
+         
+        '''
+        # The old method, kept here for reference
+        #dt = np.mean(np.diff(raw_start_in_seconds))
+        #for d in duration_list:
+        #    if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds
+        #        duration = d
+        
+        duration = np.diff(raw_start_in_seconds)[0]
+
+        return duration
+    
+    def subset_by_channels(self, channel_subset):
+        ''' Get a measurement object containing only the channels with names
+        contained in the channel_sublet list '''
+        
+        m = self.__class__() # Create an object of the same type as this one.
+        m.channels = dict([(channel, self.channels[channel]) for channel 
+                            in channel_subset])
+        m.update()
+        return m
+        
+    def subset_by_time(self, start_time, stop_time):
+    
+        if start_time > stop_time:
+            raise ValueError('Stop time should be after start time')
+            
+        if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']):
+            raise ValueError('The time interval specified is not part of the measurement')
+
+        m = self.__class__() # Create an object of the same type as this one.
+        for (channel_name, channel)  in self.channels.items():
+            m.channels[channel_name] = channel.subset_by_time(start_time, stop_time)
+        m.update()
+        return m
+    
+    def r_plot(self):
+        #Make a basic plot of the data.
+        #Should include some dictionary with params to make plot stable.
+        pass 
+    
+    def r_pdf(self):
+        # Create a pdf report using a basic plot and meta-data.
+        pass
+    
+    def save(self):
+        #Save the current state of the object to continue the analysis later.
+        pass
+        
+    def get_PT(self):
+        ''' Sets the pressure and temperature at station level .
+        The results are stored in the info dictionary.        
+        '''
+    
+        self.info['Temperature'] = 10.0
+        self.info['Pressure'] = 930.0
+        
+
+    def save_as_netcdf(self, filename):
+        """Saves the measurement in the netcdf format as required by the SCC.
+        Input: filename
+        """
+        params = self.extra_netcdf_parameters
+        needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure']
+        
+        for parameter in needed_parameters:
+            stored_value = self.info.get(parameter, None)
+            if stored_value is None:
+                raise ValueError('A value needs to be specified for %s' % parameter)
+        
+                
+        dimensions = {'points': 1,
+                'channels': 1,
+                'time': None,
+                'nb_of_time_scales': 1,
+                'scan_angles': 1,}  # Mandatory dimensions. Time bck not implemented
+
+        global_att = {'Measurement_ID': None,
+                    'RawData_Start_Date': None,
+                    'RawData_Start_Time_UT': None,
+                    'RawData_Stop_Time_UT': None,
+                    'RawBck_Start_Date': None,
+                    'RawBck_Start_Time_UT': None,
+                    'RawBck_Stop_Time_UT': None,
+                    'Sounding_File_Name': None,
+                    'LR_File_Name': None,
+                    'Overlap_File_Name': None,
+                    'Location': None,
+                    'System': None,
+                    'Latitude_degrees_north': None,
+                    'Longitude_degrees_east': None,
+                    'Altitude_meter_asl': None}
+
+        channel_variables = \
+                    {'channel_ID': (('channels', ), 'i'),
+                    'Background_Low': (('channels', ), 'd'),
+                    'Background_High': (('channels', ), 'd'),
+                    'LR_Input': (('channels', ), 'i'),
+                    'DAQ_Range': (('channels', ), 'd'), 
+                    'Depolarization_Factor': (('channels', ), 'd'), }
+                    
+        
+        channels = self.channels.keys()
+
+        input_values = dict(self.dimensions, **self.variables)
+        
+        # Add some mandatory global attributes
+        input_values['Measurement_ID'] = self.info['Measurement_ID']
+        input_values['RawData_Start_Date'] = '\'%s\'' % self.info['start_time'].strftime('%Y%m%d')
+        input_values['RawData_Start_Time_UT'] = '\'%s\'' % self.info['start_time'].strftime('%H%M%S')
+        input_values['RawData_Stop_Time_UT'] = '\'%s\'' % self.info['stop_time'].strftime('%H%M%S')
+        
+        # Add some optional global attributes
+        input_values['System'] = params.general_parameters['System']
+        input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north']
+        input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east']
+        input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl']
+        
+        # Open a netCDF4 file
+        f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file
+        
+        # Create the dimensions in the file
+        for (d,v) in dimensions.iteritems():
+            v = input_values.pop(d, v)
+            f.createDimension(d,v)
+        
+        # Create global attributes
+        for (attrib,value) in global_att.iteritems():
+            val = input_values.pop(attrib,value)
+            if val: 
+                exec('f.%s = %s' % (attrib,val))
+        
+        """ Variables """
+        # Write the values of fixes channel parameters
+        for (var,t) in channel_variables.iteritems():
+            temp_v = f.createVariable(var,t[1],t[0])
+            for (channel, n) in zip(channels, range(len(channels))):
+                temp_v[n] = params.channel_parameters[channel][var]
+        
+        # Write the id_timescale values
+        temp_id_timescale = f.createVariable('id_timescale','i',('channels',))
+        for (channel, n) in zip(channels, range(len(channels))):
+            temp_id_timescale[n] = self.variables['id_timescale'][channel]
+        
+        # Laser pointing angle
+        temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',))
+        temp_v[:] = params.general_parameters['Laser_Pointing_Angle']
+        
+        # Molecular calculation
+        temp_v = f.createVariable('Molecular_Calc','i')
+        temp_v[:] = params.general_parameters['Molecular_Calc']    
+
+        # Laser pointing angles of profiles
+        temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales'))
+        for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'],
+                                  range(len(self.variables['Raw_Data_Start_Time']))):
+            temp_v[:len(time_scale), n] = 0  # The lidar has only one laser pointing angle
+        
+        # Raw data start/stop time
+        temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales'))
+        temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales'))
+        for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'],
+                                             self.variables['Raw_Data_Stop_Time'],
+                                             range(len(self.variables['Raw_Data_Start_Time']))):
+            temp_raw_start[:len(start_time),n] = start_time
+            temp_raw_stop[:len(stop_time),n] = stop_time
+
+        #Laser shots
+        temp_v =  f.createVariable('Laser_Shots','i',('time','channels'))
+        for (channel,n) in zip(channels, range(len(channels))):
+            time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]])
+            temp_v[:time_length, n] = params.channel_parameters[channel]['Laser_Shots']
+        
+        #Raw lidar data
+        temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points'))
+        for (channel,n) in zip(channels, range(len(channels))):
+            c = self.channels[channel]
+            temp_v[:len(c.time),n, :c.points] = c.matrix
+        
+        self.add_dark_measurements_to_netcdf(f, channels)
+            
+        #Pressure at lidar station
+        temp_v = f.createVariable('Pressure_at_Lidar_Station','d')
+        temp_v[:] = self.info['Pressure']
+
+        #Temperature at lidar station
+        temp_v = f.createVariable('Temperature_at_Lidar_Station','d')
+        temp_v[:] = self.info['Temperature']   
+        
+        self.save_netcdf_extra(f)
+        f.close()
+   
+    def add_dark_measurements_to_netcdf(self, f, channels):
+
+        # Get dark measurements. If it is not given in self.dark_measurement
+        # try to get it using the get_dark_measurements method. If none is found
+        # return without adding something.
+        if self.dark_measurement is None:
+            self.dark_measurement = self.get_dark_measurements()
+        
+        if self.dark_measurement is None:
+            return
+        
+        dark_measurement = self.dark_measurement
+        
+        # Calculate the length of the time_bck dimensions
+        number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()]
+        max_number_of_profiles = np.max(number_of_profiles)
+        
+        # Create the dimension
+        f.createDimension('time_bck', max_number_of_profiles)
+         
+        # Save the dark measurement data
+        temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points'))
+        for (channel,n) in zip(channels, range(len(channels))):
+            c = dark_measurement.channels[channel]
+            temp_v[:len(c.time),n, :c.points] = c.matrix
+        
+        # Dark profile start/stop time
+        temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time','nb_of_time_scales'))
+        temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time','nb_of_time_scales'))
+        for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'],
+                                             dark_measurement.variables['Raw_Data_Stop_Time'],
+                                             range(len(dark_measurement.variables['Raw_Data_Start_Time']))):
+            temp_raw_start[:len(start_time),n] = start_time
+            temp_raw_stop[:len(stop_time),n] = stop_time
+        
+        # Dark measurement start/stop time
+        f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d')
+        f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S')
+        f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S')
+
+
+
+    def save_netcdf_extra(self, f):
+        pass
+        
+    def _gettime(self, date_str, time_str):        
+        t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S')
+        return t
+    
+    def plot(self):
+        for channel in self.channels:
+            self.channels[channel].plot(show_plot = False)
+        plt.show()
+    
+    def get_dark_measurements(self):
+        return None
+
+
+class Lidar_channel:
+
+    def __init__(self,channel_parameters):
+        c = 299792458 #Speed of light
+        self.wavelength = channel_parameters['name']
+        self.name = str(self.wavelength)
+        self.binwidth = float(channel_parameters['binwidth']) # in microseconds
+        self.data = {}
+        self.resolution = self.binwidth * c / 2
+        self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution/2.0 # Change: add half bin in the z 
+        self.points = len(channel_parameters['data'])
+        self.rc = []
+        self.duration  = 60
+        
+    def calculate_rc(self):
+        background = np.mean(self.matrix[:,4000:], axis = 1) #Calculate the background from 30000m and above
+        self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2)
+        
+    
+    def update(self):
+        self.start_time = min(self.data.keys())
+        self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration)
+        self.time = tuple(sorted(self.data.keys()))
+        sorted_data = sorted(self.data.iteritems(), key=itemgetter(0))
+        self.matrix =  np.array(map(itemgetter(1),sorted_data))
+    
+    def _nearest_dt(self,dtime):
+        margin = datetime.timedelta(seconds = 300)        
+        if ((dtime + margin) < self.start_time)| ((dtime - margin)   > self.stop_time):
+            print "Requested date not covered in this file"
+            raise
+        dt = abs(self.time - np.array(dtime))
+        dtmin = min(dt)
+
+        if dtmin > datetime.timedelta(seconds = 60):
+            print "Nearest profile more than 60 seconds away. dt = %s." % dtmin
+        ind_t = np.where(dt == dtmin)
+        ind_a= ind_t[0]
+        if len(ind_a) > 1:
+            ind_a = ind_a[0]
+        chosen_time = self.time[ind_a]
+        return chosen_time, ind_a
+    
+    def subset_by_time(self, start_time, stop_time):
+        
+        time_array = np.array(self.time)
+        condition = (time_array >= start_time) & (time_array <= stop_time)
+        
+        subset_time = time_array[condition]
+        subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time])
+        
+        #Create a list with the values needed by channel's __init__()
+        parameters_values = {'name': self.wavelength, 
+                             'binwidth': self.binwidth, 
+                             'data': subset_data[subset_time[0]],}
+
+        c = Lidar_channel(parameters_values)
+        c.data = subset_data
+        c.update()
+        return c
+        
+    
+    def profile(self,dtime, signal_type = 'rc'):
+        t, idx = self._nearest_dt(dtime)
+        if signal_type == 'rc':
+            data = self.rc
+        else:
+            data = self.matrix
+        
+        prof = data[idx,:][0]
+        return prof, t
+
+    def get_slice(self, starttime, endtime, signal_type = 'rc'):
+        if signal_type == 'rc':
+            data = self.rc
+        else:
+            data = self.matrix
+        tim = np.array(self.time)
+        starttime = self._nearest_dt(starttime)[0]
+        endtime = self._nearest_dt(endtime)[0]
+        condition = (tim >= starttime) & (tim <= endtime)
+        sl  = data[condition, :]
+        t = tim[condition]
+        return sl,t
+    
+    def av_profile(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'):
+        starttime = tim - duration/2
+        endtime = tim + duration/2
+        d,t = self.get_slice(starttime, endtime, signal_type = signal_type)
+        prof = np.mean(d, axis = 0)
+        tmin = min(t)
+        tmax = max(t)
+        tav = tmin + (tmax-tmin)/2
+        return prof,(tav, tmin,tmax)
+       
+    def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet):
+        #if filename is not None:
+        #    matplotlib.use('Agg')
+        
+        fig = plt.figure()
+        ax1 = fig.add_subplot(111)
+        self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom)
+        ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name)) 
+        
+        if filename is not None:
+            pass
+            #plt.savefig(filename)
+        else:
+            if show_plot:
+                plt.show()
+        #plt.close() ???
+    
+    def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', zoom = [0,12000,0,-1]):
+        
+        if signal_type == 'rc':
+            if len(self.rc) == 0:
+                self.calculate_rc()
+            data = self.rc
+        else:
+            data = self.matrix
+        
+        hmax_idx = self.index_at_height(zoom[1])
+        
+        ax1.set_ylabel('Altitude (km)')
+        ax1.set_xlabel('Time UTC')
+        #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064
+        #dateFormatter = mpl.dates.DateFormatter('%H.%M')
+        #hourlocator = mpl.dates.HourLocator()
+        
+        #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m')
+        #daylocator = mpl.dates.DayLocator()       
+        hourFormatter = mpl.dates.DateFormatter('%H.%M')
+        hourlocator = mpl.dates.AutoDateLocator(interval_multiples=True)
+
+        
+        #ax1.axes.xaxis.set_major_formatter(dayFormatter)
+        #ax1.axes.xaxis.set_major_locator(daylocator)
+        ax1.axes.xaxis.set_major_formatter(hourFormatter)
+        ax1.axes.xaxis.set_major_locator(hourlocator)
+
+        
+        ts1 = mpl.dates.date2num(self.start_time)
+        ts2 = mpl.dates.date2num(self.stop_time)
+        
+        
+        im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]],
+            aspect = 'auto',
+            origin = 'lower',
+            cmap = cmap,
+            #vmin = 0,
+            vmin = data[:,10:400].max() * 0.1,
+            #vmax = 1.4*10**7,
+            vmax = data[:,10:400].max() * 0.9,
+            extent = [ts1,ts2,self.z[zoom[0]]/1000.0, self.z[hmax_idx]/1000.0],
+            )
+        
+        cb1 = plt.colorbar(im1)
+        cb1.ax.set_ylabel('a.u.')
+
+    def index_at_height(self, height):
+        idx = np.array(np.abs(self.z - height).argmin())
+        if idx.size >1:
+            idx =idx[0]
+        return idx
+    
+def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID):
+    #Read the lidar files and select channels
+    temp_m = LidarClass(files)
+    m = temp_m.subset_by_channels(channels)
+    m.get_PT()
+    m.info['Measurement_ID'] = measurement_ID
+    m.save_as_netcdf(filename)
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/licel.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,207 @@
+import numpy as np
+import datetime
+from generic import BaseLidarMeasurement, Lidar_channel
+import musa_netcdf_parameters
+import musa_2009_netcdf_parameters
+
+licel_file_header_format = ['Filename',
+                            'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle',
+                            'LS1 Rate1 LS2 Rate2 DataSets', ]
+licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID'
+
+class LicelFile:
+    def __init__(self, filename):
+        self.start_time = None
+        self.stop_time = None
+        self.import_file(filename)
+        self.calculate_physical()
+        self.filename = filename
+    
+    def calculate_physical(self):
+        for channel in self.channels.itervalues():
+            channel.calculate_physical()
+      
+    def import_file(self, filename):
+        """Imports a licel file.
+        Input: filename
+        Output: object """
+
+        raw_info = {}
+        channels = {}
+        channel_info = []
+
+        f = open(filename, 'r')
+
+        #Read the first 3 lines of the header
+        raw_info = {}
+        for c1 in range(3):
+            raw_info.update(match_lines(f.readline(), licel_file_header_format[c1]))
+
+        start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime'])
+        stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime'])
+        date_format = '%d/%m/%Y %H:%M:%S'
+        self.start_time = datetime.datetime.strptime(start_string, date_format)
+        self.stop_time = datetime.datetime.strptime(stop_string, date_format)
+        self.latitude = float(raw_info['Latitude'])
+        self.longitude = float(raw_info['Longtitude'])
+        
+        # Read the rest of the header.
+        for c1 in range(int(raw_info['DataSets'])):
+            channel_info.append(match_lines(f.readline(), licel_file_channel_format))
+
+        # Check the complete header is read
+        a = f.readline()
+
+        # Import the data
+        for current_channel_info in channel_info:
+            raw_data = np.fromfile(f, 'l', int(current_channel_info['DataPoints']))
+            a = np.fromfile(f, 'b', 1)
+            b = np.fromfile(f, 'b', 1)
+            if (a[0] != 13) | (b[0] != 10):
+                print "Warning: No end of line found after record. File could be corrupt"
+            channel = LicelFileChannel(current_channel_info, raw_data)
+            
+            channel_name = channel.channel_name
+            if channel_name in channels.keys():
+                # If the analog/photon naming scheme is not enough, find a new one!
+                raise IOError('Trying to import two channels with the same name')
+            
+            channels[channel_name] = channel
+        f.close()
+
+        self.raw_info = raw_info
+        self.channels = channels
+    
+    
+
+class LicelFileChannel:
+    
+    def __init__(self, raw_info = None, raw_data = None):
+        self.raw_info = raw_info
+        self.raw_data = raw_data
+    
+    @property
+    def wavelength(self):
+        if self.raw_info is not None:
+            wave_str = self.raw_info['Wavelength']
+            wavelength = wave_str.split('.')[0]
+            return int(wavelength)
+        else:
+            return None
+            
+    @property
+    def channel_name(self):
+        '''
+        Construct the channel name adding analog photon info to avoid duplicates
+        '''
+        acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton'])
+        channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type)  
+        return channel_name
+    
+    def analog_photon_string(self, analog_photon_number):
+        if analog_photon_number == '0':
+            string = 'an'
+        else:
+            string = 'ph'
+        return string
+                
+    def calculate_physical(self):
+        data = self.raw_data
+
+        number_of_shots = float(self.raw_info['NShots'])
+        norm = data/number_of_shots
+        dz = float(self.raw_info['BinW'])
+ 
+        if self.raw_info['AnalogPhoton']=='0':
+            # If the channel is in analog mode
+            ADCb = int(self.raw_info['ADCbits'])
+            ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV
+            channel_data = norm*ADCrange/((2**ADCb)-1)
+   
+            # print ADCb, ADCRange,cdata,norm
+        else:
+            # If the channel is in photoncounting mode
+            # Frequency deduced from range resolution! (is this ok?)
+            # c = 300 # The result will be in MHZ
+            # SR  = c/(2*dz) # To account for pulse folding
+            # channel_data = norm*SR
+            # CHANGE:
+            # For the SCC the data are needed in photons
+            channel_data = norm*number_of_shots
+            #print res,c,cdata,norm
+
+        #Calculate Z
+        number_of_bins  = int(self.raw_info['DataPoints'])
+        self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)])
+        self.dz = dz
+        self.number_of_bins = number_of_bins
+        self.data = channel_data
+        
+class LicelLidarMeasurement(BaseLidarMeasurement):
+    '''
+    
+    '''
+    
+    extra_netcdf_parameters = musa_netcdf_parameters
+    
+    def import_file(self, filename):
+        if filename in self.files:
+            print "File has been imported already:" + filename
+        else:
+            current_file = LicelFile(filename)
+            for channel_name, channel in current_file.channels.items():
+                if channel_name not in self.channels:
+                    self.channels[channel_name] = Licel_channel(channel)
+                self.channels[channel_name].data[current_file.start_time] = channel.data   
+            self.files.append(current_file.filename)
+        
+    def append(self, other):
+
+        self.start_times.extend(other.start_times)
+        self.stop_times.extend(other.stop_times)
+
+        for channel_name, channel in self.channels.items():
+            channel.append(other.channels[channel_name])
+
+
+class Licel_channel(Lidar_channel):
+    
+    def __init__(self, channel_file):
+        c = 299792458.0 #Speed of light
+        self.wavelength = channel_file.wavelength
+        self.name = channel_file.channel_name
+        self.binwidth = channel_file.dz * 2 / c # in microseconds
+        self.data = {}
+        self.resolution = channel_file.dz
+        self.z = np.arange(channel_file.number_of_bins) *self.resolution + self.resolution/2.0 #Change: add half bin in the z 
+        self.points = channel_file.number_of_bins
+        self.rc = []
+        self.duration  = 60
+    
+    def append(self, other):
+        if self.info != other.info:
+            raise ValueError('Channel info are different. Data can not be combined.')
+
+        self.data = np.vstack([self.data, other.data])
+
+    def __unicode__(self):
+        return "<Licel channel: %s>" % self.info['Wavelength']
+    
+    def __str__(self):
+        return unicode(self).encode('utf-8')
+
+class Licel2009LidarMeasurement(LicelLidarMeasurement):
+    extra_netcdf_parameters = musa_2009_netcdf_parameters
+
+def match_lines(f1,f2):
+    list1 = f1.split()
+    list2 = f2.split()
+    
+    if len(list1) != len(list2):
+        print "Warning: Combining lists of different lengths."
+        print "List 1: %s" % list1
+        print "List 2: %s" % list2
+    combined = zip(list2,list1)
+    combined = dict(combined)
+    return combined
+    
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/musa.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,5 @@
+from licel import LicelLidarMeasurement
+from ciao import CiaoMixin
+
+class MusaLidarMeasurement(CiaoMixin, LicelLidarMeasurement):
+    pass
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/musa_2009_netcdf_parameters.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,97 @@
+general_parameters =  \
+{'System': '\'MUSA\'',
+ 'Laser_Pointing_Angle': 0,
+ 'Molecular_Calc': 0,
+ 'Latitude_degrees_north': 40.601039,
+ 'Longitude_degrees_east': 15.723771,
+ 'Altitude_meter_asl': 760}
+
+channel_parameters = \
+{'01064.o_an': {'channel_ID': 203,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':500.0,
+         'Depolarization_Factor': 0,},
+ '00355.o_an': {'channel_ID': 193,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0,},
+ '00355.o_ph': {'channel_ID': 194,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '00387.o_an': {'channel_ID': 195,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':20.0,
+         'Depolarization_Factor': 0,},
+ '00387.o_ph': {'channel_ID': 196,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '00532.p_an': {'channel_ID': 197,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0,},
+ '00532.p_ph': {'channel_ID': 198,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '00532.s_an': {'channel_ID': 199,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0.0441,},
+ '00532.s_ph': {'channel_ID': 200,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0.0441,},
+ '00607.o_an': {'channel_ID': 201,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':20.0,
+         'Depolarization_Factor': 0,},
+ '00607.o_ph': {'channel_ID': 202,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+         }
+
+#For testing. To be read from milos files.
+'''
+measurement_parameters = \
+{'Pressure_at_Lidar_Station': 930,
+ 'Temperature_at_Lidar_Station': 15,
+ 'Measurement_ID': '12345'}
+'''
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/musa_netcdf_parameters.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,97 @@
+general_parameters =  \
+{'System': '\'MUSA\'',
+ 'Laser_Pointing_Angle': 0,
+ 'Molecular_Calc': 0,
+ 'Latitude_degrees_north': 40.601039,
+ 'Longitude_degrees_east': 15.723771,
+ 'Altitude_meter_asl': 760}
+
+channel_parameters = \
+{'1064.o_an': {'channel_ID': 203,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':500.0,
+         'Depolarization_Factor': 0,},
+ '355.o_an': {'channel_ID': 193,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0,},
+ '355.o_ph': {'channel_ID': 194,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '387.o_an': {'channel_ID': 195,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':20.0,
+         'Depolarization_Factor': 0,},
+ '387.o_ph': {'channel_ID': 196,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532.p_an': {'channel_ID': 197,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0,},
+ '532.p_ph': {'channel_ID': 198,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532.s_an': {'channel_ID': 199,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':100.0,
+         'Depolarization_Factor': 0.0441,},
+ '532.s_ph': {'channel_ID': 200,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0.0441,},
+ '607.o_an': {'channel_ID': 201,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':20.0,
+         'Depolarization_Factor': 0,},
+ '607.o_ph': {'channel_ID': 202,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1200,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+         }
+
+#For testing. To be read from milos files.
+'''
+measurement_parameters = \
+{'Pressure_at_Lidar_Station': 930,
+ 'Temperature_at_Lidar_Station': 15,
+ 'Measurement_ID': '12345'}
+'''
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pearl.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,139 @@
+import datetime
+import os
+import glob
+
+import numpy as np
+
+from generic import BaseLidarMeasurement, Lidar_channel
+from ciao import CiaoMixin
+
+import pearl_netcdf_parameters
+from report_file import Report_file
+
+
+repository = '/mnt/storage/lidar_data/pearl/'
+
+
+class PearlLidarMeasurement(CiaoMixin, BaseLidarMeasurement):
+    
+    extra_netcdf_parameters = pearl_netcdf_parameters
+    
+    def import_file(self,filename):
+        ''' Import a pearl file. '''
+        
+        if filename in self.files:
+            print "File has been imported already:" + filename
+        else:
+            parameters, channels_dict = self.read_pearl_data(filename)
+            start_time = self._gettime(parameters['Acq_date'],parameters['Acq_start_time'])
+            
+            for channel_info in channels_dict.itervalues():
+                
+                if channel_info['name'] == '1064ALR':
+                    name = '1064'
+                    tm = start_time
+                elif channel_info['name'] == '1064BLR':
+                    name = '1064'
+                    tm = start_time + datetime.timedelta(seconds = 30)
+                else:
+                    name = channel_info['name']
+                    tm = start_time
+                if name not in self.channels:
+                    self.channels[name] = Lidar_channel(channel_info)
+                self.channels[name].data[tm] = channel_info['data']
+            self.files.append(filename)
+    
+    def read_pearl_data(self, filename):
+        '''
+        Reads a pearl file.
+        
+        Returns:
+        parameters - a dictionary of general parameters
+        channels   - a dictionary with keys the channel number and values lists
+                     [channel name, channel bin width, channel data].
+        '''
+        f = open(filename,'r') # Open the file
+        s = f.read(26) # Read the first 26 bytes
+        
+        #Get the values in a dictionary
+        parameters = {}
+        parameters['Acq_date'] = s[0:10] # First 10 bytes are the acquisition date.
+        parameters['Acq_start_time'] = s[10:20].strip() # Next 10 bytes are start time. Strip from trailing spaces.
+        parameters['Channel_no'] = np.fromstring(s[20:22], dtype = np.int16) # Next 2 bytes are the number of channels. Short integer.
+        parameters['Point_no'] = np.fromstring(s[22:26], dtype = np.int32) # Next 4 bytes are the number of points. Integer.
+        p = parameters # Just for less typing
+        
+        # Read the channel parameters
+        len = 20*p['Channel_no']    
+        s = f.read(len)
+        channels = {}
+        for (c1,n) in zip(range(0,len, 20),range(p['Channel_no'])):
+            channels[str(n)] = {'name' : s[c1+10:c1+20].strip(),
+                                'binwidth' : s[c1:c1+10].strip()}
+        
+        #Read the data
+        data = np.fromfile(f,dtype = np.float32)
+        #print filename + ': ' + str(data.size) +',' + str(p['Point_no']) +str(p['Channel_no'])
+        data = data.reshape(p['Point_no'],p['Channel_no'])
+        for ch in channels.iterkeys():
+            channels[ch]['data'] = data[:,int(ch)]
+        #Close the file
+        f.close()
+        return parameters,channels
+                
+
+def get_measurement_for_interval(start_time, stop_time):
+    ''' Searches for a pearl measurement based on a time interval     
+    '''
+    
+    correct_series = None
+    day = datetime.timedelta(hours = 24)
+    
+    if start_time > stop_time:
+            raise ValueError('Stop time should be after start time')
+    
+    
+    
+    #The list of directories based on the given time. Same, previous, Next day
+    possible_paths = [get_path(t) for t in [start_time - day, start_time, start_time + day] 
+                            if get_path(t) is not None]
+    for path in possible_paths:
+        try:
+            rf = Report_file(path)
+        except:
+            rf = None
+        
+        if rf is not None:
+            for serie in rf.series:
+                if (start_time > serie.starttime) and (stop_time < serie.endtime):
+                    correct_series = serie
+    
+    if correct_series:
+        files = correct_series.files.get('apd', []) + correct_series.files.get('mcb', [])
+        m_series = PearlLidarMeasurement(files)
+        m_subset = m_series.subset_by_time(start_time, stop_time)
+        return m_subset
+    else:
+        return None
+
+
+def get_channel(tim, channel = '1064'):
+    if channel =='1064':
+        extension = '*.apd'
+    else:
+        extension = '*.mcb'
+    
+    dirstr = get_path(tim)
+    
+    if not os.path.isdir(dirstr):
+        raise IOError('No measurement for that date (directory does not exist.).')
+        #No measurement for that date (directory does not exist.).
+    files = glob.glob(dirstr + extension)
+    m = PearlLidarMeasurement(files)
+    c = m.channels[channel]
+    return c
+    
+
+def get_path(tim):
+    dirstr = repository +  tim.strftime('%Y')+ '/' +tim.strftime('%d%m%Y') + '/'
+    return dirstr
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pearl_netcdf_parameters.py	Tue May 15 16:43:26 2012 +0200
@@ -0,0 +1,98 @@
+general_parameters =  \
+{'System': '\'PEARL\'',
+ 'Laser_Pointing_Angle': 0,
+ 'Molecular_Calc': 0,
+ 'Latitude_degrees_north': 40.601039,
+ 'Longitude_degrees_east': 15.723771,
+ 'Altitude_meter_asl': 760}
+
+channel_parameters = \
+{'355HR': {'channel_ID': 8,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '355LR': {'channel_ID': 14,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '386HR': {'channel_ID': 9,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '386LR': {'channel_ID': 15,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '407HR': {'channel_ID': 17,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532HR': {'channel_ID': 7,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532LR': {'channel_ID': 13,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532SHR': {'channel_ID': 12,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '532PHR': {'channel_ID': 11,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '607HR': {'channel_ID': 10,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 3000,
+         'LR_Input':1,
+         'DAQ_Range':0,
+         'Depolarization_Factor': 0,},
+ '1064': {'channel_ID': 6,
+         'Background_Low': 30000,
+         'Background_High': 50000,
+         'Laser_Shots': 1500,
+         'LR_Input':1,
+         'DAQ_Range':100,
+         'Depolarization_Factor': 0,},
+        
+         }
+
+#For testing. To be read from milos files.
+'''
+measurement_parameters = \
+{'Pressure_at_Lidar_Station': 930,
+ 'Temperature_at_Lidar_Station': 15,
+ 'Measurement_ID': '12345'}
+'''
+
+

mercurial