# HG changeset patch # User Iannis # Date 1571755608 -10800 # Node ID 12c895b99a3ca03565ad1a36653d7bc16be041d0 # Parent a6812046c8b3dfd8afd4d924c9360ba4a603d305 Initial thoughts on plotting lidar data on lat/lon coordinates. diff -r a6812046c8b3 -r 12c895b99a3c atmospheric_lidar/licel.py --- a/atmospheric_lidar/licel.py Mon Oct 21 15:08:18 2019 +0300 +++ b/atmospheric_lidar/licel.py Tue Oct 22 17:46:48 2019 +0300 @@ -281,9 +281,6 @@ self.start_time = local_start_time.astimezone(pytz.utc) self.stop_time = local_stop_time.astimezone(pytz.utc) - self.latitude = float(raw_info['latitude']) - self.longitude = float(raw_info['longitude']) - # Read the rest of the header. for c1 in range(int(raw_info['number_of_datasets'])): channel_info.append(self.match_lines(f.readline().decode(), self.licel_file_channel_format)) @@ -443,6 +440,8 @@ self._assign_unique_property('active', file_channel.active) self._assign_unique_property('laser_used', file_channel.laser_used) self._assign_unique_property('analog_photon_string', file_channel.analog_photon_string) + self._assign_unique_property('latitude', current_file.latitude) + self._assign_unique_property('longitude', current_file.longitude) def _assign_unique_property(self, property_name, value): diff -r a6812046c8b3 -r 12c895b99a3c atmospheric_lidar/licelv2.py --- a/atmospheric_lidar/licelv2.py Mon Oct 21 15:08:18 2019 +0300 +++ b/atmospheric_lidar/licelv2.py Tue Oct 22 17:46:48 2019 +0300 @@ -263,7 +263,7 @@ return len(self.standard_deviations) != 0 -class LicelChannelV2(LicelChannel): +class LicelV2Channel(LicelChannel): def __init__(self): @@ -272,23 +272,17 @@ self.bin_shift = None self.bin_shift_decimal_places = None - super(LicelChannelV2, self).__init__() + self.zenith_angles = [] + self.azimuth_angles = [] + + super(LicelV2Channel, self).__init__() def append_file(self, current_file, file_channel): """ Append file to the current object """ - self._assign_properties(current_file, file_channel) - - self.binwidth = self.resolution * 2. / c # in seconds - self.z = file_channel.z - - self.data[current_file.start_time] = file_channel.data - self.raw_info.append(file_channel.raw_info) - self.laser_shots.append(file_channel.number_of_shots) - self.duration.append(file_channel.duration) - self.number_of_shots.append(file_channel.number_of_shots) - self.discriminator.append(file_channel.discriminator) - self.hv.append(file_channel.hv) + super(LicelV2Channel, self).append_file(current_file, file_channel) + self.zenith_angles.append(current_file.zenith_angle) + self.azimuth_angles.append(current_file.azimuth_angle) def _assign_properties(self, current_file, file_channel): @@ -297,7 +291,7 @@ self._assign_unique_property('bin_shift', file_channel.bin_shift) self._assign_unique_property('bin_shift_decimal_places', file_channel.bin_shift_decimal_places) - super(LicelChannelV2, self)._assign_properties(current_file, file_channel) + super(LicelV2Channel, self)._assign_properties(current_file, file_channel) @property def is_analog(self): @@ -317,10 +311,10 @@ class LicelLidarMeasurementV2(LicelLidarMeasurement): file_class = LicelFileV2 - channel_class = LicelChannelV2 + channel_class = LicelV2Channel photodiode_class = PhotodiodeChannel powermeter_class = PhotodiodeChannel - standard_deviation_class = LicelChannelV2 + standard_deviation_class = LicelV2Channel def __init__(self, file_list=None, licel_timezone='UTC'): diff -r a6812046c8b3 -r 12c895b99a3c atmospheric_lidar/raymetrics.py --- a/atmospheric_lidar/raymetrics.py Mon Oct 21 15:08:18 2019 +0300 +++ b/atmospheric_lidar/raymetrics.py Tue Oct 22 17:46:48 2019 +0300 @@ -1,8 +1,8 @@ """ Code to read Raymetrics version of Licel binary files.""" import logging +import matplotlib as mpl import numpy as np -import matplotlib as mpl from matplotlib import pyplot as plt from .licel import LicelFile, LicelLidarMeasurement, LicelChannel, PhotodiodeChannel @@ -117,6 +117,59 @@ self.zenith_start = self.zenith_start_raw self.zenith_stop = self.zenith_stop_raw + def get_coordinates(self, channel_name): + """ + Calculate the lat, lon, z coordinates for each measurement point. + + Parameters + ---------- + channel_name : str + The name of the channel. Only the channel object knows about the + range resolution. + + Returns + ------- + lat : array + Latitude array + lon : array + Longitude array + z : array + Altitude array in meters + """ + R_earth = 6378137 # Earth radius in meters + + # Shortcuts to make equations cleaner + lat_center = self.latitude + lon_center = self.longitude + r = self.channels[channel_name].z + azimuth = self.azimuth_angle + zenith = self.zenith_angle + + # Convert all angles to radiants + zenith_rad = np.deg2rad(zenith)[:, np.newaxis] + azimuth_rad = np.deg2rad(azimuth)[:, np.newaxis] + + lat_center_rad = np.deg2rad(lat_center) + lon_center_rad = np.deg2rad(lon_center) + + # Generate the mesh + R, Zeniths = np.meshgrid(r, zenith_rad) + R_ground = R * np.sin(Zeniths) + Z = R * np.cos(Zeniths) + + # Equations from https://www.movable-type.co.uk/scripts/latlong.html + delta = R_ground / R_earth + lat_out_rad = np.arcsin(np.sin(lat_center_rad) * np.cos(delta) + + np.cos(lat_center_rad) * np.sin(delta) * np.cos(azimuth_rad)) + lon_out_rad = lon_center_rad + np.arctan2(np.sin(azimuth_rad) * np.sin(delta) * np.cos(lat_center_rad), + np.cos(delta) - np.sin(lat_center_rad) * np.sin(lat_out_rad)) + + # Convert back to degrees + lat_out = np.rad2deg(lat_out_rad) + lon_out = np.rad2deg(lon_out_rad) + + return lat_out, lon_out, Z + class ScanningChannel(LicelChannel): """ A class representing measurements of a specific lidar channel, during a scanning measurement. """ @@ -168,7 +221,7 @@ 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.): + 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. @@ -214,7 +267,7 @@ 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.): + 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. @@ -248,8 +301,9 @@ fig = plt.figure(figsize=figsize) ax1 = fig.add_subplot(111) - 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) + 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) @@ -260,7 +314,8 @@ 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): + 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. @@ -299,9 +354,10 @@ 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) - + 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'), @@ -316,7 +372,6 @@ 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, box=False): @@ -356,7 +411,7 @@ 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/1000., 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, @@ -417,7 +472,7 @@ 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 + padding = 0.5 # km if box: ax1.set_xlim(-z_max / 1000. - padding, z_max / 1000. + padding) @@ -516,7 +571,6 @@ # ax1.xaxis.set_ticks_position('bottom') # ax1.yaxis.set_ticks_position('left') - @staticmethod def _polar_to_ground(z, azimuth, zenith): """ @@ -587,6 +641,53 @@ return x, y + def get_coordinates(self,): + """ + Calculate the lat, lon, z coordinates for each measurement point. + + Returns + ------- + lat : array + Latitude array + lon : array + Longitude array + z : array + Altitude array in meters + """ + R_earth = 6378137 # Earth radius in meters + + # Shortcuts to make equations cleaner + lat_center = self.latitude + lon_center = self.longitude + r = self.z + azimuth = self.azimuth_angles + zenith = self.zenith_angles + + # Convert all angles to radiants + zenith_rad = np.deg2rad(zenith)[:, np.newaxis] + azimuth_rad = np.deg2rad(azimuth)[:, np.newaxis] + + lat_center_rad = np.deg2rad(lat_center) + lon_center_rad = np.deg2rad(lon_center) + + # Generate the mesh + R, Zeniths = np.meshgrid(r, zenith_rad) + R_ground = R * np.sin(Zeniths) + Z = R * np.cos(Zeniths) + + # Equations from https://www.movable-type.co.uk/scripts/latlong.html + delta = R_ground / R_earth + lat_out_rad = np.arcsin(np.sin(lat_center_rad) * np.cos(delta) + + np.cos(lat_center_rad) * np.sin(delta) * np.cos(azimuth_rad)) + lon_out_rad = lon_center_rad + np.arctan2(np.sin(azimuth_rad) * np.sin(delta) * np.cos(lat_center_rad), + np.cos(delta) - np.sin(lat_center_rad) * np.sin(lat_out_rad)) + + # Convert back to degrees + lat_out = np.rad2deg(lat_out_rad) + lon_out = np.rad2deg(lon_out_rad) + + return lat_out, lon_out, Z + class ScanningLidarMeasurement(LicelLidarMeasurement): """ A class representing a scanning measurement set. @@ -645,4 +746,4 @@ class FixedPointingLidarMeasurement(LicelLidarMeasurement): file_class = FixedPointingFile - channel_class = FixedPointingChannel \ No newline at end of file + channel_class = FixedPointingChannel