Mercurial > public > atmospheric_lidar / changeset
changeset
Update to integrate automatic cloud masking in licel convertion.
Fri, 23 Feb 2018 14:45:36 +0200
- author
- Iannis <i.binietoglou@impworks.gr>
- date
- Fri, 23 Feb 2018 14:45:36 +0200
- changeset 128
- 168da625489b
- parent 127
- 6051a9bae77b
- child 129
- 6f902a45b83d
Update to integrate automatic cloud masking in licel convertion.
--- a/atmospheric_lidar/generic.py Fri Feb 23 08:58:50 2018 +0200 +++ b/atmospheric_lidar/generic.py Fri Feb 23 14:45:36 2018 +0200 @@ -443,6 +443,7 @@ # Write the custom subclass variables: for variable in self._get_custom_variables(channel_names): + logger.debug("Saving variable {0}".format(variable['name'])) temp_v = f.createVariable(variable["name"], variable["type"], variable["dimensions"]) temp_v[:] = variable["values"]
--- a/atmospheric_lidar/licel.py Fri Feb 23 08:58:50 2018 +0200 +++ b/atmospheric_lidar/licel.py Fri Feb 23 14:45:36 2018 +0200 @@ -12,6 +12,136 @@ c = 299792458.0 # Speed of light +class LicelFileChannel: + """ A class representing a single channel found in a single Licel file.""" + + def __init__(self, raw_info, raw_data, duration, use_id_as_name=False): + """ + This is run when creating a new object. + + Parameters + ---------- + raw_info : dict + A dictionary containing raw channel information. + raw_data : dict + An array with raw channel data. + duration : float + Duration of the file, in seconds + use_id_as_name : bool + If True, the transient digitizer name (e.g. BT0) is used as a channel + name. If False, a more descriptive name is used (e.g. '01064.o_an'). + """ + self.raw_info = raw_info + self.raw_data = raw_data + self.duration = duration + self.use_id_as_name = use_id_as_name + self.adcbits = int(raw_info['ADCbits']) + self.active = int(raw_info['Active']) + self.analog_photon = raw_info['AnalogPhoton'] + self.bin_width = float(raw_info['BinW']) + self.data_points = int(raw_info['DataPoints']) + + self.hv = float(raw_info['HV']) + self.id = raw_info['ID'] + self.laser_user = int(raw_info['LaserUsed']) + self.number_of_shots = int(raw_info['NShots']) + self.wavelength_str = raw_info['Wavelength'] + + if self.is_analog: + self.discriminator = float(raw_info['Discriminator']) * 1000 # Analog range in mV + else: + self.discriminator = float(raw_info['Discriminator']) + + @property + def wavelength(self): + """ Property describing the nominal wavelength of the channel. + + Returns + ------- + : int or None + The integer value describing the wavelength. If no raw_info have been provided, + returns None. + """ + wavelength = self.wavelength_str.split('.')[0] + return int(wavelength) + + @property + def channel_name(self): + """ + Construct the channel name adding analog photon info to avoid duplicates + + If use_id_as_name is True, the channel name will be the transient digitizer ID (e.g. BT01). + This could be useful if the lidar system has multiple telescopes, so the descriptive name is + not unique. + + Returns + ------- + channel_name : str + The channel name + """ + if self.use_id_as_name: + channel_name = self.id + else: + acquisition_type = self.analog_photon_string + channel_name = "%s_%s" % (self.wavelength_str, acquisition_type) + return channel_name + + @property + def analog_photon_string(self): + """ Convert the analog/photon flag found in the Licel file to a proper sting. + + Returns + ------- + string : str + 'an' or 'ph' string, for analog or photon-counting channel respectively. + """ + if self.analog_photon == '0': + string = 'an' + else: + string = 'ph' + return string + + def calculate_physical(self): + """ Calculate physically-meaningful data from raw channel data: + + * In case of analog signals, the data are converted to mV. + * In case of photon counting signals, data are stored as number of photons. + + In addition, some ancillary variables are also calculated (z, dz, number_of_bins). + """ + data = self.raw_data + + norm = data / float(self.number_of_shots) + dz = float(self.raw_info['BinW']) + + if self.raw_info['AnalogPhoton'] == '0': + # If the channel is in analog mode + ADCrange = self.discriminator # Discriminator value already in mV + channel_data = norm * ADCrange / ((2 ** self.adcbits) - 1) # TODO: Check this. Agrees with Licel docs, but differs from their LabView code. + + # print ADCb, ADCRange,cdata,norm + else: + # If the channel is in photoncounting mode + # Frequency deduced from range resolution! (is this ok?) + # c = 300 # The result will be in MHZ + # SR = c/(2*dz) # To account for pulse folding + # channel_data = norm*SR + # CHANGE: + # For the SCC the data are needed in photons + channel_data = norm * self.number_of_shots + # print res,c,cdata,norm + + # Calculate Z + + self.z = np.array([dz * bin_number + dz / 2.0 for bin_number in range(self.data_points)]) + self.dz = dz + self.data = channel_data + + @property + def is_analog(self): + return self.analog_photon == '0' + + class LicelFile: """ A class representing a single binary Licel file. """ @@ -26,7 +156,7 @@ def __init__(self, file_path, use_id_as_name=False, licel_timezone="UTC"): """ This is run when creating a new object. - + Parameters ---------- file_path : str @@ -36,7 +166,7 @@ name. If False, a more descriptive name is used (e.g. '01064.o_an'). licel_timezone : str The timezone of dates found in the Licel files. Should match the available - timezones in the TZ database. + timezones in the TZ database. """ self.filename = file_path self.use_id_as_name = use_id_as_name @@ -53,7 +183,7 @@ def _import_file(self, file_path): """ Read the content of the Licel file. - + Parameters ---------- file_path : str @@ -90,8 +220,8 @@ self.channels = channels def read_header(self, f): - """ Read the header of an open Licel file. - + """ Read the header of an open Licel file. + Parameters ---------- f : file-like object @@ -163,8 +293,8 @@ return raw_dict def duration(self): - """ Return the duration of the file. - + """ Return the duration of the file. + Returns ------- : float @@ -187,135 +317,75 @@ return combined -class LicelFileChannel: - """ A class representing a single channel found in a single Licel file.""" - - def __init__(self, raw_info, raw_data, duration, use_id_as_name=False): - """ - This is run when creating a new object. +class LicelChannel(LidarChannel): - Parameters - ---------- - raw_info : dict - A dictionary containing raw channel information. - raw_data : dict - An array with raw channel data. - duration : float - Duration of the file, in seconds - use_id_as_name : bool - If True, the transient digitizer name (e.g. BT0) is used as a channel - name. If False, a more descriptive name is used (e.g. '01064.o_an'). - """ - self.raw_info = raw_info - self.raw_data = raw_data - self.duration = duration - self.use_id_as_name = use_id_as_name - self.adcbits = int(raw_info['ADCbits']) - self.active = int(raw_info['Active']) - self.analog_photon = raw_info['AnalogPhoton'] - self.bin_width = float(raw_info['BinW']) - self.data_points = int(raw_info['DataPoints']) + def __init__(self): + self.name = None + self.resolution = None + self.points = None + self.wavelength = None + self.laser_used = None - self.hv = float(raw_info['HV']) - self.id = raw_info['ID'] - self.laser_user = int(raw_info['LaserUsed']) - self.number_of_shots = int(raw_info['NShots']) - self.wavelength_str = raw_info['Wavelength'] - - if self.is_analog: - self.discriminator = float(raw_info['Discriminator']) * 1000 # Analog range in mV - else: - self.discriminator = float(raw_info['Discriminator']) + self.rc = [] + self.raw_info = [] + self.laser_shots = [] + self.duration = [] + self.number_of_shots = [] + self.discriminator = [] + self.hv = [] + self.data = {} - @property - def wavelength(self): - """ Property describing the nominal wavelength of the channel. - - Returns - ------- - : int or None - The integer value describing the wavelength. If no raw_info have been provided, - returns None. - """ - wavelength = self.wavelength_str.split('.')[0] - return int(wavelength) + def append_file(self, file_start_time, file_channel): + """ Append file to the current object """ - @property - def channel_name(self): - """ - Construct the channel name adding analog photon info to avoid duplicates + self._assign_properties(file_channel) + + self.binwidth = self.resolution * 2. / c # in seconds + self.z = file_channel.z - If use_id_as_name is True, the channel name will be the transient digitizer ID (e.g. BT01). - This could be useful if the lidar system has multiple telescopes, so the descriptive name is - not unique. - - Returns - ------- - channel_name : str - The channel name - """ - if self.use_id_as_name: - channel_name = self.id - else: - acquisition_type = self.analog_photon_string - channel_name = "%s_%s" % (self.wavelength_str, acquisition_type) - return channel_name - - @property - def analog_photon_string(self): - """ Convert the analog/photon flag found in the Licel file to a proper sting. - - Returns - ------- - string : str - 'an' or 'ph' string, for analog or photon-counting channel respectively. - """ - if self.analog_photon == '0': - string = 'an' - else: - string = 'ph' - return string + self.data[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) - def calculate_physical(self): - """ Calculate physically-meaningful data from raw channel data: - - * In case of analog signals, the data are converted to mV. - * In case of photon counting signals, data are stored as number of photons. - - In addition, some ancillary variables are also calculated (z, dz, number_of_bins). - """ - data = self.raw_data - - norm = data / float(self.number_of_shots) - dz = float(self.raw_info['BinW']) - - if self.raw_info['AnalogPhoton'] == '0': - # If the channel is in analog mode - ADCrange = self.discriminator # Discriminator value already in mV - channel_data = norm * ADCrange / ((2 ** self.adcbits) - 1) # TODO: Check this. Agrees with Licel docs, but differs from their LabView code. + def _assign_properties(self, file_channel): + self._assign_unique_property('name', file_channel.channel_name) + self._assign_unique_property('resolution', file_channel.dz) + self._assign_unique_property('wavelength', file_channel.wavelength) + self._assign_unique_property('points', file_channel.data_points) + self._assign_unique_property('adcbits', file_channel.adcbits) + self._assign_unique_property('active', file_channel.active) + self._assign_unique_property('laser_user', file_channel.laser_user) + self._assign_unique_property('adcbints', file_channel.adcbits) + self._assign_unique_property('analog_photon_string', file_channel.analog_photon_string) - # print ADCb, ADCRange,cdata,norm + def _assign_unique_property(self, property_name, value): + + current_value = getattr(self, property_name, None) + + if current_value is None: + setattr(self, property_name, value) else: - # If the channel is in photoncounting mode - # Frequency deduced from range resolution! (is this ok?) - # c = 300 # The result will be in MHZ - # SR = c/(2*dz) # To account for pulse folding - # channel_data = norm*SR - # CHANGE: - # For the SCC the data are needed in photons - channel_data = norm * self.number_of_shots - # print res,c,cdata,norm - - # Calculate Z - - self.z = np.array([dz * bin_number + dz / 2.0 for bin_number in range(self.data_points)]) - self.dz = dz - self.data = channel_data + if current_value != value: + raise ValueError('Cannot combine channels with different values of {0}.'.format(property_name)) @property def is_analog(self): - return self.analog_photon == '0' + return self.analog_photon_string == 'an' + + @property + def is_photon_counting(self): + return self.analog_photon_string == 'ph' + def __unicode__(self): + return "<Licel channel: %s>" % self.name + + def __str__(self): + return unicode(self).encode('utf-8') + class LicelLidarMeasurement(BaseLidarMeasurement): @@ -465,73 +535,3 @@ This requires changes in generic.py """ raise NotImplementedError("Subsetting by time, not yet implemented for Licel files.") - - -class LicelChannel(LidarChannel): - - def __init__(self): - self.name = None - self.resolution = None - self.points = None - self.wavelength = None - self.laser_used = None - - self.rc = [] - self.raw_info = [] - self.laser_shots = [] - self.duration = [] - self.number_of_shots = [] - self.discriminator = [] - self.hv = [] - self.data = {} - - def append_file(self, file_start_time, file_channel): - """ Append file to the current object """ - - self._assign_properties(file_channel) - - self.binwidth = self.resolution * 2. / c # in seconds - self.z = file_channel.z - - self.data[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) - - def _assign_properties(self, file_channel): - self._assign_unique_property('name', file_channel.channel_name) - self._assign_unique_property('resolution', file_channel.dz) - self._assign_unique_property('wavelength', file_channel.wavelength) - self._assign_unique_property('points', file_channel.data_points) - self._assign_unique_property('adcbits', file_channel.adcbits) - self._assign_unique_property('active', file_channel.active) - self._assign_unique_property('laser_user', file_channel.laser_user) - self._assign_unique_property('adcbints', file_channel.adcbits) - self._assign_unique_property('analog_photon_string', file_channel.analog_photon_string) - - def _assign_unique_property(self, property_name, value): - - current_value = getattr(self, property_name, None) - - if current_value is None: - setattr(self, property_name, value) - else: - if current_value != value: - raise ValueError('Cannot combine channels with different values of {0}.'.format(property_name)) - - @property - def is_analog(self): - return self.analog_photon_string == 'an' - - @property - def is_photon_counting(self): - return self.analog_photon_string == 'ph' - - def __unicode__(self): - return "<Licel channel: %s>" % self.name - - def __str__(self): - return unicode(self).encode('utf-8')
--- a/atmospheric_lidar/scripts/licel2scc.py Fri Feb 23 08:58:50 2018 +0200 +++ b/atmospheric_lidar/scripts/licel2scc.py Fri Feb 23 14:45:36 2018 +0200 @@ -7,6 +7,9 @@ import os import sys +from matplotlib import pyplot as plt +import yaml + from ..licel import LicelLidarMeasurement from ..__init__ import __version__ @@ -74,7 +77,27 @@ return settings -def get_cloud_free_files(CustomLidarMeasurement, files, args): +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(CustomLidarMeasurement, files, settings): logger.warning("Starting cloud mask procedure. This is an experimental feature.") try: @@ -84,15 +107,30 @@ sys.exit(1) measurement = CustomLidarMeasurement(files) - channel = measurement.channels[args.cloudmask_channel] + channel = measurement.channels[settings['channel']] cloud_mask = cloudmask.CloudMaskRaw(channel) - idxs = cloud_mask.cloud_free_periods(args.cloudfree_period, args.cloud_search_height) + + 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']: + output_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')) + output_path = os.path.join(settings['plot_directory'], output_filename) + fig, _ = cloud_mask.plot_cloudfree(idxs) + + plt.savefig(output_path) + file_list = [] for idx_min, idx_max in idxs: current_files = measurement.files[idx_min:idx_max] @@ -168,16 +206,9 @@ parser.add_argument('--licel_timezone', help="String describing the timezone according to the tz database.", default="UTC", dest="licel_timezone", ) - parser.add_argument('--cloudmask', help="Experimental feature to automatically cloud mask measurements", - default=False, action='store_true', - ) - parser.add_argument('--cloudmask_channel', help="Name of channel to apply the cloud mask.") - parser.add_argument('--cloudfree_period', type=float, help="Duration (in min) of cloud-free periods", - default="30", - ) - parser.add_argument('--cloud_search_height', type=float, help="Maximum altitude (in m) to check for clouds.", - default="12000", - ) + parser.add_argument('--cloudmask_settings', help="Experimental feature to automatically cloud mask measurements", + default="") + # 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, @@ -213,8 +244,10 @@ CustomLidarMeasurement = create_custom_class(args.parameter_file, args.id_as_name, args.temperature, args.pressure, args.licel_timezone) - if args.cloudmask: - file_lists = get_cloud_free_files(CustomLidarMeasurement, files, args) + if args.cloudmask_settings: + cloudmask_settings = read_cloudmask_settings_file(args.cloudmask_settings) + + file_lists = get_cloud_free_files(CustomLidarMeasurement, files, cloudmask_settings) for n, files in enumerate(file_lists): measurement_id, measurement_no = get_corrected_measurement_id(args, n)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/config/cloudmask_config_sample.yaml Fri Feb 23 14:45:36 2018 +0200 @@ -0,0 +1,7 @@ +channel: "01064.o_an" # Identified of the channel to use for cloud masking +cloudfree_period_min: 60 # Minimum duration of cloud-free period to convert (in minutes) +file_duration_max: 120 # Maximum duration of file to create (in minutes) +max_cloud_height: 8000 # Maximum altitude (in m) to check for clouds. +plot: True # Create plots of cloud free periods +plot_directory: ./ # Output directory for creating plots +
