Mercurial > public > atmospheric_lidar / file comparison
comparison: licel.py
licel.py
- changeset 36
- a281a26f4626
- parent 35
- b1146c96f727
- child 37
- 7c76fdbdf1a3
equal
deleted
inserted
replaced
| 35:b1146c96f727 | 36:a281a26f4626 |
|---|---|
| 1 import numpy as np | |
| 2 import datetime | |
| 3 from generic import BaseLidarMeasurement, LidarChannel | |
| 4 import musa_netcdf_parameters | |
| 5 import musa_2009_netcdf_parameters | |
| 6 | |
| 7 licel_file_header_format = ['Filename', | |
| 8 'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', # Appart from Site that is read manually | |
| 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 | |
| 13 class LicelFile: | |
| 14 | |
| 15 def __init__(self, filename): | |
| 16 self.start_time = None | |
| 17 self.stop_time = None | |
| 18 self.import_file(filename) | |
| 19 self.calculate_physical() | |
| 20 self.filename = filename | |
| 21 | |
| 22 def calculate_physical(self): | |
| 23 for channel in self.channels.itervalues(): | |
| 24 channel.calculate_physical() | |
| 25 | |
| 26 def import_file(self, filename): | |
| 27 """Imports a licel file. | |
| 28 Input: filename | |
| 29 Output: object """ | |
| 30 | |
| 31 channels = {} | |
| 32 | |
| 33 with open(filename, 'rb') as f: | |
| 34 | |
| 35 self.read_header(f) | |
| 36 | |
| 37 # Check the complete header is read | |
| 38 a = f.readline() | |
| 39 | |
| 40 # Import the data | |
| 41 for current_channel_info in self.channel_info: | |
| 42 raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints'])) | |
| 43 a = np.fromfile(f, 'b', 1) | |
| 44 b = np.fromfile(f, 'b', 1) | |
| 45 | |
| 46 if (a[0] != 13) | (b[0] != 10): | |
| 47 print "Warning: No end of line found after record. File could be corrupt" | |
| 48 channel = LicelFileChannel(current_channel_info, raw_data, self.duration()) | |
| 49 | |
| 50 channel_name = channel.channel_name | |
| 51 | |
| 52 if channel_name in channels.keys(): | |
| 53 # If the analog/photon naming scheme is not enough, find a new one! | |
| 54 raise IOError('Trying to import two channels with the same name') | |
| 55 | |
| 56 channels[channel_name] = channel | |
| 57 | |
| 58 self.channels = channels | |
| 59 | |
| 60 | |
| 61 def read_header(self, f): | |
| 62 """ Read the header of a open file f. | |
| 63 | |
| 64 Returns raw_info and channel_info. Updates some object properties. """ | |
| 65 | |
| 66 #Read the first 3 lines of the header | |
| 67 raw_info = {} | |
| 68 channel_info = [] | |
| 69 | |
| 70 # Read first line | |
| 71 raw_info['Filename'] = f.readline().strip() | |
| 72 | |
| 73 # Read second line | |
| 74 second_line = f.readline() | |
| 75 | |
| 76 # Many licel files don't follow the licel standard. Specifically, the | |
| 77 # measurement site is not always 8 characters, and can include white | |
| 78 # spaces. For this, the site name is detect everything before the first | |
| 79 # date. For efficiency, the first date is found by the first '/'. | |
| 80 # e.g. assuming a string like 'Site name 01/01/2010 ...' | |
| 81 | |
| 82 site_name = second_line.split('/')[0][:-2] | |
| 83 clean_site_name = site_name.strip() | |
| 84 raw_info['Site'] = clean_site_name | |
| 85 raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1])) | |
| 86 | |
| 87 # Read third line | |
| 88 third_line = f.readline() | |
| 89 raw_info.update(match_lines(third_line, licel_file_header_format[2])) | |
| 90 | |
| 91 # Update the object properties based on the raw info | |
| 92 start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) | |
| 93 stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) | |
| 94 date_format = '%d/%m/%Y %H:%M:%S' | |
| 95 | |
| 96 self.start_time = datetime.datetime.strptime(start_string, date_format) | |
| 97 self.stop_time = datetime.datetime.strptime(stop_string, date_format) | |
| 98 self.latitude = float(raw_info['Latitude']) | |
| 99 self.longitude = float(raw_info['Longtitude']) | |
| 100 | |
| 101 # Read the rest of the header. | |
| 102 for c1 in range(int(raw_info['DataSets'])): | |
| 103 channel_info.append(match_lines(f.readline(), licel_file_channel_format)) | |
| 104 | |
| 105 self.raw_info = raw_info | |
| 106 self.channel_info = channel_info | |
| 107 | |
| 108 def duration(self): | |
| 109 """ Return the duration of the file. """ | |
| 110 dt = self.stop_time - self.start_time | |
| 111 return dt.seconds | |
| 112 | |
| 113 | |
| 114 class LicelFileChannel: | |
| 115 | |
| 116 def __init__(self, raw_info = None, raw_data = None, duration = None): | |
| 117 self.raw_info = raw_info | |
| 118 self.raw_data = raw_data | |
| 119 self.duration = duration | |
| 120 | |
| 121 @property | |
| 122 def wavelength(self): | |
| 123 if self.raw_info is not None: | |
| 124 wave_str = self.raw_info['Wavelength'] | |
| 125 wavelength = wave_str.split('.')[0] | |
| 126 return int(wavelength) | |
| 127 else: | |
| 128 return None | |
| 129 | |
| 130 @property | |
| 131 def channel_name(self): | |
| 132 ''' | |
| 133 Construct the channel name adding analog photon info to avoid duplicates | |
| 134 ''' | |
| 135 acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton']) | |
| 136 channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type) | |
| 137 return channel_name | |
| 138 | |
| 139 def analog_photon_string(self, analog_photon_number): | |
| 140 if analog_photon_number == '0': | |
| 141 string = 'an' | |
| 142 else: | |
| 143 string = 'ph' | |
| 144 return string | |
| 145 | |
| 146 def calculate_physical(self): | |
| 147 data = self.raw_data | |
| 148 | |
| 149 number_of_shots = float(self.raw_info['NShots']) | |
| 150 norm = data / number_of_shots | |
| 151 dz = float(self.raw_info['BinW']) | |
| 152 | |
| 153 if self.raw_info['AnalogPhoton']=='0': | |
| 154 # If the channel is in analog mode | |
| 155 ADCb = int(self.raw_info['ADCbits']) | |
| 156 ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV | |
| 157 channel_data = norm*ADCrange/((2**ADCb)-1) | |
| 158 | |
| 159 # print ADCb, ADCRange,cdata,norm | |
| 160 else: | |
| 161 # If the channel is in photoncounting mode | |
| 162 # Frequency deduced from range resolution! (is this ok?) | |
| 163 # c = 300 # The result will be in MHZ | |
| 164 # SR = c/(2*dz) # To account for pulse folding | |
| 165 # channel_data = norm*SR | |
| 166 # CHANGE: | |
| 167 # For the SCC the data are needed in photons | |
| 168 channel_data = norm *number_of_shots | |
| 169 #print res,c,cdata,norm | |
| 170 | |
| 171 #Calculate Z | |
| 172 number_of_bins = int(self.raw_info['DataPoints']) | |
| 173 self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) | |
| 174 self.dz = dz | |
| 175 self.number_of_bins = number_of_bins | |
| 176 self.data = channel_data | |
| 177 | |
| 178 | |
| 179 class LicelLidarMeasurement(BaseLidarMeasurement): | |
| 180 ''' | |
| 181 | |
| 182 ''' | |
| 183 extra_netcdf_parameters = musa_netcdf_parameters | |
| 184 raw_info = {} # Keep the raw info from the files | |
| 185 durations = {} # Keep the duration of the files | |
| 186 | |
| 187 def import_file(self, filename): | |
| 188 if filename in self.files: | |
| 189 print "File has been imported already:" + filename | |
| 190 else: | |
| 191 current_file = LicelFile(filename) | |
| 192 self.raw_info[current_file.filename] = current_file.raw_info | |
| 193 self.durations[current_file.filename] = current_file.duration() | |
| 194 | |
| 195 for channel_name, channel in current_file.channels.items(): | |
| 196 if channel_name not in self.channels: | |
| 197 self.channels[channel_name] = LicelChannel(channel) | |
| 198 self.channels[channel_name].data[current_file.start_time] = channel.data | |
| 199 self.files.append(current_file.filename) | |
| 200 | |
| 201 def append(self, other): | |
| 202 | |
| 203 self.start_times.extend(other.start_times) | |
| 204 self.stop_times.extend(other.stop_times) | |
| 205 | |
| 206 for channel_name, channel in self.channels.items(): | |
| 207 channel.append(other.channels[channel_name]) | |
| 208 | |
| 209 def _get_duration(self, raw_start_in_seconds): | |
| 210 ''' Return the duration for a given time scale. If only a single | |
| 211 file is imported, then this cannot be guessed from the time difference | |
| 212 and the raw_info of the file are checked. | |
| 213 ''' | |
| 214 | |
| 215 if len(raw_start_in_seconds) == 1: # If only one file imported | |
| 216 duration = self.durations.itervalues().next() # Get the first (and only) raw_info | |
| 217 duration_sec = duration | |
| 218 else: | |
| 219 duration_sec = np.diff(raw_start_in_seconds)[0] | |
| 220 | |
| 221 return duration_sec | |
| 222 | |
| 223 | |
| 224 class LicelChannel(LidarChannel): | |
| 225 | |
| 226 def __init__(self, channel_file): | |
| 227 c = 299792458.0 # Speed of light | |
| 228 self.wavelength = channel_file.wavelength | |
| 229 self.name = channel_file.channel_name | |
| 230 self.binwidth = channel_file.dz * 2 / c # in seconds | |
| 231 self.data = {} | |
| 232 self.resolution = channel_file.dz | |
| 233 self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z | |
| 234 self.points = channel_file.number_of_bins | |
| 235 self.rc = [] | |
| 236 self.duration = channel_file.duration | |
| 237 | |
| 238 def append(self, other): | |
| 239 if self.info != other.info: | |
| 240 raise ValueError('Channel info are different. Data can not be combined.') | |
| 241 | |
| 242 self.data = np.vstack([self.data, other.data]) | |
| 243 | |
| 244 def __unicode__(self): | |
| 245 return "<Licel channel: %s>" % self.name | |
| 246 | |
| 247 def __str__(self): | |
| 248 return unicode(self).encode('utf-8') | |
| 249 | |
| 250 | |
| 251 class Licel2009LidarMeasurement(LicelLidarMeasurement): | |
| 252 extra_netcdf_parameters = musa_2009_netcdf_parameters | |
| 253 | |
| 254 | |
| 255 def match_lines(f1,f2): | |
| 256 list1 = f1.split() | |
| 257 list2 = f2.split() | |
| 258 | |
| 259 if len(list1) != len(list2): | |
| 260 print "Warning: Combining lists of different lengths." | |
| 261 print "List 1: %s" % list1 | |
| 262 print "List 2: %s" % list2 | |
| 263 combined = zip(list2,list1) | |
| 264 combined = dict(combined) | |
| 265 return combined | |
| 266 |
