Mercurial > public > atmospheric_lidar / changeset
changeset
Added vertical lidar measurement class.
Tue, 25 Sep 2018 13:30:45 +0300
- author
- Ioannis <ioannis@inoe.ro>
- date
- Tue, 25 Sep 2018 13:30:45 +0300
- changeset 158
- 5d99df1c7e81
- parent 157
- 03b470b0a05f (diff)
- parent 152
- 3e4e6472b88b (current diff)
- child 159
- c31f5388b482
Added vertical lidar measurement class.
--- a/atmospheric_lidar/diva.py Thu Sep 06 16:51:13 2018 +0300 +++ b/atmospheric_lidar/diva.py Tue Sep 25 13:30:45 2018 +0300 @@ -8,13 +8,16 @@ import datetime import os import numpy as np +import logging import pytz from .generic import BaseLidarMeasurement +logger = logging.getLogger(__name__) -class DivaOutput(BaseLidarMeasurement): + +class DivaMixin: def save_as_diva_netcdf(self, output_path, parameter_file): """ Save the current data in the 'draft' DIVA format. """ @@ -413,3 +416,48 @@ raise ValueError('Emission polarization not one of {0}: {1}'.format(choices.keys(), pol_string)) return choices[pol_string] + + +class DivaLidarMeasurement(BaseLidarMeasurement): + + def __init__(self, file_list): + """ + This is run when creating a new object. + + Parameters + ---------- + file_list : list or str + A list of the full paths to the input file(s). + """ + if isinstance(file_list, str): + file_list = [file_list, ] + super(DivaLidarMeasurement, self).__init__(file_list=file_list) + + def _import_file(self, filename): + """ Import data from a single DIVA file. """ + + 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) + self.durations[current_file.filename] = current_file.duration() + + + self._create_or_append_channel(current_file) + file_laser_shots = [] + self.laser_shots.append(file_laser_shots) + + +class DivaChannel(object): + + def __init__(self, file_name, channel_group): + """ This is run when first creating the object. + + Parameters + ---------- + file_name : str + The filename of the diva dataset. + channel_group : str + The name of the netCDF4 group that holds the channel data. + """ + self.file_name = file_name + self.channel_group = channel_group +
--- a/atmospheric_lidar/generic.py Thu Sep 06 16:51:13 2018 +0300 +++ b/atmospheric_lidar/generic.py Tue Sep 25 13:30:45 2018 +0300 @@ -34,8 +34,8 @@ Parameters ---------- - file_list : list - A list of the full paths to the input file. + file_list : list or str + A list of the full paths to the input file(s). """ self.info = {} self.dimensions = {}
--- a/atmospheric_lidar/licel.py Thu Sep 06 16:51:13 2018 +0300 +++ b/atmospheric_lidar/licel.py Tue Sep 25 13:30:45 2018 +0300 @@ -1,11 +1,13 @@ import datetime import logging import copy +import os import numpy as np import pytz -from generic import BaseLidarMeasurement, LidarChannel +from .generic import BaseLidarMeasurement, LidarChannel +from .diva import DivaMixin logger = logging.getLogger(__name__) @@ -171,7 +173,9 @@ corresponding methods directly. This is used to speed up reading files when only header information are required. """ - self.filename = file_path + self.file_path = file_path + self.file_name = os.path.basename(file_path) + self.use_id_as_name = use_id_as_name self.start_time = None self.stop_time = None @@ -186,7 +190,7 @@ channels = {} photodiodes = {} - with open(self.filename, 'rb') as f: + with open(self.file_path, 'rb') as f: self.read_header(f) @@ -200,7 +204,7 @@ b = np.fromfile(f, 'b', 1) if (a[0] != 13) | (b[0] != 10): - logger.warning("No end of line found after record. File could be corrupt: %s" % self.filename) + logger.warning("No end of line found after record. File could be corrupt: %s" % self.file_path) logger.warning('a: {0}, b: {1}.'.format(a, b)) channel = self.channel_data_class(current_channel_info, raw_data, self.duration(), @@ -295,6 +299,7 @@ site_name = second_line.split('/')[0][:-2] clean_site_name = site_name.strip() raw_info['site'] = clean_site_name + self.site = clean_site_name raw_info.update(self.match_lines(second_line[len(clean_site_name) + 1:], self.licel_file_header_format[1])) return raw_info @@ -327,7 +332,7 @@ def import_header_only(self): """ Import only the header lines, without reading the actual data.""" - with open(self.filename, 'rb') as f: + with open(self.file_path, 'rb') as f: self.read_header(f) @staticmethod @@ -452,15 +457,15 @@ 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) - self.raw_info[current_file.filename] = current_file.raw_info - self.durations[current_file.filename] = current_file.duration() + self.raw_info[current_file.file_path] = current_file.raw_info + self.durations[current_file.file_path] = current_file.duration() file_laser_shots = [] self._create_or_append_channel(current_file) self.laser_shots.append(file_laser_shots) - self.files.append(current_file.filename) + self.files.append(current_file.file_path) def _create_or_append_channel(self, current_file): @@ -600,5 +605,8 @@ for key in keys: channel = self.channels[key] print("{0:<3} {1:<10} {2:<4} {3:<10} {4:<5}".format(channel.name, channel.wavelength, - channel.analog_photon_string, channel.resolution, - channel.points)) \ No newline at end of file + channel.analog_photon_string, channel.resolution, + channel.points)) + +class LicelDivaLidarMeasurement(DivaMixin, LicelLidarMeasurement): + pass \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/scripts/licel2tc.py Tue Sep 25 13:30:45 2018 +0300 @@ -0,0 +1,290 @@ +""" Command line tool to convert Licel binary files to EARLINET telecover files. +""" +import argparse +import glob +import logging +import os +import sys + +import numpy as np +from matplotlib import pyplot as plt +import yaml + +from lidar_processing import pre_processing + +from ..licel import LicelFile +from ..__init__ import __version__ + +logger = logging.getLogger(__name__) + + +class TelecoverDataset(object): + def __init__(self, filenames, settings_path, licel_timezone='UTC'): + """ Create a telecover dataset object. + + Parameters + ---------- + filenames : list of str + List of candidate files to include in the telecover. + settings_path : str + Path to the YAML settings file. + licel_timezone : str + The timezone of dates found in the Licel files. Should match the available + timezones in the TZ database. + """ + self.settings_path = settings_path + self.filenames = filenames + self.licel_timezone = licel_timezone + + self.settings = self.read_settings(settings_path) + self.file_sites = {} + self.read_files() + + self.check_files_ok() + + @staticmethod + def read_settings(settings_path): + """ Read the configuration file. + + The file should be in YAML syntax.""" + + if not os.path.isfile(settings_path): + raise IOError("Wrong path for settings file (%s)" % settings_path) + + with open(settings_path) as yaml_file: + try: + settings = yaml.safe_load(yaml_file) + logger.debug("Read settings file(%s)" % settings_path) + except: + raise IOError("Could not parse YAML file (%s)" % settings_path) + + return settings + + def read_files(self): + for filename in self.filenames: + # Load only header + current_file = LicelFile(filename, use_id_as_name=True, licel_timezone=self.licel_timezone, import_now=False) + current_file.import_header_only() + + # Populate the dictionary linking sites with files. + site = current_file.site + + if site in self.settings['sectors'].values(): + current_file.import_file() # Import full data + if current_file.site in self.file_sites.keys(): + + self.file_sites[site].append(current_file) + else: + self.file_sites[site] = [current_file, ] + + def check_files_ok(self): + """ Check if the available files are enough to create the ASCII files. """ + # Check that all sectors are available. + sectors = self.settings['sectors'] + for sector_name, site_name in sectors.items(): + if site_name not in self.file_sites.keys(): + raise ValueError('Site name {0} corresponding to sector {1} not found in any file.'.format(site_name, + sector_name)) + + len_sector = len(self.file_sites[site_name]) + if len_sector > 1: + logger.info('More than one files ({0}) found in sector {1} (site name {2})'.format(len_sector, + sector_name, + site_name)) + + # Check that all files have the correct channels + channels = self.settings['channels'].keys() + for site_name in sectors.values(): + for current_file in self.file_sites[site_name]: + for channel in channels: + if channel not in current_file.channels.keys(): + raise ValueError( + 'File {0} does not contain the required channel {1}.'.format(current_file.file_name, + channel)) + + def save_ascii(self, output_dir='.'): + """ Save the files in the appropriate format. + + Parameters + --------- + output_dir : str + The output directory. + + Returns + ------- + file_paths : list of str + A list of output paths for the saved files. + """ + + date, z = self._get_file_metadata() + + file_paths = [] + + logger.info('Creating {0} files.'.format(len(self.settings['channels']))) + + for channel_id, channel_settings in self.settings['channels'].items(): + + output_filename = '{call_sign}_tc_{date}_{name}.dat'.format(call_sign=self.settings['call_sign'], + date=date.strftime('%Y%m%dT%H'), + name=channel_settings['short_name']) + + output_path = os.path.join(output_dir, output_filename) + + logger.info("Output file path: {0}".format(output_path)) + + header_txt = ("{s[call_sign]} ({s[site]})\r\n" + "{s[system_name]}\r\n" + "{name}\r\n" + "{date}\r\n" + "range, {sectors}").format(s=self.settings, + name=channel_settings['name'], + date=date.strftime('%d.%m.%Y, %H'), + sectors=', '.join(self.settings['sector_order'])) + + # Get altitude data + sectors = self.settings['sector_order'] + + bin_max = self.settings['bins_to_keep'] + + data = [z, ] + for sector in sectors: + site_name = self.settings['sectors'][sector] + files = self.file_sites[site_name] + channel_data = [self.pre_process_channel(f.channels[channel_id]) for f in files] + data.append(np.mean(channel_data, axis=0)) + + data = np.array(data).T[:bin_max, :] + + np.savetxt(output_path, data, fmt='%.5e', delimiter=',', newline='\r\n', header=header_txt, comments='') + + logger.info("File saved.") + + file_paths.append(output_path) + + return file_paths + + def _get_file_metadata(self): + """ Get date string and altitude array that describes the measurement. + + Currently, use the date of the first file in the North direction. + """ + first_sector = self.settings['sector_order'][0] # Could be N, or NO, or ... + first_site = self.settings['sectors'][first_sector] + first_file = self.file_sites[first_site][0] + + channel_id = self.settings['channels'].keys()[0] + channel = first_file.channels[channel_id] + + return first_file.start_time, channel.z + + def pre_process_channel(self, channel): + """ Pre-process channel according to the settings. + + The role of this method is to choose if the channel is photon or analog and call the approparite + method. + """ + settings = self.settings['channels'][channel.channel_name] + + trigger_delay = settings.get('trigger_delay', 0) + background_min = settings['background_min'] + background_max = settings['background_max'] + + if channel.is_analog: + data_rc = self.pre_process_analog(channel, trigger_delay, background_min, background_max) + else: + dead_time = settings['dead_time'] + data_rc = self.pre_process_photon(channel, dead_time, trigger_delay, background_min, background_max) + + return data_rc + + @staticmethod + def pre_process_photon(channel, dead_time, trigger_delay, background_min, background_max): + """ Pre-process photon counting signals""" + interval_ns = channel.bin_width * 1e9 * channel.number_of_shots # Interval in ns, assuming constant laser shots + + data_dt = pre_processing.correct_dead_time_nonparalyzable(channel.data, interval_ns, dead_time) + data_bs, background_mean, background_std = pre_processing.subtract_background(data_dt, background_min, + background_max) + data_tc = pre_processing.correct_trigger_delay_bins(data_bs, trigger_delay) + data_rc = pre_processing.apply_range_correction(data_tc, channel.z) + + return data_rc + + @staticmethod + def pre_process_analog(channel, trigger_delay, background_min, background_max): + + data_bs, background_mean, background_std = pre_processing.subtract_background(channel.data, + background_min, + background_max) + data_tc = pre_processing.correct_trigger_delay_bins(data_bs, trigger_delay) + data_rc = pre_processing.apply_range_correction(data_tc, channel.z) + + return data_rc + + +def convert_to_telecover(file_paths, settings_path, licel_timezone, output_dir): + """ Convert files to SCC. """ + + try: + dataset = TelecoverDataset(file_paths, settings_path, licel_timezone) + except Exception as err: + logger.error(err) + sys.exit(1) + + try: + dataset.save_ascii(output_dir) + except Exception as err: + logger.error(err) + sys.exit(2) + + +def main(): + # Define the command line argument + parser = argparse.ArgumentParser( + description="A program to convert Licel binary files to EARLIENT telecover ASCII format") + parser.add_argument("settings_file", help="The path to a parameter YAML.") + parser.add_argument("files", + help="Location of licel files. Use relative path and filename wildcards. (default './*.*')", + default="./*.*") + parser.add_argument('-o', '--output', help='Output directory.', default='.') + + parser.add_argument('--licel_timezone', help="String describing the timezone according to the tz database.", + default="UTC", dest="licel_timezone", + ) + # Verbosity settings from http://stackoverflow.com/a/20663028 + parser.add_argument('-d', '--debug', help="Print dubuging information.", action="store_const", + dest="loglevel", const=logging.DEBUG, default=logging.INFO, + ) + parser.add_argument('-s', '--silent', help="Show only warning and error messages.", action="store_const", + dest="loglevel", const=logging.WARNING + ) + parser.add_argument('--version', help="Show current version.", action='store_true') + + args = parser.parse_args() + + # Get the logger with the appropriate level + logging.basicConfig(format='%(levelname)s: %(message)s', level=args.loglevel) + logger = logging.getLogger(__name__) + + # Check for version + if args.version: + print("Version: %s" % __version__) + sys.exit(0) + + # Get a list of files to process + logger.debug("Search path: {0}".format(os.path.expanduser(args.files))) + files = glob.glob(os.path.expanduser(args.files)) + + # If not files found, exit + if len(files) == 0: + logger.error("No files found when searching for %s." % args.files) + sys.exit(1) + + # If everything OK, proceed + logger.info("Found {0} files matching {1}".format(len(files), args.files)) + + convert_to_telecover(files, args.settings_file, args.licel_timezone, args.output) + + +if __name__ == "__main__": + main()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/config/telecover_settings.yaml Tue Sep 25 13:30:45 2018 +0300 @@ -0,0 +1,26 @@ +site: Munich # e.g. Munich +call_sign: mu # according to EARLINET +system_name: POLIS # according to HOI +channels: # Link Licel channel name with HOI channel name + BT0: + name: 355 elastic, parallel, analog # For file header + short_name: 355_epa # For file name + trigger_delay: -9 # Positive: trigger delay, Negative: bin shift + background_min: 14000 # Minimum bin for background subtraction + background_max: 15000 # Maximum bin for background subtraction + BC1: + name: 355 elastic, parallel, photon-counting + short_name: 355_epp + trigger_delay: -9 # Positive: trigger delay, Negative: bin shift + background_min: 14000 # Minimum bin for background subtraction + background_max: 15000 # Maximum bin for background subtraction + dead_time: 3.8 # Dead time in ns +sectors: # Link Licel "site" value (as recorded in the file header) with telecover sector + N: 'NO' # These are key: value pairs. The key is the name of the sector, the value is the "site" name in the file + E: EA # 'NO' is a reserved word and evaluates to false. Thus it has to be written in quotes + W: WE + S: SO + N2: NO2 + D: DARK +sector_order: [N, E, W, S, N2, D] # Order for recording +bins_to_keep: 1500 # Number of bins to keep in each profile \ No newline at end of file
--- a/setup.py Thu Sep 06 16:51:13 2018 +0300 +++ b/setup.py Tue Sep 25 13:30:45 2018 +0300 @@ -58,6 +58,7 @@ ], entry_points={ 'console_scripts': ['licel2scc = atmospheric_lidar.scripts.licel2scc:main', - 'licel2scc-depol = atmospheric_lidar.scripts.licel2scc_depol:main'], + 'licel2scc-depol = atmospheric_lidar.scripts.licel2scc_depol:main', + 'licel2tc = atmospheric_lidar.scripts.licel2tc:main'], }, )
