Atmospheric Lidar Package: atmospheric_lidar/raymetrics.py@53124a886152 (annotated)

atmospheric_lidar/raymetrics.py@53124a886152 (annotated)

atmospheric_lidar/raymetrics.py

Wed, 18 Nov 2020 15:11:18 +0000

author: George Doxastakis <gdoxastakis.ee@gmail.com>
date: Wed, 18 Nov 2020 15:11:18 +0000
changeset 207: 53124a886152
parent 194: 809190df0dc8
child 209: f3bec37bfcd0
permissions: -rw-r--r--

Added custom_info field support

 """ Code to read Raymetrics version of Licel binary files."""
 import logging
 import matplotlib as mpl
 import numpy as np
 from matplotlib import pyplot as plt
 from .licel import LicelFile, LicelLidarMeasurement, LicelChannel, PhotodiodeChannel
 logger = logging.getLogger(__name__)
 class ScanningFile(LicelFile):
     """ Raymetrics is using a custom version of licel file format to store scanning lidar measurements.
     The file includes one extra line describing the scan strategy of the dataset. The extra parameters are:
     `azimuth_start`
         Start azimuth angle for the scan, relative to instrument zero position (degrees).
     `azimuth_stop`
         Stop azimuth angle for the scan, relative to instrument zero position (degrees).
     `azimuth_step`
         Step of the azimuth scan (degrees).
     `zenith_start`
         Start zenith angle for the scan, relative to *nadir* (degrees). Take care that this is actually
         nadir angle. Vertical measurements correspond to -90.
     `zenith_stop`
         Stop zenith angle for the scan, relative to *nadir* (degrees). Take care that this is actually
         nadir angle. Vertical measurements correspond to -90.
     `zenith_step`
         Step of the zenith scan (degrees).
     `azimuth_offset`
         Offset of instrument zero from North (degrees). Using this value you can convert `azimuth_start` and
         `azimuth_stop` to absolute values.
     Moreover, four new parameters are added in the second line of the file:
     `zenith_angle`
         Zenith angle of the current file. Take care that this is actually
         nadir angle. Vertical measurements correspond to -90.
     `azimuth_angle`
         Azimuth angle of the current file. Value relative to instrument zero position.
     `temperature`
         Ambient temperature (degrees C)
     `pressure`
         Ambient pressure (hPa)
     """
     # Specifications of the header lines.
     licel_file_header_format = ['filename',
                                 'start_date start_time end_date end_time altitude longitude latitude zenith_angle azimuth_angle temperature pressure custom_info',
                                 # Appart from Site that is read manually
                                 'azimuth_start azimuth_stop azimuth_step zenith_start zenith_stop zenith_step azimuth_offset',
                                 'LS1 rate_1 LS2 rate_2 number_of_datasets', ]
     # 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
     skip_scan_overview_line = False  # Skip the 3d line containing azimuth_start, stop etc. Used to overcome a bug in some output files.
     def _read_rest_of_header(self, f):
         """ Read the third and fourth row of  of the header lines.
         The first two rows are read in the licel class.
         Parameters
         ----------
         f : file
            An open file-like object.
         Returns
         -------
         raw_info : dict
            A dictionary containing all parameters of the third and fourth line of the header.
         """
         raw_info = {}
         third_line = f.readline().decode()
         self.header_lines.append(third_line.strip())
         raw_info.update(self.match_lines(third_line, self.licel_file_header_format[2]))
         fourth_line = f.readline().decode()
         self.header_lines.append(fourth_line.strip())
         raw_info.update(self.match_lines(fourth_line, self.licel_file_header_format[3]))
         return raw_info
     def _assign_properties(self):
         """ Assign scanning-specific parameters found in the header as object properties."""
         super(ScanningFile, self)._assign_properties()
         self.azimuth_angle_raw = float(self.raw_info['azimuth_angle'])
         self.temperature = float(self.raw_info['temperature'])
         self.pressure = float(self.raw_info['pressure'])
         if not self.skip_scan_overview_line:
             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_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'])
         else:
             self.azimuth_start_raw = np.nan
             self.azimuth_stop_raw = np.nan
             self.azimuth_step = np.nan
             self.zenith_start_raw = np.nan
             self.zenith_stop_raw = np.nan
             self.zenith_step = np.nan
             self.azimuth_offset = 0
         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.')
             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
         try:
             self.custom_info = self.raw_info['custom_info'].strip('"')
         except KeyError:
             self.custom_info = None
     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. """
     def __init__(self):
         super(ScanningChannel, self).__init__()
         self.azimuth_start = None
         self.azimuth_stop = None
         self.azimuth_step = None
         self.zenith_start = None
         self.zenith_stop = None
         self.zenith_step = None
         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,
         i.e. that all files belong to the same measurements set.
         Parameters
         ----------
         current_file : ScanningFile object
            A ScanningFile object being imported
         file_channel : LicelChannelData object
            A specific LicelChannelData object holding data found in the file.
         """
         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,
                  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 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)
         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 ({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, grid=(1, 4), ax1_position=(0, 0), ax1_span=2, ax2_position=(0, 2), ax2_span=2):
         """
         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 = plt.subplot2grid(grid, ax1_position, colspan=ax1_span)
         ax2 = plt.subplot2grid(grid, ax2_position, colspan=ax2_span)
         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, box=False):
         """
         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 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.
         """
         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, box=False):
         """
         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.
         """
         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.,
                         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:
                 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)
         data_min_idx = z_min_idx + first_signal_bin
         data_max_idx = z_max_idx + first_signal_bin
         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.set(adjustable='box', aspect='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)
         # 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):
         """
         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.sin(azimuth_rad)[:, np.newaxis]
         y = Z_ground * np.cos(azimuth_rad)[:, np.newaxis]
         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) - np.deg2rad(cross_section_azimuth)
         # Generate the mesh
         Z, Azimuth_differences = np.meshgrid(z, azimuth_difference_rad)
         x = Z * np.sin(zenith_rad)[:, np.newaxis] * np.cos(Azimuth_differences)
         y = Z * np.cos(zenith_rad)[:, np.newaxis]
         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 ScanningFileMissingLine(ScanningFile):
     skip_scan_overview_line = True
 class ScanningLidarMeasurement(LicelLidarMeasurement):
     """ A class representing a scanning measurement set.
     It useses `ScanningFile` and `ScanningChannel` classes for handling the data.
     """
     file_class = ScanningFile
     channel_class = ScanningChannel
     photodiode_class = PhotodiodeChannel
 class FixedPointingFile(LicelFile):
     """ Raymetrics is using a custom version of licel file format to store
     vertical lidar measurements.
     `temperature`
         Ambient temperature (degrees C)
     `pressure`
         Ambient pressure (hPa)
     """
     # Specifications of the header lines.
     licel_file_header_format = ['filename',
                                 'start_date start_time end_date end_time altitude longitude latitude zenith_angle azimuth_angle temperature pressure custom_info',
                                 # Appart from Site that is read manually
                                 'LS1 rate_1 LS2 rate_2 number_of_datasets', ]
     fix_zenith_angle = True
     def _assign_properties(self):
         """ Assign scanning-specific parameters found in the header as object properties."""
         super(FixedPointingFile, self)._assign_properties()
         self.temperature = float(self.raw_info['temperature'])
         self.pressure = float(self.raw_info['pressure'])
         self.azimuth_angle = float(self.raw_info['azimuth_angle'])
         try:
             self.custom_info = self.raw_info['custom_info'].strip('"')
         except KeyError:
             self.custom_info = None
 class FixedPointingChannel(LicelChannel):
     """ A class representing measurements of a specific lidar channel, during a fixed pointing measurement. """
     def __init__(self):
         super(FixedPointingChannel, self).__init__()
         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(FixedPointingChannel, 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)
 class FixedPointingLidarMeasurement(LicelLidarMeasurement):
     file_class = FixedPointingFile
     channel_class = FixedPointingChannel

Mercurial > public > atmospheric_lidar / annotate

atmospheric_lidar/raymetrics.py@53124a886152 (annotated)

atmospheric_lidar/raymetrics.py