Fri, 13 Dec 2013 11:41:07 +0100
Improved labels in colorplots
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 | |
binietoglou@0 | 12 | class LicelFile: |
binietoglou@0 | 13 | def __init__(self, filename): |
binietoglou@0 | 14 | self.start_time = None |
binietoglou@0 | 15 | self.stop_time = None |
binietoglou@0 | 16 | self.import_file(filename) |
binietoglou@0 | 17 | self.calculate_physical() |
binietoglou@0 | 18 | self.filename = filename |
binietoglou@0 | 19 | |
binietoglou@0 | 20 | def calculate_physical(self): |
binietoglou@0 | 21 | for channel in self.channels.itervalues(): |
binietoglou@0 | 22 | channel.calculate_physical() |
binietoglou@0 | 23 | |
binietoglou@0 | 24 | def import_file(self, filename): |
binietoglou@0 | 25 | """Imports a licel file. |
binietoglou@0 | 26 | Input: filename |
binietoglou@0 | 27 | Output: object """ |
binietoglou@0 | 28 | |
binietoglou@0 | 29 | raw_info = {} |
binietoglou@0 | 30 | channels = {} |
binietoglou@0 | 31 | channel_info = [] |
binietoglou@0 | 32 | |
ulalume3@2 | 33 | f = open(filename, 'rb') |
binietoglou@0 | 34 | |
binietoglou@0 | 35 | #Read the first 3 lines of the header |
binietoglou@0 | 36 | raw_info = {} |
binietoglou@0 | 37 | for c1 in range(3): |
binietoglou@0 | 38 | raw_info.update(match_lines(f.readline(), licel_file_header_format[c1])) |
binietoglou@0 | 39 | |
binietoglou@0 | 40 | start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) |
binietoglou@0 | 41 | stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) |
binietoglou@0 | 42 | date_format = '%d/%m/%Y %H:%M:%S' |
binietoglou@0 | 43 | self.start_time = datetime.datetime.strptime(start_string, date_format) |
binietoglou@0 | 44 | self.stop_time = datetime.datetime.strptime(stop_string, date_format) |
binietoglou@0 | 45 | self.latitude = float(raw_info['Latitude']) |
binietoglou@0 | 46 | self.longitude = float(raw_info['Longtitude']) |
binietoglou@0 | 47 | |
binietoglou@0 | 48 | # Read the rest of the header. |
binietoglou@0 | 49 | for c1 in range(int(raw_info['DataSets'])): |
binietoglou@0 | 50 | channel_info.append(match_lines(f.readline(), licel_file_channel_format)) |
binietoglou@0 | 51 | |
binietoglou@0 | 52 | # Check the complete header is read |
binietoglou@0 | 53 | a = f.readline() |
binietoglou@0 | 54 | |
binietoglou@0 | 55 | # Import the data |
binietoglou@0 | 56 | for current_channel_info in channel_info: |
binietoglou@16 | 57 | raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints'])) |
binietoglou@0 | 58 | a = np.fromfile(f, 'b', 1) |
binietoglou@0 | 59 | b = np.fromfile(f, 'b', 1) |
ulalume3@2 | 60 | |
binietoglou@0 | 61 | if (a[0] != 13) | (b[0] != 10): |
binietoglou@0 | 62 | print "Warning: No end of line found after record. File could be corrupt" |
binietoglou@0 | 63 | channel = LicelFileChannel(current_channel_info, raw_data) |
binietoglou@0 | 64 | |
binietoglou@0 | 65 | channel_name = channel.channel_name |
binietoglou@0 | 66 | if channel_name in channels.keys(): |
binietoglou@0 | 67 | # If the analog/photon naming scheme is not enough, find a new one! |
binietoglou@0 | 68 | raise IOError('Trying to import two channels with the same name') |
binietoglou@0 | 69 | |
binietoglou@0 | 70 | channels[channel_name] = channel |
binietoglou@0 | 71 | f.close() |
binietoglou@0 | 72 | |
binietoglou@0 | 73 | self.raw_info = raw_info |
binietoglou@0 | 74 | self.channels = channels |
binietoglou@0 | 75 | |
binietoglou@0 | 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 |
binietoglou@0 | 140 | |
binietoglou@0 | 141 | class LicelLidarMeasurement(BaseLidarMeasurement): |
binietoglou@0 | 142 | ''' |
binietoglou@0 | 143 | |
binietoglou@0 | 144 | ''' |
binietoglou@0 | 145 | |
binietoglou@0 | 146 | extra_netcdf_parameters = musa_netcdf_parameters |
binietoglou@0 | 147 | |
binietoglou@0 | 148 | def import_file(self, filename): |
binietoglou@0 | 149 | if filename in self.files: |
binietoglou@0 | 150 | print "File has been imported already:" + filename |
binietoglou@0 | 151 | else: |
binietoglou@0 | 152 | current_file = LicelFile(filename) |
binietoglou@0 | 153 | for channel_name, channel in current_file.channels.items(): |
binietoglou@0 | 154 | if channel_name not in self.channels: |
binietoglou@0 | 155 | self.channels[channel_name] = Licel_channel(channel) |
binietoglou@0 | 156 | self.channels[channel_name].data[current_file.start_time] = channel.data |
binietoglou@0 | 157 | self.files.append(current_file.filename) |
binietoglou@0 | 158 | |
binietoglou@0 | 159 | def append(self, other): |
binietoglou@0 | 160 | |
binietoglou@0 | 161 | self.start_times.extend(other.start_times) |
binietoglou@0 | 162 | self.stop_times.extend(other.stop_times) |
binietoglou@0 | 163 | |
binietoglou@0 | 164 | for channel_name, channel in self.channels.items(): |
binietoglou@0 | 165 | channel.append(other.channels[channel_name]) |
binietoglou@0 | 166 | |
binietoglou@0 | 167 | |
binietoglou@0 | 168 | class Licel_channel(Lidar_channel): |
binietoglou@0 | 169 | |
binietoglou@0 | 170 | def __init__(self, channel_file): |
binietoglou@0 | 171 | c = 299792458.0 #Speed of light |
binietoglou@0 | 172 | self.wavelength = channel_file.wavelength |
binietoglou@0 | 173 | self.name = channel_file.channel_name |
binietoglou@0 | 174 | self.binwidth = channel_file.dz * 2 / c # in microseconds |
binietoglou@0 | 175 | self.data = {} |
binietoglou@0 | 176 | self.resolution = channel_file.dz |
binietoglou@0 | 177 | self.z = np.arange(channel_file.number_of_bins) *self.resolution + self.resolution/2.0 #Change: add half bin in the z |
binietoglou@0 | 178 | self.points = channel_file.number_of_bins |
binietoglou@0 | 179 | self.rc = [] |
binietoglou@0 | 180 | self.duration = 60 |
binietoglou@0 | 181 | |
binietoglou@0 | 182 | def append(self, other): |
binietoglou@0 | 183 | if self.info != other.info: |
binietoglou@0 | 184 | raise ValueError('Channel info are different. Data can not be combined.') |
binietoglou@0 | 185 | |
binietoglou@0 | 186 | self.data = np.vstack([self.data, other.data]) |
binietoglou@0 | 187 | |
binietoglou@0 | 188 | def __unicode__(self): |
binietoglou@0 | 189 | return "<Licel channel: %s>" % self.info['Wavelength'] |
binietoglou@0 | 190 | |
binietoglou@0 | 191 | def __str__(self): |
binietoglou@0 | 192 | return unicode(self).encode('utf-8') |
binietoglou@0 | 193 | |
binietoglou@0 | 194 | class Licel2009LidarMeasurement(LicelLidarMeasurement): |
binietoglou@0 | 195 | extra_netcdf_parameters = musa_2009_netcdf_parameters |
binietoglou@0 | 196 | |
binietoglou@0 | 197 | def match_lines(f1,f2): |
binietoglou@0 | 198 | list1 = f1.split() |
binietoglou@0 | 199 | list2 = f2.split() |
binietoglou@0 | 200 | |
binietoglou@0 | 201 | if len(list1) != len(list2): |
binietoglou@0 | 202 | print "Warning: Combining lists of different lengths." |
binietoglou@0 | 203 | print "List 1: %s" % list1 |
binietoglou@0 | 204 | print "List 2: %s" % list2 |
binietoglou@0 | 205 | combined = zip(list2,list1) |
binietoglou@0 | 206 | combined = dict(combined) |
binietoglou@0 | 207 | return combined |
binietoglou@0 | 208 |