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: