1 import numpy as np |
1 import numpy as np |
2 import datetime |
2 import datetime |
3 from generic import BaseLidarMeasurement, Lidar_channel |
3 from generic import BaseLidarMeasurement, LidarChannel |
4 import musa_netcdf_parameters |
4 import musa_netcdf_parameters |
5 import musa_2009_netcdf_parameters |
5 import musa_2009_netcdf_parameters |
6 |
6 |
7 licel_file_header_format = ['Filename', |
7 licel_file_header_format = ['Filename', |
8 'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', |
8 'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', |
9 'LS1 Rate1 LS2 Rate2 DataSets', ] |
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' |
10 licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID' |
11 |
11 |
12 |
12 |
13 class LicelFile: |
13 class LicelFile: |
|
14 |
14 def __init__(self, filename): |
15 def __init__(self, filename): |
15 self.start_time = None |
16 self.start_time = None |
16 self.stop_time = None |
17 self.stop_time = None |
17 self.import_file(filename) |
18 self.import_file(filename) |
18 self.calculate_physical() |
19 self.calculate_physical() |
59 a = np.fromfile(f, 'b', 1) |
60 a = np.fromfile(f, 'b', 1) |
60 b = np.fromfile(f, 'b', 1) |
61 b = np.fromfile(f, 'b', 1) |
61 |
62 |
62 if (a[0] != 13) | (b[0] != 10): |
63 if (a[0] != 13) | (b[0] != 10): |
63 print "Warning: No end of line found after record. File could be corrupt" |
64 print "Warning: No end of line found after record. File could be corrupt" |
64 channel = LicelFileChannel(current_channel_info, raw_data) |
65 channel = LicelFileChannel(current_channel_info, raw_data, self.duration()) |
65 |
66 |
66 channel_name = channel.channel_name |
67 channel_name = channel.channel_name |
67 if channel_name in channels.keys(): |
68 if channel_name in channels.keys(): |
68 # If the analog/photon naming scheme is not enough, find a new one! |
69 # If the analog/photon naming scheme is not enough, find a new one! |
69 raise IOError('Trying to import two channels with the same name') |
70 raise IOError('Trying to import two channels with the same name') |
71 channels[channel_name] = channel |
72 channels[channel_name] = channel |
72 f.close() |
73 f.close() |
73 |
74 |
74 self.raw_info = raw_info |
75 self.raw_info = raw_info |
75 self.channels = channels |
76 self.channels = channels |
76 |
77 |
77 |
78 def duration(self): |
|
79 """ Return the duration of the file. """ |
|
80 dt = self.stop_time - self.start_time |
|
81 return dt.seconds |
|
82 |
|
83 |
78 class LicelFileChannel: |
84 class LicelFileChannel: |
79 |
85 |
80 def __init__(self, raw_info = None, raw_data = None): |
86 def __init__(self, raw_info = None, raw_data = None, duration = None): |
81 self.raw_info = raw_info |
87 self.raw_info = raw_info |
82 self.raw_data = raw_data |
88 self.raw_data = raw_data |
83 |
89 self.duration = duration |
|
90 |
84 @property |
91 @property |
85 def wavelength(self): |
92 def wavelength(self): |
86 if self.raw_info is not None: |
93 if self.raw_info is not None: |
87 wave_str = self.raw_info['Wavelength'] |
94 wave_str = self.raw_info['Wavelength'] |
88 wavelength = wave_str.split('.')[0] |
95 wavelength = wave_str.split('.')[0] |
108 |
115 |
109 def calculate_physical(self): |
116 def calculate_physical(self): |
110 data = self.raw_data |
117 data = self.raw_data |
111 |
118 |
112 number_of_shots = float(self.raw_info['NShots']) |
119 number_of_shots = float(self.raw_info['NShots']) |
113 norm = data/number_of_shots |
120 norm = data / number_of_shots |
114 dz = float(self.raw_info['BinW']) |
121 dz = float(self.raw_info['BinW']) |
115 |
122 |
116 if self.raw_info['AnalogPhoton']=='0': |
123 if self.raw_info['AnalogPhoton']=='0': |
117 # If the channel is in analog mode |
124 # If the channel is in analog mode |
118 ADCb = int(self.raw_info['ADCbits']) |
125 ADCb = int(self.raw_info['ADCbits']) |
126 # c = 300 # The result will be in MHZ |
133 # c = 300 # The result will be in MHZ |
127 # SR = c/(2*dz) # To account for pulse folding |
134 # SR = c/(2*dz) # To account for pulse folding |
128 # channel_data = norm*SR |
135 # channel_data = norm*SR |
129 # CHANGE: |
136 # CHANGE: |
130 # For the SCC the data are needed in photons |
137 # For the SCC the data are needed in photons |
131 channel_data = norm*number_of_shots |
138 channel_data = norm *number_of_shots |
132 #print res,c,cdata,norm |
139 #print res,c,cdata,norm |
133 |
140 |
134 #Calculate Z |
141 #Calculate Z |
135 number_of_bins = int(self.raw_info['DataPoints']) |
142 number_of_bins = int(self.raw_info['DataPoints']) |
136 self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) |
143 self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) |
141 |
148 |
142 class LicelLidarMeasurement(BaseLidarMeasurement): |
149 class LicelLidarMeasurement(BaseLidarMeasurement): |
143 ''' |
150 ''' |
144 |
151 |
145 ''' |
152 ''' |
146 |
|
147 extra_netcdf_parameters = musa_netcdf_parameters |
153 extra_netcdf_parameters = musa_netcdf_parameters |
148 raw_info = {} # Keep the raw info from the files |
154 raw_info = {} # Keep the raw info from the files |
|
155 durations = {} # Keep the duration of the files |
149 |
156 |
150 def import_file(self, filename): |
157 def import_file(self, filename): |
151 if filename in self.files: |
158 if filename in self.files: |
152 print "File has been imported already:" + filename |
159 print "File has been imported already:" + filename |
153 else: |
160 else: |
154 current_file = LicelFile(filename) |
161 current_file = LicelFile(filename) |
155 self.raw_info[current_file.filename] = current_file.raw_info |
162 self.raw_info[current_file.filename] = current_file.raw_info |
|
163 self.durations[current_file.filename] = current_file.duration() |
156 |
164 |
157 for channel_name, channel in current_file.channels.items(): |
165 for channel_name, channel in current_file.channels.items(): |
158 if channel_name not in self.channels: |
166 if channel_name not in self.channels: |
159 self.channels[channel_name] = Licel_channel(channel) |
167 self.channels[channel_name] = LicelChannel(channel) |
160 self.channels[channel_name].data[current_file.start_time] = channel.data |
168 self.channels[channel_name].data[current_file.start_time] = channel.data |
161 self.files.append(current_file.filename) |
169 self.files.append(current_file.filename) |
162 |
170 |
163 def append(self, other): |
171 def append(self, other): |
164 |
172 |
170 |
178 |
171 def _get_duration(self, raw_start_in_seconds): |
179 def _get_duration(self, raw_start_in_seconds): |
172 ''' Return the duration for a given time scale. If only a single |
180 ''' Return the duration for a given time scale. If only a single |
173 file is imported, then this cannot be guessed from the time difference |
181 file is imported, then this cannot be guessed from the time difference |
174 and the raw_info of the file are checked. |
182 and the raw_info of the file are checked. |
175 |
|
176 ''' |
183 ''' |
177 |
184 |
178 if len(raw_start_in_seconds) == 1: # If only one file imported |
185 if len(raw_start_in_seconds) == 1: # If only one file imported |
179 raw_info = self.raw_info.itervalues().next() # Get the first (and only) raw_info |
186 duration = self.durations.itervalues().next() # Get the first (and only) raw_info |
180 start_str = "%s %s" % (raw_info['StartDate'], raw_info['StartTime']) |
|
181 end_str = "%s %s" % (raw_info['EndDate'], raw_info['EndTime']) |
|
182 start_time = datetime.datetime.strptime(start_str, "%d/%m/%Y %H:%M:%S") |
|
183 end_time = datetime.datetime.strptime(end_str, "%d/%m/%Y %H:%M:%S") |
|
184 duration = end_time - start_time |
|
185 duration_sec = duration.seconds |
187 duration_sec = duration.seconds |
186 else: |
188 else: |
187 duration_sec = np.diff(raw_start_in_seconds)[0] |
189 duration_sec = np.diff(raw_start_in_seconds)[0] |
188 |
190 |
189 return duration_sec |
191 return duration_sec |
190 |
192 |
191 |
193 |
192 class Licel_channel(Lidar_channel): |
194 class LicelChannel(LidarChannel): |
193 |
195 |
194 def __init__(self, channel_file): |
196 def __init__(self, channel_file): |
195 c = 299792458.0 #Speed of light |
197 c = 299792458.0 # Speed of light |
196 self.wavelength = channel_file.wavelength |
198 self.wavelength = channel_file.wavelength |
197 self.name = channel_file.channel_name |
199 self.name = channel_file.channel_name |
198 self.binwidth = channel_file.dz * 2 / c # in microseconds |
200 self.binwidth = channel_file.dz * 2 / c # in seconds |
199 self.data = {} |
201 self.data = {} |
200 self.resolution = channel_file.dz |
202 self.resolution = channel_file.dz |
201 self.z = np.arange(channel_file.number_of_bins) *self.resolution + self.resolution/2.0 #Change: add half bin in the z |
203 self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z |
202 self.points = channel_file.number_of_bins |
204 self.points = channel_file.number_of_bins |
203 self.rc = [] |
205 self.rc = [] |
204 self.duration = 60 |
206 self.duration = channel_file.duration |
205 |
207 |
206 def append(self, other): |
208 def append(self, other): |
207 if self.info != other.info: |
209 if self.info != other.info: |
208 raise ValueError('Channel info are different. Data can not be combined.') |
210 raise ValueError('Channel info are different. Data can not be combined.') |
209 |
211 |