Mercurial > public > atmospheric_lidar / changeset
changeset
First attempt for RHI plot from scanning lidar. Axes are still reversed.
Fri, 14 Dec 2018 19:23:38 +0200
- author
- Iannis <i.binietoglou@impworks.gr>
- date
- Fri, 14 Dec 2018 19:23:38 +0200
- changeset 176
- 34866a2a1aa5
- parent 175
- 6960ccccc713
- child 177
- 580821c7487c
First attempt for RHI plot from scanning lidar. Axes are still reversed.
| 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 16:16:19 2018 +0200 +++ b/atmospheric_lidar/licel.py Fri Dec 14 19:23:38 2018 +0200 @@ -154,7 +154,10 @@ channel_data_class = LicelChannelData - def __init__(self, file_path, use_id_as_name=False, licel_timezone="UTC", import_now=True, fix_zenith_angle=False): + # If True, it corrects the old Raymetrics convention of zenith angle definition (zenith = -90 degrees) + fix_zenith_angle = False + + def __init__(self, file_path, use_id_as_name=False, licel_timezone="UTC", import_now=True): """ This is run when creating a new object. @@ -172,8 +175,6 @@ If True, the header and data are read immediately. If not, the user has to call the corresponding methods directly. This is used to speed up reading files when only header information are required. - fix_zenith_angle : bool - If True, it corrects the old Raymetrics convention of zenith angle definition (zenith = -90 degrees) """ self.file_path = file_path self.file_name = os.path.basename(file_path) @@ -182,7 +183,6 @@ self.start_time = None self.stop_time = None self.licel_timezone = licel_timezone - self.fix_zenith_angle = fix_zenith_angle if import_now: self.import_file() @@ -289,7 +289,10 @@ self.zenith_angle_raw = float(self.raw_info['zenith_angle']) + logger.debug('Fix zenith angle? %s' % self.fix_zenith_angle) + if self.fix_zenith_angle: + logger.debug('Fixing zenith angle.') self.zenith_angle = self._correct_zenith_angle(self.zenith_angle_raw) else: self.zenith_angle = self.zenith_angle_raw @@ -308,7 +311,7 @@ corrected_angle : float Corrected zenith angle. """ - corrected_angle = 90 - zenith_angle + corrected_angle = 90 + zenith_angle return corrected_angle def _read_second_header_line(self, f): @@ -465,7 +468,7 @@ channel_class = LicelChannel photodiode_class = PhotodiodeChannel - def __init__(self, file_list=None, use_id_as_name=False, licel_timezone='UTC', fix_zenith_angle=False): + def __init__(self, file_list=None, use_id_as_name=False, licel_timezone='UTC'): self.raw_info = {} # Keep the raw info from the files self.durations = {} # Keep the duration of the files self.laser_shots = [] @@ -473,7 +476,6 @@ self.use_id_as_name = use_id_as_name self.licel_timezone = licel_timezone self.photodiodes = {} - self.fix_zenith_angle = fix_zenith_angle super(LicelLidarMeasurement, self).__init__(file_list) @@ -484,7 +486,7 @@ else: logger.debug('Importing file {0}'.format(filename)) current_file = self.file_class(filename, use_id_as_name=self.use_id_as_name, - licel_timezone=self.licel_timezone, fix_zenith_angle=self.fix_zenith_angle) + licel_timezone=self.licel_timezone) self.raw_info[current_file.file_path] = current_file.raw_info self.durations[current_file.file_path] = current_file.duration()
--- a/atmospheric_lidar/raymetrics.py Fri Dec 14 16:16:19 2018 +0200 +++ b/atmospheric_lidar/raymetrics.py Fri Dec 14 19:23:38 2018 +0200 @@ -65,6 +65,8 @@ # Specifications of the channel lines in the header licel_file_channel_format = 'active analog_photon laser_used number_of_datapoints 1 HV bin_width wavelength d1 d2 d3 d4 ADCbits number_of_shots discriminator ID' + fix_zenith_angle = True + def _read_rest_of_header(self, f): """ Read the third and fourth row of of the header lines. @@ -109,6 +111,7 @@ 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) else: @@ -203,9 +206,55 @@ if show_plot: plt.show() + def plot_rhi(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_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") + + 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.): + 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): """ Draw channel data as a PPI plot. @@ -213,6 +262,10 @@ ---------- ax1 : axis object The axis object to draw. + x : array + X axis coordinates + y : array + Y axis coordinates cmap : cmap An instance of a matplotlib colormap to use. signal_type : str @@ -235,6 +288,101 @@ If True, remove noisy bins. 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) + + 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) + + 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): + """ + Draw channel data as a PPI plot. + + Parameters + ---------- + ax1 : axis object + The axis object to draw. + x : array + X axis coordinates + y : array + Y axis coordinates + cmap : cmap + An instance of a matplotlib colormap to use. + signal_type : str + If 'rc', the range corrected signal is plotted. 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. + 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)) + + 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) + + 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): + """ + Draw channel data as a PPI plot. + + Parameters + ---------- + ax1 : axis object + The axis object to draw. + x : array + X axis coordinates + y : array + Y axis coordiantes + 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. + first_signal_bin : int + First signal bin. Can be used to fix analog bin shift of Licel channels. """ if signal_type == 'rc': if len(self.rc) == 0: @@ -250,9 +398,11 @@ 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) + data_min_idx = z_min_idx + first_signal_bin + data_max_idx = z_max_idx + first_signal_bin - im1 = ax1.pcolormesh(x, y, data[:, z_min_idx:z_max_idx], cmap=cmap, vmin=vmin, vmax=vmax) + 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') @@ -302,6 +452,45 @@ return x, y + @staticmethod + def _polar_to_cross_section(z, azimuth, zenith, cross_section_azimuth): + """ + 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 + cross_section_azimuth : float + Azimuth angle of plane in degrees + + Returns + ------- + x : array + X axis in meters + y : array + Y axis in meters + """ + + zenith_rad = np.deg2rad(zenith) + + # The angle between measurements and the cross section plance + azimuth_difference_rad = np.deg2rad(azimuth) - 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] + + return x, y + + class ScanningLidarMeasurement(LicelLidarMeasurement): """ A class representing a scanning measurement set.
