Mercurial > public > atmospheric_lidar / changeset

--- 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):
 

--- 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'):
 

--- 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