Mercurial > public > atmospheric_lidar / file comparison
comparison: licel.py
licel.py
- changeset 0
- 9d2b98ecf23d
- child 2
- b11082bfcb4a
equal
deleted
inserted
replaced
| -1:000000000000 | 0:9d2b98ecf23d |
|---|---|
| 1 import numpy as np | |
| 2 import datetime | |
| 3 from generic import BaseLidarMeasurement, Lidar_channel | |
| 4 import musa_netcdf_parameters | |
| 5 import musa_2009_netcdf_parameters | |
| 6 | |
| 7 licel_file_header_format = ['Filename', | |
| 8 'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', | |
| 9 'LS1 Rate1 LS2 Rate2 DataSets', ] | |
| 10 licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID' | |
| 11 | |
| 12 class LicelFile: | |
| 13 def __init__(self, filename): | |
| 14 self.start_time = None | |
| 15 self.stop_time = None | |
| 16 self.import_file(filename) | |
| 17 self.calculate_physical() | |
| 18 self.filename = filename | |
| 19 | |
| 20 def calculate_physical(self): | |
| 21 for channel in self.channels.itervalues(): | |
| 22 channel.calculate_physical() | |
| 23 | |
| 24 def import_file(self, filename): | |
| 25 """Imports a licel file. | |
| 26 Input: filename | |
| 27 Output: object """ | |
| 28 | |
| 29 raw_info = {} | |
| 30 channels = {} | |
| 31 channel_info = [] | |
| 32 | |
| 33 f = open(filename, 'r') | |
| 34 | |
| 35 #Read the first 3 lines of the header | |
| 36 raw_info = {} | |
| 37 for c1 in range(3): | |
| 38 raw_info.update(match_lines(f.readline(), licel_file_header_format[c1])) | |
| 39 | |
| 40 start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) | |
| 41 stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) | |
| 42 date_format = '%d/%m/%Y %H:%M:%S' | |
| 43 self.start_time = datetime.datetime.strptime(start_string, date_format) | |
| 44 self.stop_time = datetime.datetime.strptime(stop_string, date_format) | |
| 45 self.latitude = float(raw_info['Latitude']) | |
| 46 self.longitude = float(raw_info['Longtitude']) | |
| 47 | |
| 48 # Read the rest of the header. | |
| 49 for c1 in range(int(raw_info['DataSets'])): | |
| 50 channel_info.append(match_lines(f.readline(), licel_file_channel_format)) | |
| 51 | |
| 52 # Check the complete header is read | |
| 53 a = f.readline() | |
| 54 | |
| 55 # Import the data | |
| 56 for current_channel_info in channel_info: | |
| 57 raw_data = np.fromfile(f, 'l', int(current_channel_info['DataPoints'])) | |
| 58 a = np.fromfile(f, 'b', 1) | |
| 59 b = np.fromfile(f, 'b', 1) | |
| 60 if (a[0] != 13) | (b[0] != 10): | |
| 61 print "Warning: No end of line found after record. File could be corrupt" | |
| 62 channel = LicelFileChannel(current_channel_info, raw_data) | |
| 63 | |
| 64 channel_name = channel.channel_name | |
| 65 if channel_name in channels.keys(): | |
| 66 # If the analog/photon naming scheme is not enough, find a new one! | |
| 67 raise IOError('Trying to import two channels with the same name') | |
| 68 | |
| 69 channels[channel_name] = channel | |
| 70 f.close() | |
| 71 | |
| 72 self.raw_info = raw_info | |
| 73 self.channels = channels | |
| 74 | |
| 75 | |
| 76 | |
| 77 class LicelFileChannel: | |
| 78 | |
| 79 def __init__(self, raw_info = None, raw_data = None): | |
| 80 self.raw_info = raw_info | |
| 81 self.raw_data = raw_data | |
| 82 | |
| 83 @property | |
| 84 def wavelength(self): | |
| 85 if self.raw_info is not None: | |
| 86 wave_str = self.raw_info['Wavelength'] | |
| 87 wavelength = wave_str.split('.')[0] | |
| 88 return int(wavelength) | |
| 89 else: | |
| 90 return None | |
| 91 | |
| 92 @property | |
| 93 def channel_name(self): | |
| 94 ''' | |
| 95 Construct the channel name adding analog photon info to avoid duplicates | |
| 96 ''' | |
| 97 acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton']) | |
| 98 channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type) | |
| 99 return channel_name | |
| 100 | |
| 101 def analog_photon_string(self, analog_photon_number): | |
| 102 if analog_photon_number == '0': | |
| 103 string = 'an' | |
| 104 else: | |
| 105 string = 'ph' | |
| 106 return string | |
| 107 | |
| 108 def calculate_physical(self): | |
| 109 data = self.raw_data | |
| 110 | |
| 111 number_of_shots = float(self.raw_info['NShots']) | |
| 112 norm = data/number_of_shots | |
| 113 dz = float(self.raw_info['BinW']) | |
| 114 | |
| 115 if self.raw_info['AnalogPhoton']=='0': | |
| 116 # If the channel is in analog mode | |
| 117 ADCb = int(self.raw_info['ADCbits']) | |
| 118 ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV | |
| 119 channel_data = norm*ADCrange/((2**ADCb)-1) | |
| 120 | |
| 121 # print ADCb, ADCRange,cdata,norm | |
| 122 else: | |
| 123 # If the channel is in photoncounting mode | |
| 124 # Frequency deduced from range resolution! (is this ok?) | |
| 125 # c = 300 # The result will be in MHZ | |
| 126 # SR = c/(2*dz) # To account for pulse folding | |
| 127 # channel_data = norm*SR | |
| 128 # CHANGE: | |
| 129 # For the SCC the data are needed in photons | |
| 130 channel_data = norm*number_of_shots | |
| 131 #print res,c,cdata,norm | |
| 132 | |
| 133 #Calculate Z | |
| 134 number_of_bins = int(self.raw_info['DataPoints']) | |
| 135 self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) | |
| 136 self.dz = dz | |
| 137 self.number_of_bins = number_of_bins | |
| 138 self.data = channel_data | |
| 139 | |
| 140 class LicelLidarMeasurement(BaseLidarMeasurement): | |
| 141 ''' | |
| 142 | |
| 143 ''' | |
| 144 | |
| 145 extra_netcdf_parameters = musa_netcdf_parameters | |
| 146 | |
| 147 def import_file(self, filename): | |
| 148 if filename in self.files: | |
| 149 print "File has been imported already:" + filename | |
| 150 else: | |
| 151 current_file = LicelFile(filename) | |
| 152 for channel_name, channel in current_file.channels.items(): | |
| 153 if channel_name not in self.channels: | |
| 154 self.channels[channel_name] = Licel_channel(channel) | |
| 155 self.channels[channel_name].data[current_file.start_time] = channel.data | |
| 156 self.files.append(current_file.filename) | |
| 157 | |
| 158 def append(self, other): | |
| 159 | |
| 160 self.start_times.extend(other.start_times) | |
| 161 self.stop_times.extend(other.stop_times) | |
| 162 | |
| 163 for channel_name, channel in self.channels.items(): | |
| 164 channel.append(other.channels[channel_name]) | |
| 165 | |
| 166 | |
| 167 class Licel_channel(Lidar_channel): | |
| 168 | |
| 169 def __init__(self, channel_file): | |
| 170 c = 299792458.0 #Speed of light | |
| 171 self.wavelength = channel_file.wavelength | |
| 172 self.name = channel_file.channel_name | |
| 173 self.binwidth = channel_file.dz * 2 / c # in microseconds | |
| 174 self.data = {} | |
| 175 self.resolution = channel_file.dz | |
| 176 self.z = np.arange(channel_file.number_of_bins) *self.resolution + self.resolution/2.0 #Change: add half bin in the z | |
| 177 self.points = channel_file.number_of_bins | |
| 178 self.rc = [] | |
| 179 self.duration = 60 | |
| 180 | |
| 181 def append(self, other): | |
| 182 if self.info != other.info: | |
| 183 raise ValueError('Channel info are different. Data can not be combined.') | |
| 184 | |
| 185 self.data = np.vstack([self.data, other.data]) | |
| 186 | |
| 187 def __unicode__(self): | |
| 188 return "<Licel channel: %s>" % self.info['Wavelength'] | |
| 189 | |
| 190 def __str__(self): | |
| 191 return unicode(self).encode('utf-8') | |
| 192 | |
| 193 class Licel2009LidarMeasurement(LicelLidarMeasurement): | |
| 194 extra_netcdf_parameters = musa_2009_netcdf_parameters | |
| 195 | |
| 196 def match_lines(f1,f2): | |
| 197 list1 = f1.split() | |
| 198 list2 = f2.split() | |
| 199 | |
| 200 if len(list1) != len(list2): | |
| 201 print "Warning: Combining lists of different lengths." | |
| 202 print "List 1: %s" % list1 | |
| 203 print "List 2: %s" % list2 | |
| 204 combined = zip(list2,list1) | |
| 205 combined = dict(combined) | |
| 206 return combined | |
| 207 |
