Mercurial > public > atmospheric_lidar / changeset
changeset
Merge from 27:6f6f2c512a1f
Sun, 23 Nov 2014 23:25:09 +0200
- author
- Iannis <ulalume3@yahoo.com>
- date
- Sun, 23 Nov 2014 23:25:09 +0200
- changeset 28
- 73337ce10473
- parent 27
- 74f7617f5356 (diff)
- parent 26
- 84f7d763deeb (current diff)
- child 30
- 28d7b0974fe6
Merge from 27:6f6f2c512a1f
| licel.py | file | annotate | diff | comparison | revisions |
--- a/generic.py Tue Jul 01 11:31:37 2014 +0300 +++ b/generic.py Sun Nov 23 23:25:09 2014 +0200 @@ -5,6 +5,7 @@ # 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 @@ -151,6 +152,20 @@ 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. @@ -390,13 +405,13 @@ return t_mean -class Lidar_channel: +class LidarChannel: - def __init__(self,channel_parameters): - c = 299792458 #Speed of light + 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.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 @@ -405,7 +420,7 @@ 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) @@ -442,16 +457,39 @@ 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]],} + + # 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 = Lidar_channel(parameters_values) + 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': @@ -512,7 +550,10 @@ 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: @@ -561,19 +602,101 @@ #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 cb_format: - cb1 = plt.colorbar(im1, format = cb_format) + if signal_type == 'rc': + if len(self.rc) == 0: + self.calculate_rc() + data = self.rc else: - cb1 = plt.colorbar(im1) - cb1.ax.set_ylabel(cmap_label) + 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() + - # 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) + 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:
--- a/licel.py Tue Jul 01 11:31:37 2014 +0300 +++ b/licel.py Sun Nov 23 23:25:09 2014 +0200 @@ -1,6 +1,6 @@ import numpy as np import datetime -from generic import BaseLidarMeasurement, Lidar_channel +from generic import BaseLidarMeasurement, LidarChannel import musa_netcdf_parameters import musa_2009_netcdf_parameters @@ -11,6 +11,7 @@ class LicelFile: + def __init__(self, filename): self.start_time = None self.stop_time = None @@ -61,7 +62,7 @@ if (a[0] != 13) | (b[0] != 10): print "Warning: No end of line found after record. File could be corrupt" - channel = LicelFileChannel(current_channel_info, raw_data) + channel = LicelFileChannel(current_channel_info, raw_data, self.duration()) channel_name = channel.channel_name if channel_name in channels.keys(): @@ -73,14 +74,20 @@ self.raw_info = raw_info self.channels = channels - - + + def duration(self): + """ Return the duration of the file. """ + dt = self.stop_time - self.start_time + return dt.seconds + + class LicelFileChannel: - def __init__(self, raw_info = None, raw_data = None): + def __init__(self, raw_info = None, raw_data = None, duration = None): self.raw_info = raw_info self.raw_data = raw_data - + self.duration = duration + @property def wavelength(self): if self.raw_info is not None: @@ -110,7 +117,7 @@ data = self.raw_data number_of_shots = float(self.raw_info['NShots']) - norm = data/number_of_shots + norm = data / number_of_shots dz = float(self.raw_info['BinW']) if self.raw_info['AnalogPhoton']=='0': @@ -128,7 +135,7 @@ # channel_data = norm*SR # CHANGE: # For the SCC the data are needed in photons - channel_data = norm*number_of_shots + channel_data = norm *number_of_shots #print res,c,cdata,norm #Calculate Z @@ -143,9 +150,9 @@ ''' ''' - extra_netcdf_parameters = musa_netcdf_parameters - raw_info = {} # Keep the raw info from the files + raw_info = {} # Keep the raw info from the files + durations = {} # Keep the duration of the files def import_file(self, filename): if filename in self.files: @@ -153,10 +160,11 @@ else: current_file = LicelFile(filename) self.raw_info[current_file.filename] = current_file.raw_info + self.durations[current_file.filename] = current_file.duration() for channel_name, channel in current_file.channels.items(): if channel_name not in self.channels: - self.channels[channel_name] = Licel_channel(channel) + self.channels[channel_name] = LicelChannel(channel) self.channels[channel_name].data[current_file.start_time] = channel.data self.files.append(current_file.filename) @@ -172,16 +180,10 @@ ''' Return the duration for a given time scale. If only a single file is imported, then this cannot be guessed from the time difference and the raw_info of the file are checked. - ''' if len(raw_start_in_seconds) == 1: # If only one file imported - raw_info = self.raw_info.itervalues().next() # Get the first (and only) raw_info - start_str = "%s %s" % (raw_info['StartDate'], raw_info['StartTime']) - end_str = "%s %s" % (raw_info['EndDate'], raw_info['EndTime']) - start_time = datetime.datetime.strptime(start_str, "%d/%m/%Y %H:%M:%S") - end_time = datetime.datetime.strptime(end_str, "%d/%m/%Y %H:%M:%S") - duration = end_time - start_time + duration = self.durations.itervalues().next() # Get the first (and only) raw_info duration_sec = duration.seconds else: duration_sec = np.diff(raw_start_in_seconds)[0] @@ -189,19 +191,19 @@ return duration_sec -class Licel_channel(Lidar_channel): +class LicelChannel(LidarChannel): def __init__(self, channel_file): - c = 299792458.0 #Speed of light + c = 299792458.0 # Speed of light self.wavelength = channel_file.wavelength self.name = channel_file.channel_name - self.binwidth = channel_file.dz * 2 / c # in microseconds + self.binwidth = channel_file.dz * 2 / c # in seconds self.data = {} self.resolution = channel_file.dz - self.z = np.arange(channel_file.number_of_bins) *self.resolution + self.resolution/2.0 #Change: add half bin in the z + self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z self.points = channel_file.number_of_bins self.rc = [] - self.duration = 60 + self.duration = channel_file.duration def append(self, other): if self.info != other.info:
--- a/pearl.py Tue Jul 01 11:31:37 2014 +0300 +++ b/pearl.py Sun Nov 23 23:25:09 2014 +0200 @@ -4,7 +4,7 @@ import numpy as np -from generic import BaseLidarMeasurement, Lidar_channel +from generic import BaseLidarMeasurement, LidarChannel from ciao import CiaoMixin import pearl_netcdf_parameters @@ -39,7 +39,7 @@ name = channel_info['name'] tm = start_time if name not in self.channels: - self.channels[name] = Lidar_channel(channel_info) + self.channels[name] = LidarChannel(channel_info) self.channels[name].data[tm] = channel_info['data'] self.files.append(filename)
