New plots for alex.

Sun, 23 Nov 2014 23:24:32 +0200

author
Iannis <ulalume3@yahoo.com>
date
Sun, 23 Nov 2014 23:24:32 +0200
changeset 27
74f7617f5356
parent 25
563705a26f85
child 28
73337ce10473

New plots for alex.

generic.py file | annotate | diff | comparison | revisions
licel.py file | annotate | diff | comparison | revisions
pearl.py file | annotate | diff | comparison | revisions
--- a/generic.py	Wed Jun 18 17:56:46 2014 +0300
+++ b/generic.py	Sun Nov 23 23:24:32 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	Wed Jun 18 17:56:46 2014 +0300
+++ b/licel.py	Sun Nov 23 23:24:32 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	Wed Jun 18 17:56:46 2014 +0300
+++ b/pearl.py	Sun Nov 23 23:24:32 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)
     

mercurial