Mercurial > public > atmospheric_lidar / changeset

--- a/atmospheric_lidar/generic.py	Wed Dec 12 18:02:39 2018 +0200
+++ b/atmospheric_lidar/generic.py	Fri Dec 14 16:11:54 2018 +0200
@@ -712,7 +712,7 @@
         duration = np.array(durations)
         return duration
 
-    def calculate_rc(self, idx_min=4000, idx_max=5000):
+    def calculate_rc(self, idx_min=-2000, idx_max=-500):
         """ Calculate range corrected signal.
         
         The background signal is estimated as the mean signal between idx_min and idx_max. 
@@ -726,9 +726,33 @@
         idx_max : int
            Maximum index to calculate background signal.
         """
-        background = np.mean(self.matrix[:, idx_min:idx_max], axis=1)  # Calculate the background from 30000m and above
+        background = np.mean(self.matrix[:, idx_min:idx_max], axis=1)
         self.rc = (self.matrix.transpose() - background).transpose() * (self.z ** 2)
 
+    def noise_mask(self, idx_min=-2000, idx_max=-500, threshold=1.):
+        """ Calculate points that are probably noise.
+
+        To calculate this, we find the max value of the background region. Then we reject all points that
+        are less than `threshold` times this value.
+        Parameters
+        ----------
+        idx_min : int
+           Minimum index to calculate background signal.
+        idx_max : int
+           Maximum index to calculate background signal.
+        threshold : float
+           Threshold value.
+        """
+        background_max = np.max(self.matrix[:, idx_min:idx_max], axis=1)
+        background_mean = np.mean(self.matrix[:, idx_min:idx_max], axis=1)
+        mean_to_max = background_max - background_mean
+
+        noise_limit = background_mean + mean_to_max * threshold
+
+        mask = self.matrix <= noise_limit[:, np.newaxis]
+
+        return mask
+
     def update(self):
         """
         Update the time parameters and data according to the raw input data. 

--- a/atmospheric_lidar/raymetrics.py	Wed Dec 12 18:02:39 2018 +0200
+++ b/atmospheric_lidar/raymetrics.py	Fri Dec 14 16:11:54 2018 +0200
@@ -1,6 +1,10 @@
 """ Code to read Raymetrics version of Licel binary files."""
 import logging
 
+import numpy as np
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+
 from .licel import LicelFile, LicelLidarMeasurement, LicelChannel, PhotodiodeChannel
 
 logger = logging.getLogger(__name__)
@@ -153,6 +157,151 @@
         self._assign_unique_property('zenith_stop', current_file.zenith_stop)
         self._assign_unique_property('zenith_step', current_file.zenith_step)
 
+    def plot_ppi(self, figsize=(8, 4), signal_type='rc', z_min=0., z_max=12000., show_plot=True,
+             cmap=plt.cm.jet, title=None, vmin=0, vmax=1.3 * 10 ** 7, mask_noise=True, noise_threshold=1.):
+        """
+        Plot a vertical project of channel data.
+
+        Parameters
+        ----------
+        figsize : tuple
+           (width, height) of the output figure (inches)
+        signal_type : str
+           If 'rc', the range corrected signal is ploted. Else, the raw signals are used.
+        z_min : float
+           Minimum z range
+        z_max : float
+           Maximum z range
+        show_plot : bool
+           If True, the show_plot command is run.
+        cmap : cmap
+           An instance of a matplotlib colormap to use.
+        z0 : float
+           The ground-level altitude. If provided the plot shows altitude above sea level.
+        title : str
+           Optional title for the plot.
+        vmin : float
+           Minimum value for the color scale.
+        vmax : float
+           Maximum value for the color scale.
+        mask_noise : bool
+           If True, remove noisy bins.
+        noise_threshold : int
+           Threshold to use in the noise masking routine.
+        """
+        fig = plt.figure(figsize=figsize)
+        ax1 = fig.add_subplot(111)
+
+        self.draw_ppi(ax1, cmap=cmap, signal_type=signal_type, z_min=z_min, z_max=z_max, vmin=vmin, vmax=vmax,
+                      mask_noise=mask_noise, noise_threshold=noise_threshold)
+
+        if title:
+            ax1.set_title(title)
+        else:
+            ax1.set_title("PPI scan")
+
+        if show_plot:
+            plt.show()
+
+    def draw_ppi(self, ax1, cmap=plt.cm.jet, signal_type='rc',
+                  z_min=0, z_max=12000., add_colorbar=True, cmap_label='a.u.', cb_format=None,
+                  vmin=0, vmax=1.3 * 10 ** 7, mask_noise=True, noise_threshold=1.):
+        """
+        Draw channel data as a PPI plot.
+
+        Parameters
+        ----------
+        ax1 : axis object
+           The axis object to draw.
+        cmap : cmap
+           An instance of a matplotlib colormap to use.
+        signal_type : str
+           If 'rc', the range corrected signal is ploted. Else, the raw signals are used.
+        z_min : float
+           Minimum z range
+        z_max : float
+           Maximum z range
+        add_colorbar : bool
+           If True, a colorbar will be added to the plot.
+        cmap_label : str
+           Label for the colorbar. Ignored if add_colorbar is False.
+        cb_format : str
+           Colorbar tick format string.
+        vmin : float
+           Minimum value for the color scale.
+        vmax : float
+           Maximum value for the color scale.
+        mask_noise : bool
+           If True, remove noisy bins.
+        noise_threshold : int
+           Threshold to use in the noise masking routine.
+        """
+        if signal_type == 'rc':
+            if len(self.rc) == 0:
+                self.calculate_rc()
+            data = self.rc
+        else:
+            data = self.matrix
+
+        if mask_noise:
+            mask = self.noise_mask(threshold=noise_threshold)
+            data = np.ma.masked_where(mask, data)
+
+        z_min_idx = self._index_at_height(z_min)
+        z_max_idx = self._index_at_height(z_max)
+
+        x, y = self._polar_to_ground(self.z[z_min_idx:z_max_idx], self.azimuth_angles, self.zenith_angles)
+
+        im1 = ax1.pcolormesh(x, y, data[:, z_min_idx:z_max_idx], cmap=cmap, vmin=vmin, vmax=vmax)
+
+        ax1.axis('equal')
+
+        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)
+
+    @staticmethod
+    def _polar_to_ground(z, azimuth, zenith):
+        """
+        Convert polar coordinates to cartesian project for a PPI scan
+
+        Parameters
+        ----------
+        z : array
+           Distance array in meters
+        azimuth : list
+           List of profile azimuth angles in degrees
+        zenith : list
+           List of profile zenith angles in degrees
+
+        Returns
+        -------
+        x : array
+           X axis in meters
+        y : array
+           Y axis in meters
+        """
+        # Generate the mesh
+        zenith_rad = np.deg2rad(zenith)
+        azimuth_rad = np.deg2rad(azimuth)
+
+        Z, Zeniths = np.meshgrid(z, zenith_rad)
+        Z_ground = Z * np.sin(Zeniths)
+
+        x = Z_ground * np.cos(azimuth_rad)[:, np.newaxis]
+        y = Z_ground * np.sin(azimuth_rad)[:, np.newaxis]
+
+        return x, y
+
 
 class ScanningLidarMeasurement(LicelLidarMeasurement):
     """ A class representing a scanning measurement set.