binietoglou@0: import numpy as np binietoglou@0: import datetime binietoglou@0: from generic import BaseLidarMeasurement, Lidar_channel binietoglou@0: import musa_netcdf_parameters binietoglou@0: import musa_2009_netcdf_parameters binietoglou@0: binietoglou@0: licel_file_header_format = ['Filename', binietoglou@0: 'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', 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: binietoglou@0: class LicelFile: 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: raw_info = {} binietoglou@0: channels = {} binietoglou@0: channel_info = [] binietoglou@0: ulalume3@2: f = open(filename, 'rb') binietoglou@0: binietoglou@0: #Read the first 3 lines of the header binietoglou@0: raw_info = {} binietoglou@0: for c1 in range(3): binietoglou@0: raw_info.update(match_lines(f.readline(), licel_file_header_format[c1])) binietoglou@0: 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@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@0: 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@0: binietoglou@0: # Check the complete header is read binietoglou@0: a = f.readline() binietoglou@0: binietoglou@0: # Import the data binietoglou@0: for current_channel_info in channel_info: binietoglou@0: raw_data = np.fromfile(f, 'l', int(current_channel_info['DataPoints'])) binietoglou@0: a = np.fromfile(f, 'b', 1) binietoglou@0: b = np.fromfile(f, 'b', 1) ulalume3@2: binietoglou@0: if (a[0] != 13) | (b[0] != 10): binietoglou@0: print "Warning: No end of line found after record. File could be corrupt" binietoglou@0: channel = LicelFileChannel(current_channel_info, raw_data) binietoglou@0: binietoglou@0: channel_name = channel.channel_name binietoglou@0: if channel_name in channels.keys(): binietoglou@0: # If the analog/photon naming scheme is not enough, find a new one! binietoglou@0: raise IOError('Trying to import two channels with the same name') binietoglou@0: binietoglou@0: channels[channel_name] = channel binietoglou@0: f.close() binietoglou@0: binietoglou@0: self.raw_info = raw_info binietoglou@0: self.channels = channels binietoglou@0: binietoglou@0: binietoglou@0: binietoglou@0: class LicelFileChannel: binietoglou@0: binietoglou@0: def __init__(self, raw_info = None, raw_data = None): binietoglou@0: self.raw_info = raw_info binietoglou@0: self.raw_data = raw_data binietoglou@0: 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']) binietoglou@0: 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 binietoglou@0: 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 binietoglou@0: binietoglou@0: class LicelLidarMeasurement(BaseLidarMeasurement): binietoglou@0: ''' binietoglou@0: binietoglou@0: ''' binietoglou@0: binietoglou@0: extra_netcdf_parameters = musa_netcdf_parameters 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) binietoglou@0: for channel_name, channel in current_file.channels.items(): binietoglou@0: if channel_name not in self.channels: binietoglou@0: self.channels[channel_name] = Licel_channel(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: binietoglou@0: binietoglou@0: class Licel_channel(Lidar_channel): binietoglou@0: binietoglou@0: def __init__(self, channel_file): binietoglou@0: c = 299792458.0 #Speed of light binietoglou@0: self.wavelength = channel_file.wavelength binietoglou@0: self.name = channel_file.channel_name binietoglou@0: self.binwidth = channel_file.dz * 2 / c # in microseconds binietoglou@0: self.data = {} binietoglou@0: self.resolution = channel_file.dz binietoglou@0: 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 = [] binietoglou@0: self.duration = 60 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): binietoglou@0: return "" % self.info['Wavelength'] binietoglou@0: binietoglou@0: def __str__(self): binietoglou@0: return unicode(self).encode('utf-8') binietoglou@0: binietoglou@0: class Licel2009LidarMeasurement(LicelLidarMeasurement): binietoglou@0: extra_netcdf_parameters = musa_2009_netcdf_parameters binietoglou@0: 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: