--- a/licel.py Mon Feb 22 16:50:03 2016 +0200
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,266 +0,0 @@
-import numpy as np
-import datetime
-from generic import BaseLidarMeasurement, LidarChannel
-import musa_netcdf_parameters
-import musa_2009_netcdf_parameters
-
-licel_file_header_format = ['Filename',
- 'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', # Appart from Site that is read manually
- '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 """
-
- channels = {}
-
- with open(filename, 'rb') as f:
-
- self.read_header(f)
-
- # Check the complete header is read
- a = f.readline()
-
- # Import the data
- for current_channel_info in self.channel_info:
- raw_data = np.fromfile(f, 'i4', 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, self.duration())
-
- 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
-
- self.channels = channels
-
-
- def read_header(self, f):
- """ Read the header of a open file f.
-
- Returns raw_info and channel_info. Updates some object properties. """
-
- #Read the first 3 lines of the header
- raw_info = {}
- channel_info = []
-
- # Read first line
- raw_info['Filename'] = f.readline().strip()
-
- # Read second line
- second_line = f.readline()
-
- # Many licel files don't follow the licel standard. Specifically, the
- # measurement site is not always 8 characters, and can include white
- # spaces. For this, the site name is detect everything before the first
- # date. For efficiency, the first date is found by the first '/'.
- # e.g. assuming a string like 'Site name 01/01/2010 ...'
-
- site_name = second_line.split('/')[0][:-2]
- clean_site_name = site_name.strip()
- raw_info['Site'] = clean_site_name
- raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1]))
-
- # Read third line
- third_line = f.readline()
- raw_info.update(match_lines(third_line, licel_file_header_format[2]))
-
- # Update the object properties based on the raw info
- 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))
-
- self.raw_info = raw_info
- self.channel_info = channel_info
-
- def duration(self):
- """ Return the duration of the file. """
- dt = self.stop_time - self.start_time
- return dt.seconds
-
-
-class LicelFileChannel:
-
- def __init__(self, raw_info = None, raw_data = None, duration = None):
- self.raw_info = raw_info
- self.raw_data = raw_data
- self.duration = duration
-
- @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
- raw_info = {} # Keep the raw info from the files
- durations = {} # Keep the duration of the files
-
- def import_file(self, filename):
- if filename in self.files:
- print "File has been imported already:" + filename
- else:
- current_file = LicelFile(filename)
- self.raw_info[current_file.filename] = current_file.raw_info
- self.durations[current_file.filename] = current_file.duration()
-
- for channel_name, channel in current_file.channels.items():
- if channel_name not in self.channels:
- self.channels[channel_name] = LicelChannel(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])
-
- def _get_duration(self, raw_start_in_seconds):
- ''' Return the duration for a given time scale. If only a single
- file is imported, then this cannot be guessed from the time difference
- and the raw_info of the file are checked.
- '''
-
- if len(raw_start_in_seconds) == 1: # If only one file imported
- duration = self.durations.itervalues().next() # Get the first (and only) raw_info
- duration_sec = duration
- else:
- duration_sec = np.diff(raw_start_in_seconds)[0]
-
- return duration_sec
-
-
-class LicelChannel(LidarChannel):
-
- 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 seconds
- 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 = channel_file.duration
-
- 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.name
-
- 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
- 
