--- a/generic.py Mon Feb 22 16:50:03 2016 +0200
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,714 +0,0 @@
-# General imports
-import datetime
-from operator import itemgetter
-
-# Science imports
-import numpy as np
-import matplotlib as mpl
-from matplotlib.ticker import ScalarFormatter
-from matplotlib import pyplot as plt
-import netCDF4 as netcdf
-
-netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4'
-
-
-class BaseLidarMeasurement():
- """ This is the general measurement object.
- It is meant to become a general measurement object
- independent of the input files.
-
- Each subclass should implement the following:
- * the import_file method.
- * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters.
-
- You can override the get_PT method to define a custom procedure to get ground temperature and pressure.
- The one implemented by default is by using the MILOS meteorological station data.
-
- """
-
- def __init__(self, filelist = None):
- self.info = {}
- self.dimensions = {}
- self.variables = {}
- self.channels = {}
- self.attributes = {}
- self.files = []
- self.dark_measurement = None
-
- if filelist:
- self.import_files(filelist)
-
- def import_files(self, filelist):
- for f in filelist:
- self.import_file(f)
- self.update()
-
- def import_file(self,filename):
- raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.')
-
- def update(self):
- '''
- Update the the info, variables and dimensions of the lidar measurement based
- on the information found in the channels.
-
- Reading of the scan_angles parameter is not implemented.
- '''
-
- # Initialize
- start_time =[]
- stop_time = []
- points = []
- all_time_scales = []
- channel_name_list = []
-
- # Get the information from all the channels
- for channel_name, channel in self.channels.items():
- channel.update()
- start_time.append(channel.start_time)
- stop_time.append(channel.stop_time)
- points.append(channel.points)
- all_time_scales.append(channel.time)
- channel_name_list.append(channel_name)
-
- # Find the unique time scales used in the channels
- time_scales = set(all_time_scales)
-
- # Update the info dictionary
- self.info['start_time'] = min(start_time)
- self.info['stop_time'] = max(stop_time)
- self.info['duration'] = self.info['stop_time'] - self.info['start_time']
-
- # Update the dimensions dictionary
- self.dimensions['points'] = max(points)
- self.dimensions['channels'] = len(self.channels)
- # self.dimensions['scan angles'] = 1
- self.dimensions['nb_of_time_scales'] = len(time_scales)
-
- # Update the variables dictionary
- # Write time scales in seconds
- raw_Data_Start_Time = []
- raw_Data_Stop_Time = []
-
- for current_time_scale in list(time_scales):
- raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time
- raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds
- raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list
- # Check if this time scale has measurements every 30 or 60 seconds.
-
- duration = self._get_duration(raw_start_in_seconds)
-
- raw_stop_in_seconds = raw_start_in_seconds + duration
- raw_Data_Stop_Time.append(raw_stop_in_seconds)
-
- self.variables['Raw_Data_Start_Time'] = raw_Data_Start_Time
- self.variables['Raw_Data_Stop_Time'] = raw_Data_Stop_Time
-
- # Make a dictionary to match time scales and channels
- channel_timescales = []
- for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales):
- # The following lines are PEARL specific. The reason they are here is not clear.
- # if channel_name =='1064BLR':
- # channel_name = '1064'
- for (ts,n) in zip(time_scales, range(len(time_scales))):
- if current_time_scale == ts:
- channel_timescales.append([channel_name,n])
- self.variables['id_timescale'] = dict(channel_timescales)
-
- def _get_duration(self, raw_start_in_seconds):
- ''' Return the duration for a given time scale. In some files (e.g. Licel) this
- can be specified from the files themselves. In others this must be guessed.
-
- '''
- # The old method, kept here for reference
- #dt = np.mean(np.diff(raw_start_in_seconds))
- #for d in duration_list:
- # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds
- # duration = d
-
- duration = np.diff(raw_start_in_seconds)[0]
-
- return duration
-
- def subset_by_channels(self, channel_subset):
- ''' Get a measurement object containing only the channels with names
- contained in the channel_sublet list '''
-
- m = self.__class__() # Create an object of the same type as this one.
- m.channels = dict([(channel, self.channels[channel]) for channel
- in channel_subset])
- m.update()
- return m
-
- def subset_by_time(self, start_time, stop_time):
-
- if start_time > stop_time:
- raise ValueError('Stop time should be after start time')
-
- if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']):
- raise ValueError('The time interval specified is not part of the measurement')
-
- m = self.__class__() # Create an object of the same type as this one.
- for (channel_name, channel) in self.channels.items():
- m.channels[channel_name] = channel.subset_by_time(start_time, stop_time)
- m.update()
- return m
-
- def subset_by_bins(self, b_min = 0, b_max = None):
- """Remove some height bins from the file. This could be needed to
- remove aquisition artifacts at the start or the end of the files.
- """
-
- m = self.__class__() # Create an object of the same type as this one.
-
- for (channel_name, channel) in self.channels.items():
- m.channels[channel_name] = channel.subset_by_bins(b_min, b_max)
-
- m.update()
-
- return m
-
- def r_plot(self):
- #Make a basic plot of the data.
- #Should include some dictionary with params to make plot stable.
- pass
-
- def r_pdf(self):
- # Create a pdf report using a basic plot and meta-data.
- pass
-
- def save(self):
- #Save the current state of the object to continue the analysis later.
- pass
-
- def get_PT(self):
- ''' Sets the pressure and temperature at station level .
- The results are stored in the info dictionary.
- '''
-
- self.info['Temperature'] = 10.0
- self.info['Pressure'] = 930.0
-
- def subtract_dark(self):
-
- if not self.dark_measurement:
- raise IOError('No dark measurements have been imported yet.')
-
- for (channel_name, dark_channel) in self.dark_measurement.channels.iteritems():
- dark_profile = dark_channel.average_profile()
- channel = self.channels[channel_name]
-
- for measurement_time, data in channel.data.iteritems():
- channel.data[measurement_time] = data - dark_profile
-
- channel.update()
-
- def save_as_netcdf(self, filename = None):
- """Saves the measurement in the netcdf format as required by the SCC.
- Input: filename. If no filename is provided <measurement_id>.nc will
- be used.
- """
- params = self.extra_netcdf_parameters
- needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure']
-
- for parameter in needed_parameters:
- stored_value = self.info.get(parameter, None)
- if stored_value is None:
- raise ValueError('A value needs to be specified for %s' % parameter)
-
- if not filename:
- filename = "%s.nc" % self.info['Measurement_ID']
-
- dimensions = {'points': 1,
- 'channels': 1,
- 'time': None,
- 'nb_of_time_scales': 1,
- 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented
-
- global_att = {'Measurement_ID': None,
- 'RawData_Start_Date': None,
- 'RawData_Start_Time_UT': None,
- 'RawData_Stop_Time_UT': None,
- 'RawBck_Start_Date': None,
- 'RawBck_Start_Time_UT': None,
- 'RawBck_Stop_Time_UT': None,
- 'Sounding_File_Name': None,
- 'LR_File_Name': None,
- 'Overlap_File_Name': None,
- 'Location': None,
- 'System': None,
- 'Latitude_degrees_north': None,
- 'Longitude_degrees_east': None,
- 'Altitude_meter_asl': None}
-
- channel_variables = \
- {'channel_ID': (('channels', ), 'i'),
- 'Background_Low': (('channels', ), 'd'),
- 'Background_High': (('channels', ), 'd'),
- 'LR_Input': (('channels', ), 'i'),
- 'DAQ_Range': (('channels', ), 'd'),
- 'Depolarization_Factor': (('channels', ), 'd'), }
-
-
- channels = self.channels.keys()
-
- input_values = dict(self.dimensions, **self.variables)
-
- # Add some mandatory global attributes
- input_values['Measurement_ID'] = self.info['Measurement_ID']
- input_values['RawData_Start_Date'] = self.info['start_time'].strftime('%Y%m%d')
- input_values['RawData_Start_Time_UT'] = self.info['start_time'].strftime('%H%M%S')
- input_values['RawData_Stop_Time_UT'] = self.info['stop_time'].strftime('%H%M%S')
-
- # Add some optional global attributes
- input_values['System'] = params.general_parameters['System']
- input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north']
- input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east']
- input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl']
-
- # Open a netCDF4 file
- f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file
-
- # Create the dimensions in the file
- for (d,v) in dimensions.iteritems():
- v = input_values.pop(d, v)
- f.createDimension(d,v)
-
- # Create global attributes
- for (attrib,value) in global_att.iteritems():
- val = input_values.pop(attrib,value)
- if val:
- setattr(f, attrib, val)
-
- """ Variables """
- # Write the values of fixes channel parameters
- for (var,t) in channel_variables.iteritems():
- temp_v = f.createVariable(var,t[1],t[0])
- for (channel, n) in zip(channels, range(len(channels))):
- temp_v[n] = params.channel_parameters[channel][var]
-
- # Write the id_timescale values
- temp_id_timescale = f.createVariable('id_timescale','i',('channels',))
- for (channel, n) in zip(channels, range(len(channels))):
- temp_id_timescale[n] = self.variables['id_timescale'][channel]
-
- # Laser pointing angle
- temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',))
- temp_v[:] = params.general_parameters['Laser_Pointing_Angle']
-
- # Molecular calculation
- temp_v = f.createVariable('Molecular_Calc','i')
- temp_v[:] = params.general_parameters['Molecular_Calc']
-
- # Laser pointing angles of profiles
- temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales'))
- for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'],
- range(len(self.variables['Raw_Data_Start_Time']))):
- temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle
-
- # Raw data start/stop time
- temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales'))
- temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales'))
- for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'],
- self.variables['Raw_Data_Stop_Time'],
- range(len(self.variables['Raw_Data_Start_Time']))):
- temp_raw_start[:len(start_time),n] = start_time
- temp_raw_stop[:len(stop_time),n] = stop_time
-
- #Laser shots
- temp_v = f.createVariable('Laser_Shots','i',('time','channels'))
- for (channel,n) in zip(channels, range(len(channels))):
- time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]])
- # Array slicing stoped working as usual ex. temp_v[:10] = 100 does not work. ??? np.ones was added.
- temp_v[:time_length, n] = np.ones(time_length) * params.channel_parameters[channel]['Laser_Shots']
-
- #Raw lidar data
- temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points'))
- for (channel,n) in zip(channels, range(len(channels))):
- c = self.channels[channel]
- temp_v[:len(c.time),n, :c.points] = c.matrix
-
- self.add_dark_measurements_to_netcdf(f, channels)
-
- #Pressure at lidar station
- temp_v = f.createVariable('Pressure_at_Lidar_Station','d')
- temp_v[:] = self.info['Pressure']
-
- #Temperature at lidar station
- temp_v = f.createVariable('Temperature_at_Lidar_Station','d')
- temp_v[:] = self.info['Temperature']
-
- self.save_netcdf_extra(f)
- f.close()
-
- def add_dark_measurements_to_netcdf(self, f, channels):
-
- # Get dark measurements. If it is not given in self.dark_measurement
- # try to get it using the get_dark_measurements method. If none is found
- # return without adding something.
- if self.dark_measurement is None:
- self.dark_measurement = self.get_dark_measurements()
-
- if self.dark_measurement is None:
- return
-
- dark_measurement = self.dark_measurement
-
- # Calculate the length of the time_bck dimensions
- number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()]
- max_number_of_profiles = np.max(number_of_profiles)
-
- # Create the dimension
- f.createDimension('time_bck', max_number_of_profiles)
-
- # Save the dark measurement data
- temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points'))
- for (channel,n) in zip(channels, range(len(channels))):
- c = dark_measurement.channels[channel]
- temp_v[:len(c.time),n, :c.points] = c.matrix
-
- # Dark profile start/stop time
- temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time_bck','nb_of_time_scales'))
- temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time_bck','nb_of_time_scales'))
- for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'],
- dark_measurement.variables['Raw_Data_Stop_Time'],
- range(len(dark_measurement.variables['Raw_Data_Start_Time']))):
- temp_raw_start[:len(start_time),n] = start_time
- temp_raw_stop[:len(stop_time),n] = stop_time
-
- # Dark measurement start/stop time
- f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d')
- f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S')
- f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S')
-
-
-
- def save_netcdf_extra(self, f):
- pass
-
- def _gettime(self, date_str, time_str):
- t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S')
- return t
-
- def plot(self):
- for channel in self.channels:
- self.channels[channel].plot(show_plot = False)
- plt.show()
-
- def get_dark_measurements(self):
- return None
-
- @property
- def mean_time(self):
- start_time = self.info['start_time']
- stop_time = self.info['stop_time']
- dt = stop_time - start_time
- t_mean = start_time + dt / 2
- return t_mean
-
-
-class LidarChannel:
-
- def __init__(self, channel_parameters):
- c = 299792458 # Speed of light
- self.wavelength = channel_parameters['name']
- self.name = str(self.wavelength)
- self.binwidth = float(channel_parameters['binwidth']) # in microseconds
- self.data = {}
- self.resolution = self.binwidth * c / 2
- self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution / 2.0 # Change: add half bin in the z
- self.points = len(channel_parameters['data'])
- self.rc = []
- self.duration = 60
-
- def calculate_rc(self, idx_min = 4000, idx_max = 5000):
- background = np.mean(self.matrix[:,idx_min:idx_max], axis = 1) # Calculate the background from 30000m and above
- self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2)
-
-
- def update(self):
- self.start_time = min(self.data.keys())
- self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration)
- self.time = tuple(sorted(self.data.keys()))
- sorted_data = sorted(self.data.iteritems(), key=itemgetter(0))
- self.matrix = np.array(map(itemgetter(1),sorted_data))
-
- def _nearest_dt(self,dtime):
- margin = datetime.timedelta(seconds = 300)
- if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time):
- print "Requested date not covered in this file"
- raise
- dt = abs(self.time - np.array(dtime))
- dtmin = min(dt)
-
- if dtmin > datetime.timedelta(seconds = 60):
- print "Nearest profile more than 60 seconds away. dt = %s." % dtmin
- ind_t = np.where(dt == dtmin)
- ind_a= ind_t[0]
- if len(ind_a) > 1:
- ind_a = ind_a[0]
- chosen_time = self.time[ind_a]
- return chosen_time, ind_a
-
- def subset_by_time(self, start_time, stop_time):
-
- time_array = np.array(self.time)
- condition = (time_array >= start_time) & (time_array <= stop_time)
-
- subset_time = time_array[condition]
- subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time])
-
- #Create a list with the values needed by channel's __init__()
- parameter_values = {'name': self.wavelength,
- 'binwidth': self.binwidth,
- 'data': subset_data[subset_time[0]],}
-
- # We should use __class__ instead of class name, so that this works properly
- # with subclasses
- # Ex: c = self.__class__(parameters_values)
- # This does not currently work with Licel files though
- c = LidarChannel(parameter_values)
- c.data = subset_data
- c.update()
- return c
-
- def subset_by_bins(self, b_min = 0, b_max = None):
- """Remove some height bins from the file. This could be needed to
- remove aquisition artifacts at the start or the end of the files.
- """
-
- subset_data = {}
-
- for ctime, cdata in self.data.items():
- subset_data[ctime] = cdata[b_min:b_max]
-
- #Create a list with the values needed by channel's __init__()
- parameters_values = {'name': self.wavelength,
- 'binwidth': self.binwidth,
- 'data': subset_data[subset_data.keys()[0]],} # We just need any array with the correct length
-
- c = LidarChannel(parameters_values)
- c.data = subset_data
- c.update()
- return c
-
- def profile(self,dtime, signal_type = 'rc'):
- t, idx = self._nearest_dt(dtime)
- if signal_type == 'rc':
- data = self.rc
- else:
- data = self.matrix
-
- prof = data[idx,:][0]
- return prof, t
-
- def get_slice(self, starttime, endtime, signal_type = 'rc'):
- if signal_type == 'rc':
- data = self.rc
- else:
- data = self.matrix
- tim = np.array(self.time)
- starttime = self._nearest_dt(starttime)[0]
- endtime = self._nearest_dt(endtime)[0]
- condition = (tim >= starttime) & (tim <= endtime)
- sl = data[condition, :]
- t = tim[condition]
- return sl,t
-
- def profile_for_duration(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'):
- """ Calculates the profile around a specific time (tim). """
- starttime = tim - duration/2
- endtime = tim + duration/2
- d,t = self.get_slice(starttime, endtime, signal_type = signal_type)
- prof = np.mean(d, axis = 0)
- tmin = min(t)
- tmax = max(t)
- tav = tmin + (tmax-tmin)/2
- return prof,(tav, tmin,tmax)
-
- def average_profile(self):
- """ Return the average profile (NOT range corrected) for all the duration of the measurement. """
- prof = np.mean(self.matrix, axis = 0)
- return prof
-
- def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet, z0 = None, title = None, vmin = 0, vmax = 1.3 * 10 ** 7):
- #if filename is not None:
- # matplotlib.use('Agg')
-
- fig = plt.figure()
- ax1 = fig.add_subplot(111)
- self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom, z0 = z0, vmin = vmin, vmax = vmax)
-
- if title:
- ax1.set_title(title)
- else:
- ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name))
-
- if filename is not None:
- pass
- #plt.savefig(filename)
- else:
- if show_plot:
- plt.show()
- #plt.close() ???
-
- def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc',
- zoom = [0,12000,0,-1], z0 = None,
- add_colorbar = True, cmap_label = 'a.u.', cb_format = None,
- vmin = 0, vmax = 1.3 * 10 ** 7):
-
- if signal_type == 'rc':
- if len(self.rc) == 0:
- self.calculate_rc()
- data = self.rc
- else:
- data = self.matrix
-
- hmax_idx = self.index_at_height(zoom[1])
-
- # If z0 is given, then the plot is a.s.l.
- if z0:
- ax1.set_ylabel('Altitude a.s.l. [km]')
- else:
- ax1.set_ylabel('Altitude a.g.l. [km]')
- z0 = 0
-
- ax1.set_xlabel('Time UTC [hh:mm]')
- #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064
- #dateFormatter = mpl.dates.DateFormatter('%H.%M')
- #hourlocator = mpl.dates.HourLocator()
-
- #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m')
- #daylocator = mpl.dates.DayLocator()
- hourFormatter = mpl.dates.DateFormatter('%H:%M')
- hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True)
-
-
- #ax1.axes.xaxis.set_major_formatter(dayFormatter)
- #ax1.axes.xaxis.set_major_locator(daylocator)
- ax1.axes.xaxis.set_major_formatter(hourFormatter)
- ax1.axes.xaxis.set_major_locator(hourlocator)
-
-
- ts1 = mpl.dates.date2num(self.start_time)
- ts2 = mpl.dates.date2num(self.stop_time)
-
-
- im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]],
- aspect = 'auto',
- origin = 'lower',
- cmap = cmap,
- vmin = vmin,
- #vmin = data[:,10:400].max() * 0.1,
- vmax = vmax,
- #vmax = data[:,10:400].max() * 0.9,
- extent = [ts1,ts2,self.z[zoom[0]]/1000.0 + z0/1000., self.z[hmax_idx]/1000.0 + z0/1000.],
- )
-
- if add_colorbar:
- if cb_format:
- cb1 = plt.colorbar(im1, format = cb_format)
- else:
- cb1 = plt.colorbar(im1)
- cb1.ax.set_ylabel(cmap_label)
-
- # Make the ticks of the colorbar smaller, two points smaller than the default font size
- cb_font_size = mpl.rcParams['font.size'] - 2
- for ticklabels in cb1.ax.get_yticklabels():
- ticklabels.set_fontsize(cb_font_size)
- cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size)
-
-
- def draw_plot_new(self, ax1, cmap = plt.cm.jet, signal_type = 'rc',
- zoom = [0, 12000, 0, None], z0 = None,
- add_colorbar = True, cmap_label = 'a.u.',
- cb_format = None, power_limits = (-2, 2),
- date_labels = False,
- vmin = 0, vmax = 1.3 * 10 ** 7):
-
- if signal_type == 'rc':
- if len(self.rc) == 0:
- self.calculate_rc()
- data = self.rc
- else:
- data = self.matrix
-
- hmax_idx = self.index_at_height(zoom[1])
- hmin_idx = self.index_at_height(zoom[0])
-
- # If z0 is given, then the plot is a.s.l.
- if z0:
- ax1.set_ylabel('Altitude a.s.l. [km]')
- else:
- ax1.set_ylabel('Altitude a.g.l. [km]')
- z0 = 0
-
- ax1.set_xlabel('Time UTC [hh:mm]')
- #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064
- #dateFormatter = mpl.dates.DateFormatter('%H.%M')
- #hourlocator = mpl.dates.HourLocator()
-
-
- if date_labels:
- dayFormatter = mpl.dates.DateFormatter('%H:%M\n%d/%m/%Y')
- daylocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True)
- ax1.axes.xaxis.set_major_formatter(dayFormatter)
- ax1.axes.xaxis.set_major_locator(daylocator)
- else:
- hourFormatter = mpl.dates.DateFormatter('%H:%M')
- hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True)
- ax1.axes.xaxis.set_major_formatter(hourFormatter)
- ax1.axes.xaxis.set_major_locator(hourlocator)
-
-
- # Get the values of the time axis
- dt = datetime.timedelta(seconds = self.duration)
-
- time_cut = self.time[zoom[2]:zoom[3]]
- time_last = time_cut[-1] + dt # The last element needed for pcolormesh
-
- time_all = time_cut + (time_last,)
-
- t_axis = mpl.dates.date2num(time_all)
-
- # Get the values of the z axis
- z_cut = self.z[hmin_idx:hmax_idx] - self.resolution / 2.
- z_last = z_cut[-1] + self.resolution
-
- z_axis = np.append(z_cut, z_last) / 1000. + z0 / 1000. # Convert to km
-
- # Plot
- im1 = ax1.pcolormesh(t_axis, z_axis, data.T[hmin_idx:hmax_idx, zoom[2]:zoom[3]],
- cmap = cmap,
- vmin = vmin,
- vmax = vmax,
- )
-
- if add_colorbar:
- if cb_format:
- cb1 = plt.colorbar(im1, format = cb_format)
- else:
- cb_formatter = ScalarFormatter()
- cb_formatter.set_powerlimits(power_limits)
- cb1 = plt.colorbar(im1, format = cb_formatter)
- cb1.ax.set_ylabel(cmap_label)
-
- # Make the ticks of the colorbar smaller, two points smaller than the default font size
- cb_font_size = mpl.rcParams['font.size'] - 2
- for ticklabels in cb1.ax.get_yticklabels():
- ticklabels.set_fontsize(cb_font_size)
- cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size)
-
- def index_at_height(self, height):
- idx = np.array(np.abs(self.z - height).argmin())
- if idx.size > 1:
- idx =idx[0]
- return idx
-
-def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID):
- #Read the lidar files and select channels
- temp_m = LidarClass(files)
- m = temp_m.subset_by_channels(channels)
- m.get_PT()
- m.info['Measurement_ID'] = measurement_ID
- m.save_as_netcdf(filename)
-
