Mercurial > public > atmospheric_lidar / changeset
changeset
First ideas for licel-to-telecover converter.
Fri, 14 Sep 2018 16:37:20 +0300
- author
- Iannis <i.binietoglou@impworks.gr>
- date
- Fri, 14 Sep 2018 16:37:20 +0300
- changeset 156
- 1e18b2a416ad
- parent 155
- 4a596849c721
- child 157
- 03b470b0a05f
First ideas for licel-to-telecover converter.
| atmospheric_lidar/licel.py | file | annotate | diff | comparison | revisions | |
| atmospheric_lidar/scripts/licel2qa.py | file | annotate | diff | comparison | revisions | |
| atmospheric_lidar/scripts/licel2tc.py | file | annotate | diff | comparison | revisions | |
| config/telecover_settings.yaml | file | annotate | diff | comparison | revisions |
--- a/atmospheric_lidar/licel.py Thu Sep 13 14:15:34 2018 +0300 +++ b/atmospheric_lidar/licel.py Fri Sep 14 16:37:20 2018 +0300 @@ -1,6 +1,7 @@ import datetime import logging import copy +import os import numpy as np import pytz @@ -172,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 @@ -187,7 +190,7 @@ channels = {} photodiodes = {} - with open(self.filename, 'rb') as f: + with open(self.file_path, 'rb') as f: self.read_header(f) @@ -201,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(), @@ -296,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 @@ -328,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 @@ -453,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):
--- a/atmospheric_lidar/scripts/licel2qa.py Thu Sep 13 14:15:34 2018 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,280 +0,0 @@ -""" Command line tool to convert Licel binary files to EARLINET telecover files. -""" -import argparse -import glob -import importlib -import logging -import os -import sys - -from matplotlib import pyplot as plt -import yaml - -from ..licel import LicelLidarMeasurement -from ..__init__ import __version__ - -logger = logging.getLogger(__name__) - - -class TelecoverMeasurement(LicelLidarMeasurement): - - def __init__(self, file_list, use_id_as_name, licel_timezone, telecover_settings): - self.telecover_settings = telecover_settings - super(TelecoverMeasurement, self).__init__(file_list, use_id_as_name, licel_timezone=licel_timezone) - - -def create_custom_class(custom_netcdf_parameter_path, use_id_as_name=False, temperature=25., pressure=1020., - licel_timezone='UTC'): - """ This funtion creates a custom LicelLidarMeasurement subclass, - based on the input provided by the users. - - Parameters - ---------- - custom_netcdf_parameter_path : str - The path to the custom channels parameters. - use_id_as_name : bool - Defines if channels names are descriptive or transient digitizer IDs. - temperature : float - The ground temperature in degrees C (default 25.0). - pressure : float - The ground pressure in hPa (default: 1020.0). - licel_timezone : str - String describing the timezone according to the tz database. - - Returns - ------- - CustomLidarMeasurement: - A custom sub-class of LicelLidarMeasurement - """ - logger.debug('Reading parameter files: %s' % custom_netcdf_parameter_path) - custom_netcdf_parameters = read_settings_file(custom_netcdf_parameter_path) - - class CustomLidarMeasurement(LicelLidarMeasurement): - extra_netcdf_parameters = custom_netcdf_parameters - - def __init__(self, file_list=None): - super(CustomLidarMeasurement, self).__init__(file_list, use_id_as_name, licel_timezone=licel_timezone) - - def set_PT(self): - ''' Sets the pressure and temperature at station level. This is used if molecular_calc parameter is - set to 0 (i.e. use US Standard atmosphere). - - The results are stored in the info dictionary. - ''' - - self.info['Temperature'] = temperature - self.info['Pressure'] = pressure - - return CustomLidarMeasurement - - -def read_settings_file(settings_path): - """ Read the settings file. - - The file should contain python code.""" - if not os.path.isfile(settings_path): - logging.error("The provided settings path does not correspond to a file.") - sys.exit(1) - - dirname, basename = os.path.split(settings_path) - sys.path.append(dirname) - - module_name, _ = os.path.splitext(basename) - settings = importlib.import_module(module_name) - return settings - - -def read_cloudmask_settings_file(settings_file_path): - """ Read the configuration file. - - The file should be in YAML syntax.""" - - if not os.path.isfile(settings_file_path): - logging.error("Wrong path for cloudmask settings file (%s)" % settings_file_path) - sys.exit(1) - - with open(settings_file_path) as yaml_file: - try: - settings = yaml.load(yaml_file) - logging.debug("Read cloudmask settings file(%s)" % settings_file_path) - except: - logging.error("Could not parse YAML file (%s)" % settings_file_path) - sys.exit(1) - - return settings - - -def get_cloud_free_files(LidarMeasurementClass, files, settings): - """ Find cloud free periods in the given files. - - Depending on the provided settings, it could create plots of cloud mask and - selected cloud-free periods. - - Parameters - ---------- - LidarMeasurementClass : class - Class used to read the files. - files : list - A list of raw licel file paths. - settings : dict - A dictionary of cloud masking settings. - - Returns - ------- - file_list : list of lists - A list of lists containing paths to cloud-free files. - """ - logger.warning("Starting cloud mask procedure. This is an experimental feature.") - - try: - from cloudmask import cloudmask # Import here until we setup a proper installation procedure - except ImportError: - logger.error("Cloud mask module could not be loaded. Please install manually.") - sys.exit(1) - - measurement = LidarMeasurementClass(files) - channel = measurement.channels[settings['channel']] - cloud_mask = cloudmask.CloudMaskRaw(channel) - - idxs = cloud_mask.cloud_free_periods(settings['cloudfree_period_min'], - settings['file_duration_max'], - settings['max_cloud_height']) - - logger.debug('Cloud free indices: {0}'.format(idxs)) - - if len(idxs) == 0: # If no cloud-free period found - logger.info('No cloud free period found. Nothing converted.') - sys.exit(1) - - logger.info("{0} cloud free period(s) found.".format(len(idxs))) - - if settings['plot']: - # Plot cloud free periods - cloudfree_filename = "cloudfree_{0}_{1}_{2}.png".format(channel.wavelength, - channel.start_time.strftime('%Y%m%d_%H%M%S'), - channel.stop_time.strftime('%Y%m%d_%H%M%S')) - cloudfree_path = os.path.join(settings['plot_directory'], cloudfree_filename) - fig, _ = cloud_mask.plot_cloudfree(idxs) - - plt.savefig(cloudfree_path) - plt.close() - - # Plot cloud mask - cloudmask_filename = "cloudmask_{0}_{1}_{2}.png".format(channel.wavelength, - channel.start_time.strftime('%Y%m%d_%H%M%S'), - channel.stop_time.strftime('%Y%m%d_%H%M%S')) - cloudmask_path = os.path.join(settings['plot_directory'], cloudmask_filename) - - fig, _ = cloud_mask.plot_mask() - - plt.savefig(cloudmask_path) - plt.close() - - file_list = [] - for idx_min, idx_max in idxs: - current_files = measurement.files[idx_min:idx_max] - file_list.append(current_files) - - return file_list - - -def get_corrected_measurement_id(args, n): - """ Correct the provided measurement id, in case of multiple cloud-free periods. """ - if args.measurement_id is not None: - order = int(args.measurement_id[-2:]) - new_no = order + n - measurement_id = args.measurement_id[:-2] + str(new_no) - measurement_no = args.measurement_number # The same - else: - measurement_no = str(int(args.measurement_number) + n).zfill(2) - measurement_id = None - - return measurement_id, measurement_no - - -def convert_to_scc(CustomLidarMeasurement, files, dark_pattern, measurement_id, measurement_number): - """ Convert files to SCC. """ - measurement = CustomLidarMeasurement(files) - # Get a list of files containing dark measurements - if dark_pattern != "": - dark_files = glob.glob(dark_pattern) - - if dark_files: - logger.debug("Using %s as dark measurements files!" % ', '.join(dark_files)) - measurement.dark_measurement = CustomLidarMeasurement(dark_files) - else: - logger.warning( - 'No dark measurement files found when searching for %s. Will not use any dark measurements.' % dark_pattern) - try: - measurement = measurement.subset_by_scc_channels() - except ValueError as err: - logger.error(err) - sys.exit(1) - - # Save the netcdf - logger.info("Saving netcdf") - measurement.set_measurement_id(measurement_id, measurement_number) - measurement.save_as_SCC_netcdf() - logger.info("Created file %s" % measurement.scc_filename) - - -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("-i", '--id_as_name', - help="Use transient digitizer ids as channel names, instead of descriptive names", - action="store_true") - parser.add_argument("-t", "--temperature", type=float, - help="The temperature (in C) at lidar level, required if using US Standard atmosphere", - default="25") - parser.add_argument("-p", "--pressure", type=float, - help="The pressure (in hPa) at lidar level, required if using US Standard atmosphere", - default="1020") - parser.add_argument('-D', '--dark_measurements', - help="Location of files containing dark measurements. Use relative path and filename wildcars, see 'files' parameter for example.", - default="", dest="dark_files" - ) - 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 - files = glob.glob(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), os.path.abspath(args.files))) - - CustomLidarMeasurement = create_custom_class(args.parameter_file, args.id_as_name, args.temperature, - args.pressure, args.licel_timezone) - - convert_to_telecover(CustomLidarMeasurement, files, args.dark_files, args.measurement_id, args.measurement_number)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/scripts/licel2tc.py Fri Sep 14 16:37:20 2018 +0300 @@ -0,0 +1,293 @@ +""" Command line tool to convert Licel binary files to EARLINET telecover files. +""" +import argparse +import glob +import importlib +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: + current_file = LicelFile(filename, use_id_as_name=True, licel_timezone=self.licel_timezone) + + # Populate the dictionary linking sites with files. + site = current_file.site + 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 = [] + 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) + + 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='') + 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_dealy', 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(files, dark_pattern, measurement_id, measurement_number): + """ Convert files to SCC. """ + measurement = CustomLidarMeasurement(files) + # Get a list of files containing dark measurements + if dark_pattern != "": + dark_files = glob.glob(dark_pattern) + + if dark_files: + logger.debug("Using %s as dark measurements files!" % ', '.join(dark_files)) + measurement.dark_measurement = CustomLidarMeasurement(dark_files) + else: + logger.warning( + 'No dark measurement files found when searching for %s. Will not use any dark measurements.' % dark_pattern) + try: + measurement = measurement.subset_by_scc_channels() + except ValueError as err: + logger.error(err) + sys.exit(1) + + # Save the netcdf + logger.info("Saving netcdf") + measurement.set_measurement_id(measurement_id, measurement_number) + measurement.save_as_SCC_netcdf() + logger.info("Created file %s" % measurement.scc_filename) + + +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("-i", '--id_as_name', + help="Use transient digitizer ids as channel names, instead of descriptive names", + action="store_true") + parser.add_argument("-t", "--temperature", type=float, + help="The temperature (in C) at lidar level, required if using US Standard atmosphere", + default="25") + parser.add_argument("-p", "--pressure", type=float, + help="The pressure (in hPa) at lidar level, required if using US Standard atmosphere", + default="1020") + parser.add_argument('-D', '--dark_measurements', + help="Location of files containing dark measurements. Use relative path and filename wildcars, see 'files' parameter for example.", + default="", dest="dark_files" + ) + 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 + files = glob.glob(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), os.path.abspath(args.files))) + + logger.error('Nothing implemented yet!')
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/config/telecover_settings.yaml Fri Sep 14 16:37:20 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
