Mercurial > public > atmospheric_lidar / changeset
changeset
Azimuth and elevation corrections, and some plots.
Thu, 10 Jan 2019 17:16:39 +0200
- author
- Ioannis <ioannis@inoe.ro>
- date
- Thu, 10 Jan 2019 17:16:39 +0200
- changeset 177
- 580821c7487c
- parent 176
- 34866a2a1aa5
- child 178
- 804313c093a8
Azimuth and elevation corrections, and some plots.
| atmospheric_lidar/licel.py | file | annotate | diff | comparison | revisions | |
| atmospheric_lidar/raymetrics.py | file | annotate | diff | comparison | revisions |
--- a/atmospheric_lidar/licel.py Fri Dec 14 19:23:38 2018 +0200 +++ b/atmospheric_lidar/licel.py Thu Jan 10 17:16:39 2019 +0200 @@ -186,6 +186,8 @@ if import_now: self.import_file() + else: + self.import_header_only() def import_file(self): """ Read the header info and data of the Licel file.
--- a/atmospheric_lidar/raymetrics.py Fri Dec 14 19:23:38 2018 +0200 +++ b/atmospheric_lidar/raymetrics.py Thu Jan 10 17:16:39 2019 +0200 @@ -94,7 +94,7 @@ 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']) @@ -105,10 +105,9 @@ 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.") - # TODO: Apply azimuth offset correction. - self.azimuth_start = self.azimuth_start_raw - self.azimuth_stop = self.azimuth_stop_raw + self.azimuth_angle = (self.azimuth_angle_raw + self.azimuth_offset) % 360 + self.azimuth_start = (self.azimuth_start_raw + self.azimuth_offset) % 360 + self.azimuth_stop = (self.azimuth_stop_raw + self.azimuth_offset) % 360 if self.fix_zenith_angle: logger.debug('Fixing zenith start and zenith stop angles.') @@ -134,12 +133,16 @@ 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, @@ -155,9 +158,13 @@ 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, @@ -241,20 +248,78 @@ 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. @@ -291,16 +356,23 @@ 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. @@ -337,13 +409,25 @@ 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., @@ -404,7 +488,7 @@ 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: @@ -419,6 +503,20 @@ 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): """ @@ -447,8 +545,8 @@ 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] + x = Z_ground * np.sin(azimuth_rad)[:, np.newaxis] + y = Z_ground * np.cos(azimuth_rad)[:, np.newaxis] return x, y @@ -479,19 +577,17 @@ 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_cross_section * np.sin(zenith_rad)[:, np.newaxis] - y = Z_cross_section * np.cos(zenith_rad)[:, np.newaxis] + x = Z * np.sin(zenith_rad)[:, np.newaxis] * np.cos(Azimuth_differences) + y = Z * np.cos(zenith_rad)[:, np.newaxis] return x, y - class ScanningLidarMeasurement(LicelLidarMeasurement): """ A class representing a scanning measurement set.
