Tue, 12 Dec 2017 17:24:04 +0200
Licel
ioannis@55 | 1 | import datetime |
ioannis@55 | 2 | import logging |
ioannis@55 | 3 | |
binietoglou@0 | 4 | import numpy as np |
i@101 | 5 | import pytz |
ioannis@55 | 6 | |
ulalume3@27 | 7 | from generic import BaseLidarMeasurement, LidarChannel |
binietoglou@0 | 8 | |
i@101 | 9 | logger = logging.getLogger(__name__) |
i@101 | 10 | |
binietoglou@0 | 11 | licel_file_header_format = ['Filename', |
ulalume3@47 | 12 | 'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', |
ulalume3@47 | 13 | # Appart from Site that is read manually |
binietoglou@0 | 14 | 'LS1 Rate1 LS2 Rate2 DataSets', ] |
binietoglou@0 | 15 | licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID' |
binietoglou@0 | 16 | |
ioannis@22 | 17 | |
binietoglou@0 | 18 | class LicelFile: |
i@101 | 19 | def __init__(self, filename, use_id_as_name=False, licel_timezone="UTC"): |
ulalume3@47 | 20 | self.filename = filename |
ioannis@55 | 21 | self.use_id_as_name = use_id_as_name |
binietoglou@0 | 22 | self.start_time = None |
binietoglou@0 | 23 | self.stop_time = None |
i@102 | 24 | self.licel_timezone = licel_timezone |
binietoglou@0 | 25 | self.import_file(filename) |
binietoglou@0 | 26 | self.calculate_physical() |
i@102 | 27 | |
binietoglou@0 | 28 | def calculate_physical(self): |
binietoglou@0 | 29 | for channel in self.channels.itervalues(): |
binietoglou@0 | 30 | channel.calculate_physical() |
ulalume3@47 | 31 | |
binietoglou@0 | 32 | def import_file(self, filename): |
binietoglou@0 | 33 | """Imports a licel file. |
binietoglou@0 | 34 | Input: filename |
binietoglou@0 | 35 | Output: object """ |
binietoglou@0 | 36 | |
binietoglou@0 | 37 | channels = {} |
binietoglou@33 | 38 | |
binietoglou@33 | 39 | with open(filename, 'rb') as f: |
binietoglou@33 | 40 | |
binietoglou@33 | 41 | self.read_header(f) |
ulalume3@47 | 42 | |
binietoglou@33 | 43 | # Check the complete header is read |
i@101 | 44 | f.readline() |
binietoglou@0 | 45 | |
binietoglou@33 | 46 | # Import the data |
binietoglou@33 | 47 | for current_channel_info in self.channel_info: |
binietoglou@33 | 48 | raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints'])) |
binietoglou@33 | 49 | a = np.fromfile(f, 'b', 1) |
binietoglou@33 | 50 | b = np.fromfile(f, 'b', 1) |
ulalume3@47 | 51 | |
binietoglou@33 | 52 | if (a[0] != 13) | (b[0] != 10): |
ioannis@55 | 53 | logging.warning("No end of line found after record. File could be corrupt: %s" % filename) |
ioannis@55 | 54 | channel = LicelFileChannel(current_channel_info, raw_data, self.duration(), |
ioannis@55 | 55 | use_id_as_name=self.use_id_as_name) |
ulalume3@47 | 56 | |
binietoglou@33 | 57 | channel_name = channel.channel_name |
ulalume3@47 | 58 | |
binietoglou@33 | 59 | if channel_name in channels.keys(): |
binietoglou@33 | 60 | # If the analog/photon naming scheme is not enough, find a new one! |
binietoglou@33 | 61 | raise IOError('Trying to import two channels with the same name') |
ulalume3@47 | 62 | |
binietoglou@33 | 63 | channels[channel_name] = channel |
ulalume3@47 | 64 | |
binietoglou@33 | 65 | self.channels = channels |
ulalume3@47 | 66 | |
binietoglou@33 | 67 | def read_header(self, f): |
binietoglou@33 | 68 | """ Read the header of a open file f. |
binietoglou@33 | 69 | |
binietoglou@33 | 70 | Returns raw_info and channel_info. Updates some object properties. """ |
ulalume3@47 | 71 | |
ulalume3@47 | 72 | # Read the first 3 lines of the header |
binietoglou@0 | 73 | raw_info = {} |
binietoglou@33 | 74 | channel_info = [] |
ulalume3@47 | 75 | |
binietoglou@33 | 76 | # Read first line |
binietoglou@33 | 77 | raw_info['Filename'] = f.readline().strip() |
ulalume3@47 | 78 | |
binietoglou@33 | 79 | # Read second line |
binietoglou@33 | 80 | second_line = f.readline() |
ulalume3@47 | 81 | |
binietoglou@33 | 82 | # Many licel files don't follow the licel standard. Specifically, the |
binietoglou@33 | 83 | # measurement site is not always 8 characters, and can include white |
binietoglou@33 | 84 | # spaces. For this, the site name is detect everything before the first |
binietoglou@33 | 85 | # date. For efficiency, the first date is found by the first '/'. |
binietoglou@33 | 86 | # e.g. assuming a string like 'Site name 01/01/2010 ...' |
ulalume3@47 | 87 | |
binietoglou@33 | 88 | site_name = second_line.split('/')[0][:-2] |
binietoglou@33 | 89 | clean_site_name = site_name.strip() |
binietoglou@33 | 90 | raw_info['Site'] = clean_site_name |
binietoglou@33 | 91 | raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1])) |
ulalume3@47 | 92 | |
binietoglou@33 | 93 | # Read third line |
binietoglou@33 | 94 | third_line = f.readline() |
binietoglou@33 | 95 | raw_info.update(match_lines(third_line, licel_file_header_format[2])) |
ulalume3@47 | 96 | |
binietoglou@33 | 97 | # Update the object properties based on the raw info |
binietoglou@0 | 98 | start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) |
binietoglou@0 | 99 | stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) |
binietoglou@0 | 100 | date_format = '%d/%m/%Y %H:%M:%S' |
ulalume3@47 | 101 | |
i@101 | 102 | try: |
i@101 | 103 | logger.debug('Creating timezone object %s' % self.licel_timezone) |
i@101 | 104 | timezone = pytz.timezone(self.licel_timezone) |
i@101 | 105 | except: |
i@101 | 106 | raise ValueError("Cloud not create time zone object %s" % self.licel_timezone) |
i@101 | 107 | |
i@101 | 108 | # According to pytz docs, timezones do not work with default datetime constructor. |
i@101 | 109 | local_start_time = timezone.localize(datetime.datetime.strptime(start_string, date_format)) |
i@101 | 110 | local_stop_time = timezone.localize(datetime.datetime.strptime(stop_string, date_format)) |
i@101 | 111 | |
i@101 | 112 | # Only save UTC time. |
i@101 | 113 | self.start_time = local_start_time.astimezone(pytz.utc) |
i@101 | 114 | self.stop_time = local_stop_time.astimezone(pytz.utc) |
i@101 | 115 | |
binietoglou@0 | 116 | self.latitude = float(raw_info['Latitude']) |
binietoglou@0 | 117 | self.longitude = float(raw_info['Longtitude']) |
ulalume3@47 | 118 | |
binietoglou@0 | 119 | # Read the rest of the header. |
binietoglou@0 | 120 | for c1 in range(int(raw_info['DataSets'])): |
binietoglou@0 | 121 | channel_info.append(match_lines(f.readline(), licel_file_channel_format)) |
ulalume3@47 | 122 | |
binietoglou@33 | 123 | self.raw_info = raw_info |
binietoglou@33 | 124 | self.channel_info = channel_info |
ulalume3@47 | 125 | |
ulalume3@27 | 126 | def duration(self): |
ulalume3@27 | 127 | """ Return the duration of the file. """ |
ulalume3@27 | 128 | dt = self.stop_time - self.start_time |
ulalume3@27 | 129 | return dt.seconds |
ulalume3@47 | 130 | |
ulalume3@47 | 131 | |
binietoglou@0 | 132 | class LicelFileChannel: |
ulalume3@47 | 133 | def __init__(self, raw_info=None, raw_data=None, duration=None, use_id_as_name=False): |
binietoglou@0 | 134 | self.raw_info = raw_info |
binietoglou@0 | 135 | self.raw_data = raw_data |
ulalume3@27 | 136 | self.duration = duration |
ulalume3@47 | 137 | self.use_id_as_name = use_id_as_name |
ulalume3@47 | 138 | |
binietoglou@0 | 139 | @property |
binietoglou@0 | 140 | def wavelength(self): |
binietoglou@0 | 141 | if self.raw_info is not None: |
binietoglou@0 | 142 | wave_str = self.raw_info['Wavelength'] |
binietoglou@0 | 143 | wavelength = wave_str.split('.')[0] |
binietoglou@0 | 144 | return int(wavelength) |
binietoglou@0 | 145 | else: |
binietoglou@0 | 146 | return None |
ulalume3@47 | 147 | |
binietoglou@0 | 148 | @property |
binietoglou@0 | 149 | def channel_name(self): |
binietoglou@0 | 150 | ''' |
binietoglou@0 | 151 | Construct the channel name adding analog photon info to avoid duplicates |
ulalume3@47 | 152 | |
ulalume3@47 | 153 | If use_id_as_name is True, the channel name will be the transient digitizer ID (e.g. BT01). |
ulalume3@47 | 154 | This could be useful if the lidar system has multiple telescopes, so the descriptive name is |
ulalume3@47 | 155 | not unique. |
binietoglou@0 | 156 | ''' |
ulalume3@47 | 157 | if self.use_id_as_name: |
ulalume3@47 | 158 | channel_name = self.raw_info['ID'] |
ulalume3@47 | 159 | else: |
ulalume3@47 | 160 | acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton']) |
ulalume3@47 | 161 | channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type) |
binietoglou@0 | 162 | return channel_name |
ulalume3@47 | 163 | |
binietoglou@0 | 164 | def analog_photon_string(self, analog_photon_number): |
binietoglou@0 | 165 | if analog_photon_number == '0': |
binietoglou@0 | 166 | string = 'an' |
binietoglou@0 | 167 | else: |
binietoglou@0 | 168 | string = 'ph' |
binietoglou@0 | 169 | return string |
ulalume3@47 | 170 | |
binietoglou@0 | 171 | def calculate_physical(self): |
binietoglou@0 | 172 | data = self.raw_data |
binietoglou@0 | 173 | |
binietoglou@0 | 174 | number_of_shots = float(self.raw_info['NShots']) |
ulalume3@27 | 175 | norm = data / number_of_shots |
binietoglou@0 | 176 | dz = float(self.raw_info['BinW']) |
ulalume3@47 | 177 | |
ulalume3@47 | 178 | if self.raw_info['AnalogPhoton'] == '0': |
binietoglou@0 | 179 | # If the channel is in analog mode |
binietoglou@0 | 180 | ADCb = int(self.raw_info['ADCbits']) |
ulalume3@47 | 181 | ADCrange = float(self.raw_info['Discriminator']) * 1000 # Value in mV |
ulalume3@47 | 182 | channel_data = norm * ADCrange / ((2 ** ADCb) - 1) |
ulalume3@47 | 183 | |
binietoglou@0 | 184 | # print ADCb, ADCRange,cdata,norm |
binietoglou@0 | 185 | else: |
binietoglou@0 | 186 | # If the channel is in photoncounting mode |
binietoglou@0 | 187 | # Frequency deduced from range resolution! (is this ok?) |
binietoglou@0 | 188 | # c = 300 # The result will be in MHZ |
binietoglou@0 | 189 | # SR = c/(2*dz) # To account for pulse folding |
binietoglou@0 | 190 | # channel_data = norm*SR |
binietoglou@0 | 191 | # CHANGE: |
binietoglou@0 | 192 | # For the SCC the data are needed in photons |
ulalume3@47 | 193 | channel_data = norm * number_of_shots |
ulalume3@47 | 194 | # print res,c,cdata,norm |
binietoglou@0 | 195 | |
ulalume3@47 | 196 | # Calculate Z |
ulalume3@47 | 197 | number_of_bins = int(self.raw_info['DataPoints']) |
ulalume3@47 | 198 | self.z = np.array([dz * bin_number + dz / 2.0 for bin_number in range(number_of_bins)]) |
binietoglou@0 | 199 | self.dz = dz |
binietoglou@0 | 200 | self.number_of_bins = number_of_bins |
binietoglou@0 | 201 | self.data = channel_data |
ioannis@22 | 202 | |
ulalume3@47 | 203 | |
binietoglou@0 | 204 | class LicelLidarMeasurement(BaseLidarMeasurement): |
binietoglou@0 | 205 | ''' |
binietoglou@0 | 206 | |
binietoglou@0 | 207 | ''' |
ulalume3@47 | 208 | raw_info = {} # Keep the raw info from the files |
ulalume3@27 | 209 | durations = {} # Keep the duration of the files |
victor@84 | 210 | laser_shots = [] |
ulalume3@47 | 211 | |
i@101 | 212 | def __init__(self, file_list=None, use_id_as_name=False, licel_timezone='UTC'): |
ulalume3@47 | 213 | self.use_id_as_name = use_id_as_name |
i@101 | 214 | self.licel_timezone = licel_timezone |
ulalume3@92 | 215 | super(LicelLidarMeasurement, self).__init__(file_list) |
ulalume3@47 | 216 | |
ulalume3@92 | 217 | def _import_file(self, filename): |
binietoglou@0 | 218 | if filename in self.files: |
ioannis@55 | 219 | logging.warning("File has been imported already: %s" % filename) |
binietoglou@0 | 220 | else: |
i@101 | 221 | current_file = LicelFile(filename, use_id_as_name=self.use_id_as_name, licel_timezone=self.licel_timezone) |
ioannis@22 | 222 | self.raw_info[current_file.filename] = current_file.raw_info |
ulalume3@27 | 223 | self.durations[current_file.filename] = current_file.duration() |
i@104 | 224 | |
victor@84 | 225 | file_laser_shots = [] |
ulalume3@47 | 226 | |
binietoglou@0 | 227 | for channel_name, channel in current_file.channels.items(): |
binietoglou@0 | 228 | if channel_name not in self.channels: |
ulalume3@27 | 229 | self.channels[channel_name] = LicelChannel(channel) |
ulalume3@47 | 230 | self.channels[channel_name].data[current_file.start_time] = channel.data |
victor@84 | 231 | file_laser_shots.append(channel.raw_info['NShots']) |
i@104 | 232 | |
victor@84 | 233 | self.laser_shots.append(file_laser_shots) |
binietoglou@0 | 234 | self.files.append(current_file.filename) |
ulalume3@47 | 235 | |
binietoglou@0 | 236 | def append(self, other): |
binietoglou@0 | 237 | |
binietoglou@0 | 238 | self.start_times.extend(other.start_times) |
binietoglou@0 | 239 | self.stop_times.extend(other.stop_times) |
binietoglou@0 | 240 | |
binietoglou@0 | 241 | for channel_name, channel in self.channels.items(): |
binietoglou@0 | 242 | channel.append(other.channels[channel_name]) |
binietoglou@0 | 243 | |
ioannis@22 | 244 | def _get_duration(self, raw_start_in_seconds): |
ulalume3@92 | 245 | """ Return the duration for a given time scale. If only a single |
ioannis@22 | 246 | file is imported, then this cannot be guessed from the time difference |
ioannis@22 | 247 | and the raw_info of the file are checked. |
ulalume3@92 | 248 | """ |
ulalume3@47 | 249 | |
ulalume3@47 | 250 | if len(raw_start_in_seconds) == 1: # If only one file imported |
ulalume3@47 | 251 | duration = self.durations.itervalues().next() # Get the first (and only) raw_info |
ioannis@31 | 252 | duration_sec = duration |
ioannis@22 | 253 | else: |
ioannis@22 | 254 | duration_sec = np.diff(raw_start_in_seconds)[0] |
binietoglou@0 | 255 | |
ioannis@22 | 256 | return duration_sec |
i@104 | 257 | |
ulalume3@92 | 258 | def get_custom_channel_parameters(self): |
victor@84 | 259 | params = [{ |
i@104 | 260 | "name": "DAQ_Range", |
i@104 | 261 | "dimensions": ('channels',), |
i@104 | 262 | "type": 'd', |
i@104 | 263 | "values": [self.channels[x].raw_info['Discriminator'] for x in self.channels.keys()] |
i@104 | 264 | }, { |
i@104 | 265 | "name": "LR_Input", |
i@104 | 266 | "dimensions": ('channels',), |
i@104 | 267 | "type": 'i', |
i@104 | 268 | "values": [self.channels[x].raw_info['LaserUsed'] for x in self.channels.keys()] |
i@104 | 269 | }, { |
i@104 | 270 | "name": "Laser_Shots", |
i@104 | 271 | "dimensions": ('time', 'channels',), |
i@104 | 272 | "type": 'i', |
i@104 | 273 | "values": self.laser_shots |
i@104 | 274 | }, |
victor@84 | 275 | ] |
i@104 | 276 | |
victor@84 | 277 | return params |
i@104 | 278 | |
ulalume3@92 | 279 | def get_custom_general_parameters(self): |
victor@87 | 280 | params = [{ |
i@104 | 281 | "name": "Altitude_meter_asl", |
i@104 | 282 | "value": self.raw_info[self.files[0]]["Altitude"] |
i@104 | 283 | }, { |
i@104 | 284 | "name": "Latitude_degrees_north", |
i@104 | 285 | "value": self.raw_info[self.files[0]]["Latitude"] |
i@104 | 286 | }, { |
i@104 | 287 | "name": "Longitude_degrees_east", |
i@104 | 288 | "value": self.raw_info[self.files[0]]["Longtitude"] |
i@104 | 289 | }, |
victor@87 | 290 | ] |
victor@87 | 291 | |
victor@87 | 292 | return params |
ioannis@22 | 293 | |
ulalume3@47 | 294 | |
ulalume3@27 | 295 | class LicelChannel(LidarChannel): |
binietoglou@0 | 296 | def __init__(self, channel_file): |
binietoglou@33 | 297 | c = 299792458.0 # Speed of light |
binietoglou@0 | 298 | self.wavelength = channel_file.wavelength |
binietoglou@0 | 299 | self.name = channel_file.channel_name |
ulalume3@27 | 300 | self.binwidth = channel_file.dz * 2 / c # in seconds |
binietoglou@0 | 301 | self.data = {} |
binietoglou@0 | 302 | self.resolution = channel_file.dz |
ulalume3@47 | 303 | self.z = np.arange( |
ulalume3@47 | 304 | channel_file.number_of_bins) * self.resolution + self.resolution / 2.0 # Change: add half bin in the z |
binietoglou@0 | 305 | self.points = channel_file.number_of_bins |
binietoglou@0 | 306 | self.rc = [] |
ulalume3@47 | 307 | self.duration = channel_file.duration |
victor@84 | 308 | self.raw_info = channel_file.raw_info |
i@104 | 309 | |
binietoglou@0 | 310 | def append(self, other): |
binietoglou@0 | 311 | if self.info != other.info: |
binietoglou@0 | 312 | raise ValueError('Channel info are different. Data can not be combined.') |
binietoglou@0 | 313 | |
binietoglou@0 | 314 | self.data = np.vstack([self.data, other.data]) |
binietoglou@0 | 315 | |
binietoglou@0 | 316 | def __unicode__(self): |
ioannis@26 | 317 | return "<Licel channel: %s>" % self.name |
ulalume3@47 | 318 | |
binietoglou@0 | 319 | def __str__(self): |
binietoglou@0 | 320 | return unicode(self).encode('utf-8') |
binietoglou@0 | 321 | |
ioannis@22 | 322 | |
ulalume3@47 | 323 | def match_lines(f1, f2): |
binietoglou@0 | 324 | list1 = f1.split() |
binietoglou@0 | 325 | list2 = f2.split() |
ulalume3@47 | 326 | |
binietoglou@0 | 327 | if len(list1) != len(list2): |
ioannis@55 | 328 | logging.debug("Channel parameter list has different length from licel specifications.") |
ioannis@55 | 329 | logging.debug("List 1: %s" % list1) |
ioannis@55 | 330 | logging.debug("List 2: %s" % list2) |
ulalume3@47 | 331 | combined = zip(list2, list1) |
binietoglou@0 | 332 | combined = dict(combined) |
binietoglou@0 | 333 | return combined |