Atmospheric Lidar Package: comparison atmospheric

-:34866a2a1aa5
+:580821c7487c
 return raw_info
 def _assign_properties(self):
 """ Assign scanning-specific parameters found in the header as object properties."""
 super(ScanningFile, self)._assign_properties()
-self.azimuth_angle = float(self.raw_info['azimuth_angle'])
+self.azimuth_angle_raw = float(self.raw_info['azimuth_angle'])
 self.temperature = float(self.raw_info['temperature'])
 self.pressure = float(self.raw_info['pressure'])
 self.azimuth_start_raw = float(self.raw_info['azimuth_start'])
 self.azimuth_stop_raw = float(self.raw_info['azimuth_stop'])
 self.azimuth_step = float(self.raw_info['azimuth_step'])
 self.zenith_start_raw = float(self.raw_info['zenith_start'])
 self.zenith_stop_raw = float(self.raw_info['zenith_stop'])
 self.zenith_step = float(self.raw_info['zenith_step'])
 self.azimuth_offset = float(self.raw_info['azimuth_offset'])
-logger.warning("Azimuth offset correction not applied.")
+self.azimuth_angle = (self.azimuth_angle_raw + self.azimuth_offset) % 360
-# TODO: Apply azimuth offset correction.
+self.azimuth_start = (self.azimuth_start_raw + self.azimuth_offset) % 360
-self.azimuth_start = self.azimuth_start_raw
+self.azimuth_stop = (self.azimuth_stop_raw + self.azimuth_offset) % 360
-self.azimuth_stop = self.azimuth_stop_raw
 if self.fix_zenith_angle:
 logger.debug('Fixing zenith start and zenith stop angles.')
 self.zenith_start = self._correct_zenith_angle(self.zenith_start_raw)
 self.zenith_stop = self._correct_zenith_angle(self.zenith_stop_raw)
 self.zenith_stop = None
 self.zenith_step = None
 self.azimuth_offset = None
 self.zenith_angles = []
 self.azimuth_angles = []
+self.temperature = []
+self.pressure = []
 def append_file(self, current_file, file_channel):
 """ Keep track of scanning-specific variable properties of each file. """
 super(ScanningChannel, self).append_file(current_file, file_channel)
 self.zenith_angles.append(current_file.zenith_angle)
 self.azimuth_angles.append(current_file.azimuth_angle)
+self.temperature.append(current_file.temperature)
+self.pressure.append(current_file.pressure)
 def _assign_properties(self, current_file, file_channel):
 """ Assign scanning-specific properties as object properties. Check that these are unique,
 i.e. that all files belong to the same measurements set.
 A specific LicelChannelData object holding data found in the file.
 """
 super(ScanningChannel, self)._assign_properties(current_file, file_channel)
 self._assign_unique_property('azimuth_start', current_file.azimuth_start)
 self._assign_unique_property('azimuth_stop', current_file.azimuth_stop)
+self._assign_unique_property('azimuth_start_raw', current_file.azimuth_start_raw)
+self._assign_unique_property('azimuth_stop_raw', current_file.azimuth_stop_raw)
 self._assign_unique_property('azimuth_step', current_file.azimuth_step)
 self._assign_unique_property('zenith_start', current_file.zenith_start)
 self._assign_unique_property('zenith_stop', current_file.zenith_stop)
+self._assign_unique_property('zenith_start_raw', current_file.zenith_start_raw)
+self._assign_unique_property('zenith_stop_raw', current_file.zenith_stop_raw)
 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.):
 """
 Threshold to use in the noise masking routine.
 """
 fig = plt.figure(figsize=figsize)
 ax1 = fig.add_subplot(111)
-self.draw_rhi(ax1, cmap=cmap, signal_type=signal_type, z_min=z_min, z_max=z_max, vmin=vmin, vmax=vmax,
+projection_angle = self.draw_rhi(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("RHI scan")
+ax1.set_title("RHI scan ({0}$^\circ$)".format(projection_angle))
 if show_plot:
 plt.show()
+def plot_scan(self, figsize=(8, 4), signal_type='rc', z_min=0., z_max=12000., show_plot=True,
+cmap=plt.cm.jet, vmin=0, vmax=1.3 * 10 ** 7, mask_noise=True, noise_threshold=1., cb_format='%.0e', box=False):
+"""
+Plot data as RHI and PPI scans.
+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(121)
+ax2 = fig.add_subplot(122)
+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, add_colorbar=False, cb_format=cb_format,
+box=box)
+projection_angle = self.draw_rhi(ax2, 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, cb_format=cb_format, box=box)
+fig.suptitle("Channel {0}: {1} - {2}".format(self.name,
+self.start_time.strftime('%Y%m%dT%H%M'),
+self.stop_time.strftime('%Y%m%dT%H%M')))
+ax1.set_title('PPI')
+ax2.set_title("RHI ({0:.1f}$^\circ$)".format(projection_angle))
+plt.tight_layout()
+plt.subplots_adjust(top=0.85)
+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., first_signal_bin=0):
+vmin=0, vmax=1.3 * 10 ** 7, mask_noise=True, noise_threshold=1., first_signal_bin=0, box=False):
 """
 Draw channel data as a PPI plot.
 Parameters
 ----------
 noise_threshold : int
 Threshold to use in the noise masking routine.
 first_signal_bin : int
 First signal bin. Can be used to fix analog bin shift of Licel channels.
 """
-x, y = self._polar_to_ground(self.z, self.azimuth_angles, self.zenith_angles)
+x, y = self._polar_to_ground(self.z/1000., self.azimuth_angles, self.zenith_angles)
 self.draw_projection(ax1, x, y, cmap=cmap, signal_type=signal_type,
 z_min=z_min, z_max=z_max, add_colorbar=add_colorbar, cmap_label=cmap_label,
 cb_format=cb_format, vmin=vmin, vmax=vmax, mask_noise=mask_noise,
 noise_threshold=noise_threshold, first_signal_bin=first_signal_bin)
+if box:
+ax1.set_xlim(-z_max / 1000., z_max / 1000.)
+ax1.set_ylim(-z_max / 1000., z_max / 1000.)
+ax1.set_ylabel('South-North (km)')
+ax1.set_xlabel('West-East (km)')
 def draw_rhi(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., first_signal_bin=0):
+vmin=0, vmax=1.3 * 10 ** 7, mask_noise=True, noise_threshold=1., first_signal_bin=0, box=False):
 """
 Draw channel data as a PPI plot.
 Parameters
 ----------
 noise_threshold : int
 Threshold to use in the noise masking routine.
 first_signal_bin : int
 First signal bin. Can be used to fix analog bin shift of Licel channels.
 """
-x, y = self._polar_to_cross_section(self.z, self.azimuth_angles, self.zenith_angles, np.mean(self.azimuth_angles))
+projection_angle = np.mean(self.azimuth_angles)
+x, y = self._polar_to_cross_section(self.z / 1000., self.azimuth_angles, self.zenith_angles, projection_angle)
 self.draw_projection(ax1, x, y, cmap=cmap, signal_type=signal_type,
 z_min=z_min, z_max=z_max, add_colorbar=add_colorbar, cmap_label=cmap_label,
 cb_format=cb_format, vmin=vmin, vmax=vmax, mask_noise=mask_noise,
 noise_threshold=noise_threshold, first_signal_bin=first_signal_bin)
+padding = 0.5 # km
+if box:
+ax1.set_xlim(-z_max / 1000. - padding, z_max / 1000. + padding)
+ax1.set_ylim(-padding, z_max / 1000. + padding)
+ax1.set_xlabel('Distance (km)')
+ax1.set_ylabel('Height a.l. (km)')
+return projection_angle
 def draw_projection(self, ax1, x, y, 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.,
 first_signal_bin=0):
 data_max_idx = z_max_idx + first_signal_bin
 im1 = ax1.pcolormesh(x[:, z_min_idx:z_max_idx], y[:, z_min_idx:z_max_idx], data[:, data_min_idx:data_max_idx],
 cmap=cmap, vmin=vmin, vmax=vmax)
-ax1.axis('equal')
+ax1.set(adjustable='box', aspect='equal')
 if add_colorbar:
 if cb_format:
 cb1 = plt.colorbar(im1, format=cb_format)
 else:
 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)
+# Centered axis in center: https://stackoverflow.com/a/31558968
+# Move left y-axis and bottim x-axis to centre, passing through (0,0)
+# ax1.spines['left'].set_position('center')
+# ax1.spines['bottom'].set_position('center')
+#
+# # Eliminate upper and right axes
+# ax1.spines['right'].set_color('none')
+# ax1.spines['top'].set_color('none')
+#
+# # Show ticks in the left and lower axes only
+# ax1.xaxis.set_ticks_position('bottom')
+# ax1.yaxis.set_ticks_position('left')
 @staticmethod
 def _polar_to_ground(z, azimuth, zenith):
 """
 Convert polar coordinates to cartesian project for a PPI scan
 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]
+x = Z_ground * np.sin(azimuth_rad)[:, np.newaxis]
-y = Z_ground * np.sin(azimuth_rad)[:, np.newaxis]
+y = Z_ground * np.cos(azimuth_rad)[:, np.newaxis]
 return x, y
 @staticmethod
 def _polar_to_cross_section(z, azimuth, zenith, cross_section_azimuth):
 """
 zenith_rad = np.deg2rad(zenith)
 # The angle between measurements and the cross section plance
-azimuth_difference_rad = np.deg2rad(azimuth) - cross_section_azimuth
+azimuth_difference_rad = np.deg2rad(azimuth) - np.deg2rad(cross_section_azimuth)
 # Generate the mesh
 Z, Azimuth_differences = np.meshgrid(z, azimuth_difference_rad)
-Z_cross_section = Z * np.sin(Azimuth_differences)
+x = Z * np.sin(zenith_rad)[:, np.newaxis] * np.cos(Azimuth_differences)
-x = Z_cross_section * np.sin(zenith_rad)[:, np.newaxis]
+y = Z * np.cos(zenith_rad)[:, np.newaxis]
-y = Z_cross_section * np.cos(zenith_rad)[:, np.newaxis]
 return x, y
 class ScanningLidarMeasurement(LicelLidarMeasurement):
 """ A class representing a scanning measurement set.

comparison: atmospheric_lidar/raymetrics.py

atmospheric_lidar/raymetrics.py

Mercurial > public > atmospheric_lidar / file comparison

comparison: atmospheric_lidar/raymetrics.py

atmospheric_lidar/raymetrics.py