Mercurial > public > atmospheric_lidar / changeset
changeset
Merge.
Fri, 14 Dec 2018 16:16:19 +0200
- author
- Iannis <i.binietoglou@impworks.gr>
- date
- Fri, 14 Dec 2018 16:16:19 +0200
- changeset 175
- 6960ccccc713
- parent 171
- 322117fa8391 (current diff)
- parent 174
- 16cf4d961d02 (diff)
- child 176
- 34866a2a1aa5
Merge.
--- a/atmospheric_lidar/generic.py Thu Nov 29 22:26:22 2018 +0200 +++ b/atmospheric_lidar/generic.py Fri Dec 14 16:16:19 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. @@ -1116,6 +1140,7 @@ 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
--- a/atmospheric_lidar/licel.py Thu Nov 29 22:26:22 2018 +0200 +++ b/atmospheric_lidar/licel.py Fri Dec 14 16:16:19 2018 +0200 @@ -154,7 +154,7 @@ channel_data_class = LicelChannelData - def __init__(self, file_path, use_id_as_name=False, licel_timezone="UTC", import_now=True): + def __init__(self, file_path, use_id_as_name=False, licel_timezone="UTC", import_now=True, fix_zenith_angle=False): """ This is run when creating a new object. @@ -168,10 +168,12 @@ licel_timezone : str The timezone of dates found in the Licel files. Should match the available timezones in the TZ database. - import_file : bool + import_now : bool 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) @@ -180,6 +182,7 @@ 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() @@ -283,7 +286,30 @@ self.altitude = float(self.raw_info['altitude']) self.longitude = float(self.raw_info['longitude']) self.latitude = float(self.raw_info['latitude']) - self.zenith_angle = float(self.raw_info['zenith_angle']) + + self.zenith_angle_raw = float(self.raw_info['zenith_angle']) + + if self.fix_zenith_angle: + self.zenith_angle = self._correct_zenith_angle(self.zenith_angle_raw) + else: + self.zenith_angle = self.zenith_angle_raw + + @staticmethod + def _correct_zenith_angle(zenith_angle): + """ Correct zenith angle from Raymetrics convention (zenith = -90 degrees). + + Parameters + ---------- + zenith_angle : float + Zenith angle in Raymetrics convention. + + Returns + ------- + corrected_angle : float + Corrected zenith angle. + """ + corrected_angle = 90 - zenith_angle + return corrected_angle def _read_second_header_line(self, f): """ Read the second line of a licel file. """ @@ -439,7 +465,7 @@ channel_class = LicelChannel photodiode_class = PhotodiodeChannel - def __init__(self, file_list=None, use_id_as_name=False, licel_timezone='UTC'): + def __init__(self, file_list=None, use_id_as_name=False, licel_timezone='UTC', fix_zenith_angle=False): self.raw_info = {} # Keep the raw info from the files self.durations = {} # Keep the duration of the files self.laser_shots = [] @@ -447,6 +473,7 @@ 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) @@ -456,7 +483,8 @@ logger.warning("File has been imported already: %s" % filename) 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) + 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) 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 Thu Nov 29 22:26:22 2018 +0200 +++ b/atmospheric_lidar/raymetrics.py Fri Dec 14 16:16:19 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__) @@ -88,17 +92,29 @@ 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['altitude']) + self.azimuth_angle = float(self.raw_info['azimuth_angle']) self.temperature = float(self.raw_info['temperature']) self.pressure = float(self.raw_info['pressure']) - self.azimuth_start = float(self.raw_info['azimuth_start']) - self.azimuth_stop = float(self.raw_info['azimuth_stop']) + 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 = float(self.raw_info['zenith_start']) - self.zenith_stop = float(self.raw_info['zenith_stop']) + 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.") + # TODO: Apply azimuth offset correction. + self.azimuth_start = self.azimuth_start_raw + self.azimuth_stop = self.azimuth_stop_raw + + if self.fix_zenith_angle: + self.zenith_start = self._correct_zenith_angle(self.zenith_start_raw) + self.zenith_stop = self._correct_zenith_angle(self.zenith_stop_raw) + else: + self.zenith_start = self.zenith_start_raw + self.zenith_stop = self.zenith_stop_raw + class ScanningChannel(LicelChannel): """ A class representing measurements of a specific lidar channel, during a scanning measurement. """ @@ -141,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.
