licel.py

Mon, 24 Nov 2014 11:35:17 +0200

author
Iannis <ioannis@inoe.ro>
date
Mon, 24 Nov 2014 11:35:17 +0200
changeset 29
870fc8f65eeb
parent 26
84f7d763deeb
child 28
73337ce10473
permissions
-rw-r--r--

Modiifed RALI paramters

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

mercurial