generic.py

changeset 0
9d2b98ecf23d
child 1
82b144ee09b2
equal deleted inserted replaced
-1:000000000000 0:9d2b98ecf23d
1 # General imports
2 import datetime
3 from operator import itemgetter
4
5 # Science imports
6 import numpy as np
7 import matplotlib as mpl
8 from matplotlib import pyplot as plt
9 import netCDF4 as netcdf
10
11 # CNR-IMAA specific imports
12 import milos
13
14
15 netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4'
16
17 class BaseLidarMeasurement():
18 """ This is the general measurement object.
19 It is meant to become a general measurement object
20 independent of the input files.
21
22 Each subclass should implement the following:
23 * the import_file method.
24 * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters.
25
26 You can override the get_PT method to define a custom procedure to get ground temperature and pressure.
27 The one implemented by default is by using the MILOS meteorological station data.
28
29 """
30
31 def __init__(self, filelist= None):
32 self.info = {}
33 self.dimensions = {}
34 self.variables = {}
35 self.channels = {}
36 self.attributes = {}
37 self.files = []
38 self.dark_measurement = None
39 if filelist:
40 self.import_files(filelist)
41
42 def import_files(self,filelist):
43 for f in filelist:
44 self.import_file(f)
45 self.update()
46
47 def import_file(self,filename):
48 raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.')
49
50 def update(self):
51 '''
52 Update the the info, variables and dimensions of the lidar measurement based
53 on the information found in the channels.
54
55 Reading of the scan_angles parameter is not implemented.
56 '''
57
58 # Initialize
59 start_time =[]
60 stop_time = []
61 points = []
62 all_time_scales = []
63 channel_name_list = []
64
65 # Get the information from all the channels
66 for channel_name, channel in self.channels.items():
67 channel.update()
68 start_time.append(channel.start_time)
69 stop_time.append(channel.stop_time)
70 points.append(channel.points)
71 all_time_scales.append(channel.time)
72 channel_name_list.append(channel_name)
73
74 # Find the unique time scales used in the channels
75 time_scales = set(all_time_scales)
76
77 # Update the info dictionary
78 self.info['start_time'] = min(start_time)
79 self.info['stop_time'] = max(stop_time)
80 self.info['duration'] = self.info['stop_time'] - self.info['start_time']
81
82 # Update the dimensions dictionary
83 self.dimensions['points'] = max(points)
84 self.dimensions['channels'] = len(self.channels)
85 # self.dimensions['scan angles'] = 1
86 self.dimensions['nb_of_time_scales'] = len(time_scales)
87
88 # Update the variables dictionary
89 # Write time scales in seconds
90 raw_Data_Start_Time = []
91 raw_Data_Stop_Time = []
92
93 for current_time_scale in list(time_scales):
94 raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time
95 raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds
96 raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list
97 # Check if this time scale has measurements every 30 or 60 seconds.
98
99 duration = self._get_duration(raw_start_in_seconds)
100
101 raw_stop_in_seconds = raw_start_in_seconds + duration
102 raw_Data_Stop_Time.append(raw_stop_in_seconds)
103
104 self.variables['Raw_Data_Start_Time']= raw_Data_Start_Time
105 self.variables['Raw_Data_Stop_Time']= raw_Data_Stop_Time
106
107 # Make a dictionary to match time scales and channels
108 channel_timescales = []
109 for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales):
110 # The following lines are PEARL specific. The reason they are here is not clear.
111 # if channel_name =='1064BLR':
112 # channel_name = '1064'
113 for (ts,n) in zip(time_scales, range(len(time_scales))):
114 if current_time_scale == ts:
115 channel_timescales.append([channel_name,n])
116 self.variables['id_timescale'] = dict(channel_timescales)
117
118 def _get_duration(self, raw_start_in_seconds):
119 ''' Return the duration for a given time scale. In some files (ex. Licel) this
120 can be specified from the files themselves. In others this must be guessed.
121
122 '''
123 # The old method, kept here for reference
124 #dt = np.mean(np.diff(raw_start_in_seconds))
125 #for d in duration_list:
126 # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds
127 # duration = d
128
129 duration = np.diff(raw_start_in_seconds)[0]
130
131 return duration
132
133 def subset_by_channels(self, channel_subset):
134 ''' Get a measurement object containing only the channels with names
135 contained in the channel_sublet list '''
136
137 m = self.__class__() # Create an object of the same type as this one.
138 m.channels = dict([(channel, self.channels[channel]) for channel
139 in channel_subset])
140 m.update()
141 return m
142
143 def subset_by_time(self, start_time, stop_time):
144
145 if start_time > stop_time:
146 raise ValueError('Stop time should be after start time')
147
148 if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']):
149 raise ValueError('The time interval specified is not part of the measurement')
150
151 m = self.__class__() # Create an object of the same type as this one.
152 for (channel_name, channel) in self.channels.items():
153 m.channels[channel_name] = channel.subset_by_time(start_time, stop_time)
154 m.update()
155 return m
156
157 def r_plot(self):
158 #Make a basic plot of the data.
159 #Should include some dictionary with params to make plot stable.
160 pass
161
162 def r_pdf(self):
163 # Create a pdf report using a basic plot and meta-data.
164 pass
165
166 def save(self):
167 #Save the current state of the object to continue the analysis later.
168 pass
169
170 def get_PT(self):
171 ''' Sets the pressure and temperature at station level .
172 The results are stored in the info dictionary.
173 '''
174
175 self.info['Temperature'] = 10.0
176 self.info['Pressure'] = 930.0
177
178
179 def save_as_netcdf(self, filename):
180 """Saves the measurement in the netcdf format as required by the SCC.
181 Input: filename
182 """
183 params = self.extra_netcdf_parameters
184 needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure']
185
186 for parameter in needed_parameters:
187 stored_value = self.info.get(parameter, None)
188 if stored_value is None:
189 raise ValueError('A value needs to be specified for %s' % parameter)
190
191
192 dimensions = {'points': 1,
193 'channels': 1,
194 'time': None,
195 'nb_of_time_scales': 1,
196 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented
197
198 global_att = {'Measurement_ID': None,
199 'RawData_Start_Date': None,
200 'RawData_Start_Time_UT': None,
201 'RawData_Stop_Time_UT': None,
202 'RawBck_Start_Date': None,
203 'RawBck_Start_Time_UT': None,
204 'RawBck_Stop_Time_UT': None,
205 'Sounding_File_Name': None,
206 'LR_File_Name': None,
207 'Overlap_File_Name': None,
208 'Location': None,
209 'System': None,
210 'Latitude_degrees_north': None,
211 'Longitude_degrees_east': None,
212 'Altitude_meter_asl': None}
213
214 channel_variables = \
215 {'channel_ID': (('channels', ), 'i'),
216 'Background_Low': (('channels', ), 'd'),
217 'Background_High': (('channels', ), 'd'),
218 'LR_Input': (('channels', ), 'i'),
219 'DAQ_Range': (('channels', ), 'd'),
220 'Depolarization_Factor': (('channels', ), 'd'), }
221
222
223 channels = self.channels.keys()
224
225 input_values = dict(self.dimensions, **self.variables)
226
227 # Add some mandatory global attributes
228 input_values['Measurement_ID'] = self.info['Measurement_ID']
229 input_values['RawData_Start_Date'] = '\'%s\'' % self.info['start_time'].strftime('%Y%m%d')
230 input_values['RawData_Start_Time_UT'] = '\'%s\'' % self.info['start_time'].strftime('%H%M%S')
231 input_values['RawData_Stop_Time_UT'] = '\'%s\'' % self.info['stop_time'].strftime('%H%M%S')
232
233 # Add some optional global attributes
234 input_values['System'] = params.general_parameters['System']
235 input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north']
236 input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east']
237 input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl']
238
239 # Open a netCDF4 file
240 f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file
241
242 # Create the dimensions in the file
243 for (d,v) in dimensions.iteritems():
244 v = input_values.pop(d, v)
245 f.createDimension(d,v)
246
247 # Create global attributes
248 for (attrib,value) in global_att.iteritems():
249 val = input_values.pop(attrib,value)
250 if val:
251 exec('f.%s = %s' % (attrib,val))
252
253 """ Variables """
254 # Write the values of fixes channel parameters
255 for (var,t) in channel_variables.iteritems():
256 temp_v = f.createVariable(var,t[1],t[0])
257 for (channel, n) in zip(channels, range(len(channels))):
258 temp_v[n] = params.channel_parameters[channel][var]
259
260 # Write the id_timescale values
261 temp_id_timescale = f.createVariable('id_timescale','i',('channels',))
262 for (channel, n) in zip(channels, range(len(channels))):
263 temp_id_timescale[n] = self.variables['id_timescale'][channel]
264
265 # Laser pointing angle
266 temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',))
267 temp_v[:] = params.general_parameters['Laser_Pointing_Angle']
268
269 # Molecular calculation
270 temp_v = f.createVariable('Molecular_Calc','i')
271 temp_v[:] = params.general_parameters['Molecular_Calc']
272
273 # Laser pointing angles of profiles
274 temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales'))
275 for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'],
276 range(len(self.variables['Raw_Data_Start_Time']))):
277 temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle
278
279 # Raw data start/stop time
280 temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales'))
281 temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales'))
282 for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'],
283 self.variables['Raw_Data_Stop_Time'],
284 range(len(self.variables['Raw_Data_Start_Time']))):
285 temp_raw_start[:len(start_time),n] = start_time
286 temp_raw_stop[:len(stop_time),n] = stop_time
287
288 #Laser shots
289 temp_v = f.createVariable('Laser_Shots','i',('time','channels'))
290 for (channel,n) in zip(channels, range(len(channels))):
291 time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]])
292 temp_v[:time_length, n] = params.channel_parameters[channel]['Laser_Shots']
293
294 #Raw lidar data
295 temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points'))
296 for (channel,n) in zip(channels, range(len(channels))):
297 c = self.channels[channel]
298 temp_v[:len(c.time),n, :c.points] = c.matrix
299
300 self.add_dark_measurements_to_netcdf(f, channels)
301
302 #Pressure at lidar station
303 temp_v = f.createVariable('Pressure_at_Lidar_Station','d')
304 temp_v[:] = self.info['Pressure']
305
306 #Temperature at lidar station
307 temp_v = f.createVariable('Temperature_at_Lidar_Station','d')
308 temp_v[:] = self.info['Temperature']
309
310 self.save_netcdf_extra(f)
311 f.close()
312
313 def add_dark_measurements_to_netcdf(self, f, channels):
314
315 # Get dark measurements. If it is not given in self.dark_measurement
316 # try to get it using the get_dark_measurements method. If none is found
317 # return without adding something.
318 if self.dark_measurement is None:
319 self.dark_measurement = self.get_dark_measurements()
320
321 if self.dark_measurement is None:
322 return
323
324 dark_measurement = self.dark_measurement
325
326 # Calculate the length of the time_bck dimensions
327 number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()]
328 max_number_of_profiles = np.max(number_of_profiles)
329
330 # Create the dimension
331 f.createDimension('time_bck', max_number_of_profiles)
332
333 # Save the dark measurement data
334 temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points'))
335 for (channel,n) in zip(channels, range(len(channels))):
336 c = dark_measurement.channels[channel]
337 temp_v[:len(c.time),n, :c.points] = c.matrix
338
339 # Dark profile start/stop time
340 temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time','nb_of_time_scales'))
341 temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time','nb_of_time_scales'))
342 for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'],
343 dark_measurement.variables['Raw_Data_Stop_Time'],
344 range(len(dark_measurement.variables['Raw_Data_Start_Time']))):
345 temp_raw_start[:len(start_time),n] = start_time
346 temp_raw_stop[:len(stop_time),n] = stop_time
347
348 # Dark measurement start/stop time
349 f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d')
350 f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S')
351 f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S')
352
353
354
355 def save_netcdf_extra(self, f):
356 pass
357
358 def _gettime(self, date_str, time_str):
359 t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S')
360 return t
361
362 def plot(self):
363 for channel in self.channels:
364 self.channels[channel].plot(show_plot = False)
365 plt.show()
366
367 def get_dark_measurements(self):
368 return None
369
370
371 class Lidar_channel:
372
373 def __init__(self,channel_parameters):
374 c = 299792458 #Speed of light
375 self.wavelength = channel_parameters['name']
376 self.name = str(self.wavelength)
377 self.binwidth = float(channel_parameters['binwidth']) # in microseconds
378 self.data = {}
379 self.resolution = self.binwidth * c / 2
380 self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution/2.0 # Change: add half bin in the z
381 self.points = len(channel_parameters['data'])
382 self.rc = []
383 self.duration = 60
384
385 def calculate_rc(self):
386 background = np.mean(self.matrix[:,4000:], axis = 1) #Calculate the background from 30000m and above
387 self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2)
388
389
390 def update(self):
391 self.start_time = min(self.data.keys())
392 self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration)
393 self.time = tuple(sorted(self.data.keys()))
394 sorted_data = sorted(self.data.iteritems(), key=itemgetter(0))
395 self.matrix = np.array(map(itemgetter(1),sorted_data))
396
397 def _nearest_dt(self,dtime):
398 margin = datetime.timedelta(seconds = 300)
399 if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time):
400 print "Requested date not covered in this file"
401 raise
402 dt = abs(self.time - np.array(dtime))
403 dtmin = min(dt)
404
405 if dtmin > datetime.timedelta(seconds = 60):
406 print "Nearest profile more than 60 seconds away. dt = %s." % dtmin
407 ind_t = np.where(dt == dtmin)
408 ind_a= ind_t[0]
409 if len(ind_a) > 1:
410 ind_a = ind_a[0]
411 chosen_time = self.time[ind_a]
412 return chosen_time, ind_a
413
414 def subset_by_time(self, start_time, stop_time):
415
416 time_array = np.array(self.time)
417 condition = (time_array >= start_time) & (time_array <= stop_time)
418
419 subset_time = time_array[condition]
420 subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time])
421
422 #Create a list with the values needed by channel's __init__()
423 parameters_values = {'name': self.wavelength,
424 'binwidth': self.binwidth,
425 'data': subset_data[subset_time[0]],}
426
427 c = Lidar_channel(parameters_values)
428 c.data = subset_data
429 c.update()
430 return c
431
432
433 def profile(self,dtime, signal_type = 'rc'):
434 t, idx = self._nearest_dt(dtime)
435 if signal_type == 'rc':
436 data = self.rc
437 else:
438 data = self.matrix
439
440 prof = data[idx,:][0]
441 return prof, t
442
443 def get_slice(self, starttime, endtime, signal_type = 'rc'):
444 if signal_type == 'rc':
445 data = self.rc
446 else:
447 data = self.matrix
448 tim = np.array(self.time)
449 starttime = self._nearest_dt(starttime)[0]
450 endtime = self._nearest_dt(endtime)[0]
451 condition = (tim >= starttime) & (tim <= endtime)
452 sl = data[condition, :]
453 t = tim[condition]
454 return sl,t
455
456 def av_profile(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'):
457 starttime = tim - duration/2
458 endtime = tim + duration/2
459 d,t = self.get_slice(starttime, endtime, signal_type = signal_type)
460 prof = np.mean(d, axis = 0)
461 tmin = min(t)
462 tmax = max(t)
463 tav = tmin + (tmax-tmin)/2
464 return prof,(tav, tmin,tmax)
465
466 def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet):
467 #if filename is not None:
468 # matplotlib.use('Agg')
469
470 fig = plt.figure()
471 ax1 = fig.add_subplot(111)
472 self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom)
473 ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name))
474
475 if filename is not None:
476 pass
477 #plt.savefig(filename)
478 else:
479 if show_plot:
480 plt.show()
481 #plt.close() ???
482
483 def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', zoom = [0,12000,0,-1]):
484
485 if signal_type == 'rc':
486 if len(self.rc) == 0:
487 self.calculate_rc()
488 data = self.rc
489 else:
490 data = self.matrix
491
492 hmax_idx = self.index_at_height(zoom[1])
493
494 ax1.set_ylabel('Altitude (km)')
495 ax1.set_xlabel('Time UTC')
496 #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064
497 #dateFormatter = mpl.dates.DateFormatter('%H.%M')
498 #hourlocator = mpl.dates.HourLocator()
499
500 #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m')
501 #daylocator = mpl.dates.DayLocator()
502 hourFormatter = mpl.dates.DateFormatter('%H.%M')
503 hourlocator = mpl.dates.AutoDateLocator(interval_multiples=True)
504
505
506 #ax1.axes.xaxis.set_major_formatter(dayFormatter)
507 #ax1.axes.xaxis.set_major_locator(daylocator)
508 ax1.axes.xaxis.set_major_formatter(hourFormatter)
509 ax1.axes.xaxis.set_major_locator(hourlocator)
510
511
512 ts1 = mpl.dates.date2num(self.start_time)
513 ts2 = mpl.dates.date2num(self.stop_time)
514
515
516 im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]],
517 aspect = 'auto',
518 origin = 'lower',
519 cmap = cmap,
520 #vmin = 0,
521 vmin = data[:,10:400].max() * 0.1,
522 #vmax = 1.4*10**7,
523 vmax = data[:,10:400].max() * 0.9,
524 extent = [ts1,ts2,self.z[zoom[0]]/1000.0, self.z[hmax_idx]/1000.0],
525 )
526
527 cb1 = plt.colorbar(im1)
528 cb1.ax.set_ylabel('a.u.')
529
530 def index_at_height(self, height):
531 idx = np.array(np.abs(self.z - height).argmin())
532 if idx.size >1:
533 idx =idx[0]
534 return idx
535
536 def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID):
537 #Read the lidar files and select channels
538 temp_m = LidarClass(files)
539 m = temp_m.subset_by_channels(channels)
540 m.get_PT()
541 m.info['Measurement_ID'] = measurement_ID
542 m.save_as_netcdf(filename)
543

mercurial