Atmospheric Lidar Package: comparison generic.py

-:870fc8f65eeb
+:28d7b0974fe6
 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'
 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.
 dt = stop_time - start_time
 t_mean = start_time + dt / 2
 return t_mean
-class Lidar_channel:
+class LidarChannel:
-def __init__(self,channel_parameters):
+def __init__(self, channel_parameters):
-c = 299792458 #Speed of light
+c = 299792458  # Speed of light
 self.wavelength = channel_parameters['name']
 self.name = str(self.wavelength)
-self.binwidth = float(channel_parameters['binwidth']) # in microseconds
+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):
-background = np.mean(self.matrix[:,4000:], axis = 1) #Calculate the background from 30000m and above
+background = np.mean(self.matrix[:,4000:], 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())
 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__()
-parameters_values = {'name': self.wavelength,
+parameter_values = {'name': self.wavelength,
 'binwidth': self.binwidth,
 'data': subset_data[subset_time[0]],}
-c = Lidar_channel(parameters_values)
+# 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:
 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, cmap_label = 'a.u.',  cb_format = None, vmin = 0, vmax = 1.3 * 10 ** 7):
+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
 #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 cb_format:
+if add_colorbar:
-cb1 = plt.colorbar(im1, format = cb_format)
+if cb_format:
-else:
+cb1 = plt.colorbar(im1, format = cb_format)
-cb1 = plt.colorbar(im1)
+else:
-cb1.ax.set_ylabel(cmap_label)
+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
+# Make the ticks of the colorbar smaller, two points smaller than the default font size
-for ticklabels in cb1.ax.get_yticklabels():
+cb_font_size = mpl.rcParams['font.size'] - 2
-ticklabels.set_fontsize(cb_font_size)
+for ticklabels in cb1.ax.get_yticklabels():
-cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size)
+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

Mercurial > public > atmospheric_lidar / file comparison

comparison: generic.py

generic.py