binietoglou@0: import numpy as np
binietoglou@0: import datetime
ulalume3@27: from generic import BaseLidarMeasurement, LidarChannel
binietoglou@0: import musa_netcdf_parameters
binietoglou@0: import musa_2009_netcdf_parameters
binietoglou@0: 
binietoglou@0: licel_file_header_format = ['Filename',
binietoglou@33:                             'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle',  # Appart from Site that is read manually
binietoglou@0:                             'LS1 Rate1 LS2 Rate2 DataSets', ]
binietoglou@0: licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID'
binietoglou@0: 
ioannis@22: 
binietoglou@0: class LicelFile:
ulalume3@27:     
binietoglou@0:     def __init__(self, filename):
binietoglou@0:         self.start_time = None
binietoglou@0:         self.stop_time = None
binietoglou@0:         self.import_file(filename)
binietoglou@0:         self.calculate_physical()
binietoglou@0:         self.filename = filename
binietoglou@0:     
binietoglou@0:     def calculate_physical(self):
binietoglou@0:         for channel in self.channels.itervalues():
binietoglou@0:             channel.calculate_physical()
binietoglou@0:       
binietoglou@0:     def import_file(self, filename):
binietoglou@0:         """Imports a licel file.
binietoglou@0:         Input: filename
binietoglou@0:         Output: object """
binietoglou@0: 
binietoglou@0:         channels = {}
binietoglou@33: 
binietoglou@33:         with open(filename, 'rb') as f:
binietoglou@33: 
binietoglou@33:             self.read_header(f)
binietoglou@33:             
binietoglou@33:             # Check the complete header is read
binietoglou@33:             a = f.readline()
binietoglou@0: 
binietoglou@33:             # Import the data
binietoglou@33:             for current_channel_info in self.channel_info:
binietoglou@33:                 raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints']))
binietoglou@33:                 a = np.fromfile(f, 'b', 1)
binietoglou@33:                 b = np.fromfile(f, 'b', 1)
binietoglou@33:                 
binietoglou@33:                 if (a[0] != 13) | (b[0] != 10):
binietoglou@33:                     print "Warning: No end of line found after record. File could be corrupt"
binietoglou@33:                 channel = LicelFileChannel(current_channel_info, raw_data, self.duration())
binietoglou@33:                 
binietoglou@33:                 channel_name = channel.channel_name
binietoglou@33:                 
binietoglou@33:                 if channel_name in channels.keys():
binietoglou@33:                     # If the analog/photon naming scheme is not enough, find a new one!
binietoglou@33:                     raise IOError('Trying to import two channels with the same name')
binietoglou@33:                 
binietoglou@33:                 channels[channel_name] = channel
binietoglou@33:                 
binietoglou@33:         self.channels = channels
binietoglou@33:     
binietoglou@33:     
binietoglou@33:     def read_header(self, f):
binietoglou@33:         """ Read the header of a open file f. 
binietoglou@33:         
binietoglou@33:         Returns raw_info and channel_info. Updates some object properties. """
binietoglou@33:         
binietoglou@0:         #Read the first 3 lines of the header
binietoglou@0:         raw_info = {}
binietoglou@33:         channel_info = []
binietoglou@33:         
binietoglou@33:         # Read first line
binietoglou@33:         raw_info['Filename'] = f.readline().strip()
binietoglou@33:         
binietoglou@33:         # Read second line
binietoglou@33:         second_line = f.readline()
binietoglou@33:         
binietoglou@33:         # Many licel files don't follow the licel standard. Specifically, the
binietoglou@33:         # measurement site is not always 8 characters, and can include white
binietoglou@33:         # spaces. For this, the site name is detect everything before the first 
binietoglou@33:         # date. For efficiency, the first date is found by the first '/'.
binietoglou@33:         # e.g. assuming a string like 'Site name 01/01/2010 ...' 
binietoglou@33:         
binietoglou@33:         site_name = second_line.split('/')[0][:-2]
binietoglou@33:         clean_site_name = site_name.strip()
binietoglou@33:         raw_info['Site'] = clean_site_name
binietoglou@33:         raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1]))
binietoglou@33:         
binietoglou@33:         # Read third line
binietoglou@33:         third_line = f.readline()
binietoglou@33:         raw_info.update(match_lines(third_line, licel_file_header_format[2]))
binietoglou@33:         
binietoglou@33:         # Update the object properties based on the raw info
binietoglou@0:         start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime'])
binietoglou@0:         stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime'])
binietoglou@0:         date_format = '%d/%m/%Y %H:%M:%S'
binietoglou@33:         
binietoglou@0:         self.start_time = datetime.datetime.strptime(start_string, date_format)
binietoglou@0:         self.stop_time = datetime.datetime.strptime(stop_string, date_format)
binietoglou@0:         self.latitude = float(raw_info['Latitude'])
binietoglou@0:         self.longitude = float(raw_info['Longtitude'])
binietoglou@33:             
binietoglou@0:         # Read the rest of the header.
binietoglou@0:         for c1 in range(int(raw_info['DataSets'])):
binietoglou@0:             channel_info.append(match_lines(f.readline(), licel_file_channel_format))
binietoglou@33:         
binietoglou@33:         self.raw_info = raw_info
binietoglou@33:         self.channel_info = channel_info
ioannis@21:             
ulalume3@27:     def duration(self):
ulalume3@27:         """ Return the duration of the file. """
ulalume3@27:         dt = self.stop_time - self.start_time
ulalume3@27:         return dt.seconds
ulalume3@27:         
ulalume3@27:         
binietoglou@0: class LicelFileChannel:
binietoglou@0:     
ulalume3@27:     def __init__(self, raw_info = None, raw_data = None, duration = None):
binietoglou@0:         self.raw_info = raw_info
binietoglou@0:         self.raw_data = raw_data
ulalume3@27:         self.duration = duration
ulalume3@27:         
binietoglou@0:     @property
binietoglou@0:     def wavelength(self):
binietoglou@0:         if self.raw_info is not None:
binietoglou@0:             wave_str = self.raw_info['Wavelength']
binietoglou@0:             wavelength = wave_str.split('.')[0]
binietoglou@0:             return int(wavelength)
binietoglou@0:         else:
binietoglou@0:             return None
binietoglou@0:             
binietoglou@0:     @property
binietoglou@0:     def channel_name(self):
binietoglou@0:         '''
binietoglou@0:         Construct the channel name adding analog photon info to avoid duplicates
binietoglou@0:         '''
binietoglou@0:         acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton'])
binietoglou@0:         channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type)  
binietoglou@0:         return channel_name
binietoglou@0:     
binietoglou@0:     def analog_photon_string(self, analog_photon_number):
binietoglou@0:         if analog_photon_number == '0':
binietoglou@0:             string = 'an'
binietoglou@0:         else:
binietoglou@0:             string = 'ph'
binietoglou@0:         return string
binietoglou@0:                 
binietoglou@0:     def calculate_physical(self):
binietoglou@0:         data = self.raw_data
binietoglou@0: 
binietoglou@0:         number_of_shots = float(self.raw_info['NShots'])
ulalume3@27:         norm = data / number_of_shots
binietoglou@0:         dz = float(self.raw_info['BinW'])
binietoglou@0:  
binietoglou@0:         if self.raw_info['AnalogPhoton']=='0':
binietoglou@0:             # If the channel is in analog mode
binietoglou@0:             ADCb = int(self.raw_info['ADCbits'])
binietoglou@0:             ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV
binietoglou@0:             channel_data = norm*ADCrange/((2**ADCb)-1)
binietoglou@0:    
binietoglou@0:             # print ADCb, ADCRange,cdata,norm
binietoglou@0:         else:
binietoglou@0:             # If the channel is in photoncounting mode
binietoglou@0:             # Frequency deduced from range resolution! (is this ok?)
binietoglou@0:             # c = 300 # The result will be in MHZ
binietoglou@0:             # SR  = c/(2*dz) # To account for pulse folding
binietoglou@0:             # channel_data = norm*SR
binietoglou@0:             # CHANGE:
binietoglou@0:             # For the SCC the data are needed in photons
ulalume3@27:             channel_data = norm *number_of_shots
binietoglou@0:             #print res,c,cdata,norm
binietoglou@0: 
binietoglou@0:         #Calculate Z
binietoglou@0:         number_of_bins  = int(self.raw_info['DataPoints'])
binietoglou@0:         self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)])
binietoglou@0:         self.dz = dz
binietoglou@0:         self.number_of_bins = number_of_bins
binietoglou@0:         self.data = channel_data
ioannis@22: 
binietoglou@0:         
binietoglou@0: class LicelLidarMeasurement(BaseLidarMeasurement):
binietoglou@0:     '''
binietoglou@0:     
binietoglou@0:     '''
binietoglou@0:     extra_netcdf_parameters = musa_netcdf_parameters
ulalume3@27:     raw_info = {}   # Keep the raw info from the files
ulalume3@27:     durations = {}  # Keep the duration of the files
binietoglou@0:     
binietoglou@0:     def import_file(self, filename):
binietoglou@0:         if filename in self.files:
binietoglou@0:             print "File has been imported already:" + filename
binietoglou@0:         else:
binietoglou@0:             current_file = LicelFile(filename)
ioannis@22:             self.raw_info[current_file.filename] = current_file.raw_info
ulalume3@27:             self.durations[current_file.filename] = current_file.duration()
ioannis@22:             
binietoglou@0:             for channel_name, channel in current_file.channels.items():
binietoglou@0:                 if channel_name not in self.channels:
ulalume3@27:                     self.channels[channel_name] = LicelChannel(channel)
binietoglou@0:                 self.channels[channel_name].data[current_file.start_time] = channel.data   
binietoglou@0:             self.files.append(current_file.filename)
binietoglou@0:         
binietoglou@0:     def append(self, other):
binietoglou@0: 
binietoglou@0:         self.start_times.extend(other.start_times)
binietoglou@0:         self.stop_times.extend(other.stop_times)
binietoglou@0: 
binietoglou@0:         for channel_name, channel in self.channels.items():
binietoglou@0:             channel.append(other.channels[channel_name])
binietoglou@0: 
ioannis@22:     def _get_duration(self, raw_start_in_seconds):
ioannis@22:         ''' Return the duration for a given time scale. If only a single
ioannis@22:         file is imported, then this cannot be guessed from the time difference
ioannis@22:         and the raw_info of the file are checked.
ioannis@22:         '''
ioannis@22:         
ioannis@22:         if len(raw_start_in_seconds) == 1: # If only one file imported
ulalume3@27:             duration = self.durations.itervalues().next() # Get the first (and only) raw_info
ioannis@31:             duration_sec = duration
ioannis@22:         else:
ioannis@22:             duration_sec = np.diff(raw_start_in_seconds)[0]
binietoglou@0: 
ioannis@22:         return duration_sec
ioannis@22: 
ioannis@22:        
ulalume3@27: class LicelChannel(LidarChannel):
binietoglou@0:     
binietoglou@0:     def __init__(self, channel_file):
binietoglou@33:         c = 299792458.0  # Speed of light
binietoglou@0:         self.wavelength = channel_file.wavelength
binietoglou@0:         self.name = channel_file.channel_name
ulalume3@27:         self.binwidth = channel_file.dz * 2 / c  # in seconds 
binietoglou@0:         self.data = {}
binietoglou@0:         self.resolution = channel_file.dz
ulalume3@27:         self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z 
binietoglou@0:         self.points = channel_file.number_of_bins
binietoglou@0:         self.rc = []
ulalume3@27:         self.duration  = channel_file.duration
binietoglou@0:     
binietoglou@0:     def append(self, other):
binietoglou@0:         if self.info != other.info:
binietoglou@0:             raise ValueError('Channel info are different. Data can not be combined.')
binietoglou@0: 
binietoglou@0:         self.data = np.vstack([self.data, other.data])
binietoglou@0: 
binietoglou@0:     def __unicode__(self):
ioannis@26:         return "<Licel channel: %s>" % self.name
binietoglou@0:     
binietoglou@0:     def __str__(self):
binietoglou@0:         return unicode(self).encode('utf-8')
binietoglou@0: 
ioannis@22: 
binietoglou@0: class Licel2009LidarMeasurement(LicelLidarMeasurement):
binietoglou@0:     extra_netcdf_parameters = musa_2009_netcdf_parameters
binietoglou@0: 
ioannis@22: 
binietoglou@0: def match_lines(f1,f2):
binietoglou@0:     list1 = f1.split()
binietoglou@0:     list2 = f2.split()
binietoglou@0:     
binietoglou@0:     if len(list1) != len(list2):
binietoglou@0:         print "Warning: Combining lists of different lengths."
binietoglou@0:         print "List 1: %s" % list1
binietoglou@0:         print "List 2: %s" % list2
binietoglou@0:     combined = zip(list2,list1)
binietoglou@0:     combined = dict(combined)
binietoglou@0:     return combined
binietoglou@0: