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