binietoglou@0: # General imports binietoglou@0: import datetime binietoglou@0: from operator import itemgetter binietoglou@0: binietoglou@0: # Science imports binietoglou@0: import numpy as np binietoglou@0: import matplotlib as mpl binietoglou@0: from matplotlib import pyplot as plt binietoglou@0: import netCDF4 as netcdf binietoglou@0: binietoglou@1: netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4' binietoglou@0: binietoglou@0: binietoglou@0: class BaseLidarMeasurement(): binietoglou@0: """ This is the general measurement object. binietoglou@0: It is meant to become a general measurement object binietoglou@0: independent of the input files. binietoglou@0: binietoglou@0: Each subclass should implement the following: binietoglou@0: * the import_file method. binietoglou@0: * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters. binietoglou@0: binietoglou@0: You can override the get_PT method to define a custom procedure to get ground temperature and pressure. binietoglou@0: The one implemented by default is by using the MILOS meteorological station data. binietoglou@0: binietoglou@0: """ binietoglou@0: binietoglou@0: def __init__(self, filelist= None): binietoglou@0: self.info = {} binietoglou@0: self.dimensions = {} binietoglou@0: self.variables = {} binietoglou@0: self.channels = {} binietoglou@0: self.attributes = {} binietoglou@0: self.files = [] binietoglou@0: self.dark_measurement = None binietoglou@0: if filelist: binietoglou@0: self.import_files(filelist) binietoglou@0: binietoglou@0: def import_files(self,filelist): binietoglou@0: for f in filelist: binietoglou@0: self.import_file(f) binietoglou@0: self.update() binietoglou@0: binietoglou@0: def import_file(self,filename): binietoglou@0: raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.') binietoglou@0: binietoglou@0: def update(self): binietoglou@0: ''' binietoglou@0: Update the the info, variables and dimensions of the lidar measurement based binietoglou@0: on the information found in the channels. binietoglou@0: binietoglou@0: Reading of the scan_angles parameter is not implemented. binietoglou@0: ''' binietoglou@0: binietoglou@0: # Initialize binietoglou@0: start_time =[] binietoglou@0: stop_time = [] binietoglou@0: points = [] binietoglou@0: all_time_scales = [] binietoglou@0: channel_name_list = [] binietoglou@0: binietoglou@0: # Get the information from all the channels binietoglou@0: for channel_name, channel in self.channels.items(): binietoglou@0: channel.update() binietoglou@0: start_time.append(channel.start_time) binietoglou@0: stop_time.append(channel.stop_time) binietoglou@0: points.append(channel.points) binietoglou@0: all_time_scales.append(channel.time) binietoglou@0: channel_name_list.append(channel_name) binietoglou@0: binietoglou@0: # Find the unique time scales used in the channels binietoglou@0: time_scales = set(all_time_scales) binietoglou@0: binietoglou@0: # Update the info dictionary binietoglou@0: self.info['start_time'] = min(start_time) binietoglou@0: self.info['stop_time'] = max(stop_time) binietoglou@0: self.info['duration'] = self.info['stop_time'] - self.info['start_time'] binietoglou@0: binietoglou@0: # Update the dimensions dictionary binietoglou@0: self.dimensions['points'] = max(points) binietoglou@0: self.dimensions['channels'] = len(self.channels) binietoglou@0: # self.dimensions['scan angles'] = 1 binietoglou@0: self.dimensions['nb_of_time_scales'] = len(time_scales) binietoglou@0: binietoglou@0: # Update the variables dictionary binietoglou@0: # Write time scales in seconds binietoglou@0: raw_Data_Start_Time = [] binietoglou@0: raw_Data_Stop_Time = [] binietoglou@0: binietoglou@0: for current_time_scale in list(time_scales): binietoglou@0: raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time binietoglou@0: raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds binietoglou@0: raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list binietoglou@0: # Check if this time scale has measurements every 30 or 60 seconds. binietoglou@0: binietoglou@0: duration = self._get_duration(raw_start_in_seconds) binietoglou@0: binietoglou@0: raw_stop_in_seconds = raw_start_in_seconds + duration binietoglou@0: raw_Data_Stop_Time.append(raw_stop_in_seconds) binietoglou@0: binietoglou@0: self.variables['Raw_Data_Start_Time']= raw_Data_Start_Time binietoglou@0: self.variables['Raw_Data_Stop_Time']= raw_Data_Stop_Time binietoglou@0: binietoglou@0: # Make a dictionary to match time scales and channels binietoglou@0: channel_timescales = [] binietoglou@0: for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales): binietoglou@0: # The following lines are PEARL specific. The reason they are here is not clear. binietoglou@0: # if channel_name =='1064BLR': binietoglou@0: # channel_name = '1064' binietoglou@0: for (ts,n) in zip(time_scales, range(len(time_scales))): binietoglou@0: if current_time_scale == ts: binietoglou@0: channel_timescales.append([channel_name,n]) binietoglou@0: self.variables['id_timescale'] = dict(channel_timescales) binietoglou@0: binietoglou@0: def _get_duration(self, raw_start_in_seconds): binietoglou@0: ''' Return the duration for a given time scale. In some files (ex. Licel) this binietoglou@0: can be specified from the files themselves. In others this must be guessed. binietoglou@0: binietoglou@0: ''' binietoglou@0: # The old method, kept here for reference binietoglou@0: #dt = np.mean(np.diff(raw_start_in_seconds)) binietoglou@0: #for d in duration_list: binietoglou@0: # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds binietoglou@0: # duration = d binietoglou@0: binietoglou@0: duration = np.diff(raw_start_in_seconds)[0] binietoglou@0: binietoglou@0: return duration binietoglou@0: binietoglou@0: def subset_by_channels(self, channel_subset): binietoglou@0: ''' Get a measurement object containing only the channels with names binietoglou@0: contained in the channel_sublet list ''' binietoglou@0: binietoglou@0: m = self.__class__() # Create an object of the same type as this one. binietoglou@0: m.channels = dict([(channel, self.channels[channel]) for channel binietoglou@0: in channel_subset]) binietoglou@0: m.update() binietoglou@0: return m binietoglou@0: binietoglou@0: def subset_by_time(self, start_time, stop_time): binietoglou@0: binietoglou@0: if start_time > stop_time: binietoglou@0: raise ValueError('Stop time should be after start time') binietoglou@0: binietoglou@0: if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']): binietoglou@0: raise ValueError('The time interval specified is not part of the measurement') binietoglou@0: binietoglou@0: m = self.__class__() # Create an object of the same type as this one. binietoglou@0: for (channel_name, channel) in self.channels.items(): binietoglou@0: m.channels[channel_name] = channel.subset_by_time(start_time, stop_time) binietoglou@0: m.update() binietoglou@0: return m binietoglou@0: binietoglou@0: def r_plot(self): binietoglou@0: #Make a basic plot of the data. binietoglou@0: #Should include some dictionary with params to make plot stable. binietoglou@0: pass binietoglou@0: binietoglou@0: def r_pdf(self): binietoglou@0: # Create a pdf report using a basic plot and meta-data. binietoglou@0: pass binietoglou@0: binietoglou@0: def save(self): binietoglou@0: #Save the current state of the object to continue the analysis later. binietoglou@0: pass binietoglou@0: binietoglou@0: def get_PT(self): binietoglou@0: ''' Sets the pressure and temperature at station level . binietoglou@0: The results are stored in the info dictionary. binietoglou@0: ''' binietoglou@0: binietoglou@0: self.info['Temperature'] = 10.0 binietoglou@0: self.info['Pressure'] = 930.0 binietoglou@0: binietoglou@0: binietoglou@0: def save_as_netcdf(self, filename): binietoglou@0: """Saves the measurement in the netcdf format as required by the SCC. binietoglou@0: Input: filename binietoglou@0: """ binietoglou@0: params = self.extra_netcdf_parameters binietoglou@0: needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure'] binietoglou@0: binietoglou@0: for parameter in needed_parameters: binietoglou@0: stored_value = self.info.get(parameter, None) binietoglou@0: if stored_value is None: binietoglou@0: raise ValueError('A value needs to be specified for %s' % parameter) binietoglou@0: binietoglou@0: binietoglou@0: dimensions = {'points': 1, binietoglou@0: 'channels': 1, binietoglou@0: 'time': None, binietoglou@0: 'nb_of_time_scales': 1, binietoglou@0: 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented binietoglou@0: binietoglou@0: global_att = {'Measurement_ID': None, binietoglou@0: 'RawData_Start_Date': None, binietoglou@0: 'RawData_Start_Time_UT': None, binietoglou@0: 'RawData_Stop_Time_UT': None, binietoglou@0: 'RawBck_Start_Date': None, binietoglou@0: 'RawBck_Start_Time_UT': None, binietoglou@0: 'RawBck_Stop_Time_UT': None, binietoglou@0: 'Sounding_File_Name': None, binietoglou@0: 'LR_File_Name': None, binietoglou@0: 'Overlap_File_Name': None, binietoglou@0: 'Location': None, binietoglou@0: 'System': None, binietoglou@0: 'Latitude_degrees_north': None, binietoglou@0: 'Longitude_degrees_east': None, binietoglou@0: 'Altitude_meter_asl': None} binietoglou@0: binietoglou@0: channel_variables = \ binietoglou@0: {'channel_ID': (('channels', ), 'i'), binietoglou@0: 'Background_Low': (('channels', ), 'd'), binietoglou@0: 'Background_High': (('channels', ), 'd'), binietoglou@0: 'LR_Input': (('channels', ), 'i'), binietoglou@0: 'DAQ_Range': (('channels', ), 'd'), binietoglou@0: 'Depolarization_Factor': (('channels', ), 'd'), } binietoglou@0: binietoglou@0: binietoglou@0: channels = self.channels.keys() binietoglou@0: binietoglou@0: input_values = dict(self.dimensions, **self.variables) binietoglou@0: binietoglou@0: # Add some mandatory global attributes binietoglou@0: input_values['Measurement_ID'] = self.info['Measurement_ID'] binietoglou@0: input_values['RawData_Start_Date'] = '\'%s\'' % self.info['start_time'].strftime('%Y%m%d') binietoglou@0: input_values['RawData_Start_Time_UT'] = '\'%s\'' % self.info['start_time'].strftime('%H%M%S') binietoglou@0: input_values['RawData_Stop_Time_UT'] = '\'%s\'' % self.info['stop_time'].strftime('%H%M%S') binietoglou@0: binietoglou@0: # Add some optional global attributes binietoglou@0: input_values['System'] = params.general_parameters['System'] binietoglou@0: input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north'] binietoglou@0: input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east'] binietoglou@0: input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl'] binietoglou@0: binietoglou@0: # Open a netCDF4 file binietoglou@0: f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file binietoglou@0: binietoglou@0: # Create the dimensions in the file binietoglou@0: for (d,v) in dimensions.iteritems(): binietoglou@0: v = input_values.pop(d, v) binietoglou@0: f.createDimension(d,v) binietoglou@0: binietoglou@0: # Create global attributes binietoglou@0: for (attrib,value) in global_att.iteritems(): binietoglou@0: val = input_values.pop(attrib,value) binietoglou@0: if val: binietoglou@0: exec('f.%s = %s' % (attrib,val)) binietoglou@0: binietoglou@0: """ Variables """ binietoglou@0: # Write the values of fixes channel parameters binietoglou@0: for (var,t) in channel_variables.iteritems(): binietoglou@0: temp_v = f.createVariable(var,t[1],t[0]) binietoglou@0: for (channel, n) in zip(channels, range(len(channels))): binietoglou@0: temp_v[n] = params.channel_parameters[channel][var] binietoglou@0: binietoglou@0: # Write the id_timescale values binietoglou@0: temp_id_timescale = f.createVariable('id_timescale','i',('channels',)) binietoglou@0: for (channel, n) in zip(channels, range(len(channels))): binietoglou@0: temp_id_timescale[n] = self.variables['id_timescale'][channel] binietoglou@0: binietoglou@0: # Laser pointing angle binietoglou@0: temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',)) binietoglou@0: temp_v[:] = params.general_parameters['Laser_Pointing_Angle'] binietoglou@0: binietoglou@0: # Molecular calculation binietoglou@0: temp_v = f.createVariable('Molecular_Calc','i') binietoglou@0: temp_v[:] = params.general_parameters['Molecular_Calc'] binietoglou@0: binietoglou@0: # Laser pointing angles of profiles binietoglou@0: temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales')) binietoglou@0: for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'], binietoglou@0: range(len(self.variables['Raw_Data_Start_Time']))): binietoglou@0: temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle binietoglou@0: binietoglou@0: # Raw data start/stop time binietoglou@0: temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales')) binietoglou@0: temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales')) binietoglou@0: for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'], binietoglou@0: self.variables['Raw_Data_Stop_Time'], binietoglou@0: range(len(self.variables['Raw_Data_Start_Time']))): binietoglou@0: temp_raw_start[:len(start_time),n] = start_time binietoglou@0: temp_raw_stop[:len(stop_time),n] = stop_time binietoglou@0: binietoglou@0: #Laser shots binietoglou@0: temp_v = f.createVariable('Laser_Shots','i',('time','channels')) binietoglou@0: for (channel,n) in zip(channels, range(len(channels))): binietoglou@0: time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]]) binietoglou@0: temp_v[:time_length, n] = params.channel_parameters[channel]['Laser_Shots'] binietoglou@0: binietoglou@0: #Raw lidar data binietoglou@0: temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points')) binietoglou@0: for (channel,n) in zip(channels, range(len(channels))): binietoglou@0: c = self.channels[channel] binietoglou@0: temp_v[:len(c.time),n, :c.points] = c.matrix binietoglou@0: binietoglou@0: self.add_dark_measurements_to_netcdf(f, channels) binietoglou@0: binietoglou@0: #Pressure at lidar station binietoglou@0: temp_v = f.createVariable('Pressure_at_Lidar_Station','d') binietoglou@0: temp_v[:] = self.info['Pressure'] binietoglou@0: binietoglou@0: #Temperature at lidar station binietoglou@0: temp_v = f.createVariable('Temperature_at_Lidar_Station','d') binietoglou@0: temp_v[:] = self.info['Temperature'] binietoglou@0: binietoglou@0: self.save_netcdf_extra(f) binietoglou@0: f.close() binietoglou@0: binietoglou@0: def add_dark_measurements_to_netcdf(self, f, channels): binietoglou@0: binietoglou@0: # Get dark measurements. If it is not given in self.dark_measurement binietoglou@0: # try to get it using the get_dark_measurements method. If none is found binietoglou@0: # return without adding something. binietoglou@0: if self.dark_measurement is None: binietoglou@0: self.dark_measurement = self.get_dark_measurements() binietoglou@0: binietoglou@0: if self.dark_measurement is None: binietoglou@0: return binietoglou@0: binietoglou@0: dark_measurement = self.dark_measurement binietoglou@0: binietoglou@0: # Calculate the length of the time_bck dimensions binietoglou@0: number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()] binietoglou@0: max_number_of_profiles = np.max(number_of_profiles) binietoglou@0: binietoglou@0: # Create the dimension binietoglou@0: f.createDimension('time_bck', max_number_of_profiles) binietoglou@0: binietoglou@0: # Save the dark measurement data binietoglou@0: temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points')) binietoglou@0: for (channel,n) in zip(channels, range(len(channels))): binietoglou@0: c = dark_measurement.channels[channel] binietoglou@0: temp_v[:len(c.time),n, :c.points] = c.matrix binietoglou@0: binietoglou@0: # Dark profile start/stop time binietoglou@0: temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time','nb_of_time_scales')) binietoglou@0: temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time','nb_of_time_scales')) binietoglou@0: for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'], binietoglou@0: dark_measurement.variables['Raw_Data_Stop_Time'], binietoglou@0: range(len(dark_measurement.variables['Raw_Data_Start_Time']))): binietoglou@0: temp_raw_start[:len(start_time),n] = start_time binietoglou@0: temp_raw_stop[:len(stop_time),n] = stop_time binietoglou@0: binietoglou@0: # Dark measurement start/stop time binietoglou@0: f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d') binietoglou@0: f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S') binietoglou@0: f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S') binietoglou@0: binietoglou@0: binietoglou@0: binietoglou@0: def save_netcdf_extra(self, f): binietoglou@0: pass binietoglou@0: binietoglou@0: def _gettime(self, date_str, time_str): binietoglou@0: t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S') binietoglou@0: return t binietoglou@0: binietoglou@0: def plot(self): binietoglou@0: for channel in self.channels: binietoglou@0: self.channels[channel].plot(show_plot = False) binietoglou@0: plt.show() binietoglou@0: binietoglou@0: def get_dark_measurements(self): binietoglou@0: return None binietoglou@0: binietoglou@0: binietoglou@0: class Lidar_channel: binietoglou@0: binietoglou@0: def __init__(self,channel_parameters): binietoglou@0: c = 299792458 #Speed of light binietoglou@0: self.wavelength = channel_parameters['name'] binietoglou@0: self.name = str(self.wavelength) binietoglou@0: self.binwidth = float(channel_parameters['binwidth']) # in microseconds binietoglou@0: self.data = {} binietoglou@0: self.resolution = self.binwidth * c / 2 binietoglou@0: self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution/2.0 # Change: add half bin in the z binietoglou@0: self.points = len(channel_parameters['data']) binietoglou@0: self.rc = [] binietoglou@0: self.duration = 60 binietoglou@0: binietoglou@0: def calculate_rc(self): binietoglou@0: background = np.mean(self.matrix[:,4000:], axis = 1) #Calculate the background from 30000m and above binietoglou@0: self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2) binietoglou@0: binietoglou@0: binietoglou@0: def update(self): binietoglou@0: self.start_time = min(self.data.keys()) binietoglou@0: self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration) binietoglou@0: self.time = tuple(sorted(self.data.keys())) binietoglou@0: sorted_data = sorted(self.data.iteritems(), key=itemgetter(0)) binietoglou@0: self.matrix = np.array(map(itemgetter(1),sorted_data)) binietoglou@0: binietoglou@0: def _nearest_dt(self,dtime): binietoglou@0: margin = datetime.timedelta(seconds = 300) binietoglou@0: if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time): binietoglou@0: print "Requested date not covered in this file" binietoglou@0: raise binietoglou@0: dt = abs(self.time - np.array(dtime)) binietoglou@0: dtmin = min(dt) binietoglou@0: binietoglou@0: if dtmin > datetime.timedelta(seconds = 60): binietoglou@0: print "Nearest profile more than 60 seconds away. dt = %s." % dtmin binietoglou@0: ind_t = np.where(dt == dtmin) binietoglou@0: ind_a= ind_t[0] binietoglou@0: if len(ind_a) > 1: binietoglou@0: ind_a = ind_a[0] binietoglou@0: chosen_time = self.time[ind_a] binietoglou@0: return chosen_time, ind_a binietoglou@0: binietoglou@0: def subset_by_time(self, start_time, stop_time): binietoglou@0: binietoglou@0: time_array = np.array(self.time) binietoglou@0: condition = (time_array >= start_time) & (time_array <= stop_time) binietoglou@0: binietoglou@0: subset_time = time_array[condition] binietoglou@0: subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time]) binietoglou@0: binietoglou@0: #Create a list with the values needed by channel's __init__() binietoglou@0: parameters_values = {'name': self.wavelength, binietoglou@0: 'binwidth': self.binwidth, binietoglou@0: 'data': subset_data[subset_time[0]],} binietoglou@0: binietoglou@0: c = Lidar_channel(parameters_values) binietoglou@0: c.data = subset_data binietoglou@0: c.update() binietoglou@0: return c binietoglou@0: binietoglou@0: binietoglou@0: def profile(self,dtime, signal_type = 'rc'): binietoglou@0: t, idx = self._nearest_dt(dtime) binietoglou@0: if signal_type == 'rc': binietoglou@0: data = self.rc binietoglou@0: else: binietoglou@0: data = self.matrix binietoglou@0: binietoglou@0: prof = data[idx,:][0] binietoglou@0: return prof, t binietoglou@0: binietoglou@0: def get_slice(self, starttime, endtime, signal_type = 'rc'): binietoglou@0: if signal_type == 'rc': binietoglou@0: data = self.rc binietoglou@0: else: binietoglou@0: data = self.matrix binietoglou@0: tim = np.array(self.time) binietoglou@0: starttime = self._nearest_dt(starttime)[0] binietoglou@0: endtime = self._nearest_dt(endtime)[0] binietoglou@0: condition = (tim >= starttime) & (tim <= endtime) binietoglou@0: sl = data[condition, :] binietoglou@0: t = tim[condition] binietoglou@0: return sl,t binietoglou@0: binietoglou@0: def av_profile(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'): binietoglou@0: starttime = tim - duration/2 binietoglou@0: endtime = tim + duration/2 binietoglou@0: d,t = self.get_slice(starttime, endtime, signal_type = signal_type) binietoglou@0: prof = np.mean(d, axis = 0) binietoglou@0: tmin = min(t) binietoglou@0: tmax = max(t) binietoglou@0: tav = tmin + (tmax-tmin)/2 binietoglou@0: return prof,(tav, tmin,tmax) binietoglou@0: binietoglou@0: def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet): binietoglou@0: #if filename is not None: binietoglou@0: # matplotlib.use('Agg') binietoglou@0: binietoglou@0: fig = plt.figure() binietoglou@0: ax1 = fig.add_subplot(111) binietoglou@0: self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom) binietoglou@0: ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name)) binietoglou@0: binietoglou@0: if filename is not None: binietoglou@0: pass binietoglou@0: #plt.savefig(filename) binietoglou@0: else: binietoglou@0: if show_plot: binietoglou@0: plt.show() binietoglou@0: #plt.close() ??? binietoglou@0: binietoglou@0: def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', zoom = [0,12000,0,-1]): binietoglou@0: binietoglou@0: if signal_type == 'rc': binietoglou@0: if len(self.rc) == 0: binietoglou@0: self.calculate_rc() binietoglou@0: data = self.rc binietoglou@0: else: binietoglou@0: data = self.matrix binietoglou@0: binietoglou@0: hmax_idx = self.index_at_height(zoom[1]) binietoglou@0: binietoglou@0: ax1.set_ylabel('Altitude (km)') binietoglou@0: ax1.set_xlabel('Time UTC') binietoglou@0: #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 binietoglou@0: #dateFormatter = mpl.dates.DateFormatter('%H.%M') binietoglou@0: #hourlocator = mpl.dates.HourLocator() binietoglou@0: binietoglou@0: #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m') binietoglou@0: #daylocator = mpl.dates.DayLocator() binietoglou@0: hourFormatter = mpl.dates.DateFormatter('%H.%M') binietoglou@0: hourlocator = mpl.dates.AutoDateLocator(interval_multiples=True) binietoglou@0: binietoglou@0: binietoglou@0: #ax1.axes.xaxis.set_major_formatter(dayFormatter) binietoglou@0: #ax1.axes.xaxis.set_major_locator(daylocator) binietoglou@0: ax1.axes.xaxis.set_major_formatter(hourFormatter) binietoglou@0: ax1.axes.xaxis.set_major_locator(hourlocator) binietoglou@0: binietoglou@0: binietoglou@0: ts1 = mpl.dates.date2num(self.start_time) binietoglou@0: ts2 = mpl.dates.date2num(self.stop_time) binietoglou@0: binietoglou@0: binietoglou@0: im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]], binietoglou@0: aspect = 'auto', binietoglou@0: origin = 'lower', binietoglou@0: cmap = cmap, binietoglou@0: #vmin = 0, binietoglou@0: vmin = data[:,10:400].max() * 0.1, binietoglou@0: #vmax = 1.4*10**7, binietoglou@0: vmax = data[:,10:400].max() * 0.9, binietoglou@0: extent = [ts1,ts2,self.z[zoom[0]]/1000.0, self.z[hmax_idx]/1000.0], binietoglou@0: ) binietoglou@0: binietoglou@0: cb1 = plt.colorbar(im1) binietoglou@0: cb1.ax.set_ylabel('a.u.') binietoglou@0: binietoglou@0: def index_at_height(self, height): binietoglou@0: idx = np.array(np.abs(self.z - height).argmin()) binietoglou@0: if idx.size >1: binietoglou@0: idx =idx[0] binietoglou@0: return idx binietoglou@0: binietoglou@0: def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID): binietoglou@0: #Read the lidar files and select channels binietoglou@0: temp_m = LidarClass(files) binietoglou@0: m = temp_m.subset_by_channels(channels) binietoglou@0: m.get_PT() binietoglou@0: m.info['Measurement_ID'] = measurement_ID binietoglou@0: m.save_as_netcdf(filename) binietoglou@0: