--- a/atmospheric_lidar/diva.py Thu Jan 04 14:47:36 2018 +0200
+++ b/atmospheric_lidar/diva.py Sat Jan 06 10:08:46 2018 +0200
@@ -4,20 +4,23 @@
separately not to interfere with normal development.
"""
import netCDF4 as netcdf
-import pyyaml
+import yaml
import datetime
import os
+import numpy as np
+
+import pytz
from .generic import BaseLidarMeasurement
class DivaOutput(BaseLidarMeasurement):
- def save_as_diva(self, output_path, parameter_file):
+ def save_as_diva_netcdf(self, output_path, parameter_file):
""" Save the current data in the 'draft' DIVA format. """
with open(parameter_file, 'r') as f:
- parameters = pyyaml.load(f)
+ parameters = yaml.load(f)
global_parameters = parameters['global_parameters'] # Shortcut
channels = parameters['channels']
@@ -96,28 +99,29 @@
for channel_name, channel_parameters in channels.iteritems():
if channel_name not in self.channels.keys():
- raise ValueError('Channle name not one of {0}: {1}'.format(self.channels.keys()), channel_name)
+ raise ValueError('Channel name not one of {0}: {1}'.format(self.channels.keys(), channel_name))
channel = self.channels[channel_name]
- group_name = "channel_{0}".format(channel_name) # Give channels groups a standard name
+ group_name = "channel_{0}".format(channel_name.replace('.', '_')) # Give channels groups a standard name
+
g = f.createGroup(group_name)
g.long_name = channel_parameters['long_name']
- g.detector_manufacturer = channel_parameters['detector_manifacturer'] # Optional
+ g.detector_manufacturer = channel_parameters['detector_manufacturer'] # Optional
g.detector_model = channel_parameters['detector_model']
g.daq_manufacturer = channel_parameters['daq_manufacturer']
g.daq_model = channel_parameters['daq_model']
# Dimensions
g.createDimension('profile', size=None) # Infinite dimension
- g.createDimension('range', len(self.range))
+ g.createDimension('range', len(channel.z))
# Variables
- name = g.createVariable('channel_id', dimensions=('name_strlen',))
+ name = g.createVariable('channel_id', 'c', dimensions=('name_strlen',))
name.cf_role = 'timeseries_id'
name.long_name = 'channel identification'
- laser_rep_rate = g.createVariable('laser_repetition_rate')
+ laser_rep_rate = g.createVariable('laser_repetition_rate', 'f4')
laser_rep_rate.long_name = 'nominal laser repetition rate'
laser_rep_rate.units = 'Hz'
@@ -166,44 +170,42 @@
polarizer_angle.units = 'degree'
polarizer_angle.comments = 'Optional'
- dead_time_model = g.createVariable('dead_time_model', datatype='b')
- dead_time_model.long_name = 'optimal dead time model of detection system'
- dead_time_model.flag_values = '0b 1b 2b'
- dead_time_model.flag_meanings = 'paralyzable non_paralyzable other'
+ if not channel.is_analog:
+ dead_time_model = g.createVariable('dead_time_model', datatype='b')
+ dead_time_model.long_name = 'optimal dead time model of detection system'
+ dead_time_model.flag_values = '0b 1b 2b'
+ dead_time_model.flag_meanings = 'paralyzable non_paralyzable other'
- dead_time = g.createVariable('dead_time', datatype='f8')
- dead_time.long_name = 'dead time value'
- dead_time.units = 'ns'
- dead_time.comment = "Manufacturer. Source of the value."
+ dead_time = g.createVariable('dead_time', datatype='f8')
+ dead_time.long_name = 'dead time value'
+ dead_time.units = 'ns'
+ dead_time.comment = "Manufacturer. Source of the value."
bin_length = g.createVariable('bin_length', datatype='f4')
bin_length.long_name = "time duration of each bin"
bin_length.units = 'ns'
- pulse_delay = g.createVariable('pulse_delay', datatype='f4')
- pulse_delay.long_name = "time difference between channel trigger and pulse emission"
- pulse_delay.units = 'ns'
- pulse_delay.comment = "Positive indicates pre-trigger channel. Negative a channle trigger delay."
+ if channel.is_analog:
+ adc_bits = g.createVariable('adc_bits', datatype='i4')
+ adc_bits.long_name = 'analog-to-digital converter bits'
+ adc_bits.coordinates = "time"
detector_voltage = g.createVariable('detector_voltage', datatype='f4', dimensions=('profile',), zlib=True)
detector_voltage.long_name = 'detector voltage'
detector_voltage.units = 'V'
detector_voltage.coordinates = "time"
- discriminator = g.createVariable('discriminator', datatype='f8', dimensions=('profiles',))
- discriminator.long_name = 'discriminator level'
- discriminator.units = ''
+ if channel.is_photon_counting:
+ discriminator = g.createVariable('discriminator', datatype='f8', dimensions=('profile',))
+ discriminator.long_name = 'discriminator level'
+ discriminator.units = ''
- adc_range = g.createVariable('adc_range', datatype='f4', dimensions=('profile',),
- zlib=True)
- adc_range.long_name = 'analog-to-digital converter range'
- adc_range.units = 'mV'
- adc_range.coordinates = "time"
-
- adc_bits = g.createVariable('adc_bits', datatype='i4', dimensions=('profile',),
- zlib=True)
- adc_bits.long_name = 'analog-to-digital converter bits'
- adc_bits.coordinates = "time"
+ if channel.is_analog:
+ adc_range = g.createVariable('adc_range', datatype='f4', dimensions=('profile',),
+ zlib=True)
+ adc_range.long_name = 'analog-to-digital converter range'
+ adc_range.units = 'mV'
+ adc_range.coordinates = "time"
pulses = g.createVariable('pulses', datatype='i4', dimensions=('profile',),
zlib=True)
@@ -214,10 +216,10 @@
nd_filter.long_name = "neutral density filter optical depth "
nd_filter.coordinates = "time"
- emission_delay = g.createVariable('emission_delay', datatype='f4')
- emission_delay.long_name = "pulse emission difference from channel trigger"
- emission_delay.units = 'ns'
- emission_delay.comments = 'Positive values for pre-trigger systems. Negative for trigger delay.'
+ trigger_delay = g.createVariable('trigger_delay', datatype='f4')
+ trigger_delay.long_name = "chanenl trigger difference from pulse emission"
+ trigger_delay.units = 'ns'
+ trigger_delay.comments = 'Negative values for pre-trigger systems.'
time = g.createVariable('time', datatype='f8', dimensions=('profile',),
zlib=True)
@@ -225,29 +227,36 @@
time.units = "seconds since 1970-01-01 00:00:00"
time.standard_name = "time"
time.bounds = "time_bnds"
- g.createVariable('time_bnds', datatype='f8', dimensions=('profile', 'nv'), zlib=True)
+ time_bounds = g.createVariable('time_bnds', datatype='f8', dimensions=('profile', 'nv'), zlib=True)
- bin_time = g.createVariable('bin_time', datatype='f4', dimensions=('range',),
- zlib=True) # Use name 'bin_range' to avoid conflict with built-in range
- bin_time.long_name = 'bin start time since trigger'
+ bin_time = g.createVariable('bin_time', datatype='f4', dimensions=('range',), zlib=True)
+ bin_time.long_name = 'bin start time since channel trigger'
bin_time.units = "ns"
- signal = g.createVariable('signal', datatype='i4', dimensions=('profile', 'range'),
+ if channel.is_analog:
+ signal_units = 'mV'
+ signal_datatype = 'f8'
+ else:
+ signal_units = 'counts'
+ signal_datatype = 'i8'
+
+ signal = g.createVariable('signal', datatype=signal_datatype, dimensions=('profile', 'range'),
zlib=True)
signal.long_name = 'signal'
- signal.units = channel_parameters['signal_units']
+ signal.units = signal_units
signal.coordinates = "time"
signal.ancillary_variables = "signal_stddev"
# If measured
- signal_stddev = g.createVariable('signal', datatype='i4', dimensions=('profile', 'range'),
+ signal_stddev = g.createVariable('signal_stddev', datatype=signal_datatype, dimensions=('profile', 'range'),
zlib=True)
signal_stddev.long_name = 'signal standard deviation'
- signal_stddev.units = channel_parameters['signal_units']
+ signal_stddev.units = signal_units
signal_stddev.coordinates = "time"
+ signal_stddev.comments = "Only if measured. Should be removed if not."
# Assign variables
- name[:] = channel_name
+ name[:len(channel_name)] = channel_name
laser_rep_rate[:] = channel_parameters['laser_repetition_rate']
emission_energy[:] = channel_parameters['emission_energy']
emission_pol[:] = self._emission_pol_flag(channel_parameters['emission_polarization'])
@@ -256,14 +265,34 @@
detection_mode[:] = self._detection_mode_flag(channel_parameters['detection_mode'])
detection_fwhm[:] = channel_parameters['filter_fwhm']
detection_pol[:] = self._detection_pol_flag(channel_parameters['detection_polarization'])
- polarizer_angle[:] = channel_parameters['polarizer_angle'] # For now, assumed constant.
- dead_time_model[:] = self._deadtime_model_flag(channel_parameters['dead_time_model'])
- dead_time[:] = channel_parameters['dead_time']
+ polarizer_angle[:] = channel_parameters['polarizer_angle'] * np.ones(len(channel.time)) # For now, assumed constant.
+
+ if channel.is_photon_counting:
+ dead_time_model[:] = self._deadtime_model_flag(channel_parameters['dead_time_model'])
+ dead_time[:] = channel_parameters['dead_time']
+
bin_length[:] = channel_parameters['bin_length']
- pulse_delay[:] = channel_parameters['pulse_delay']
+ trigger_delay[:] = channel_parameters['trigger_delay']
+
+ detector_voltage[:] = channel.hv
+
+ if channel.is_analog:
+ adc_range[:] = channel.discriminator
+ adc_bits[:] = channel.adcbits
+ else:
+ discriminator[:] = channel.discriminator
- detector_voltage[:] = channel
+ pulses[:] = channel.laser_shots
+ epoch = datetime.datetime(2017, 1, 1, tzinfo=pytz.utc)
+ seconds_since_epoch = [(t - epoch).total_seconds() for t in channel.time]
+ time[:] = seconds_since_epoch
+ time_bounds[:, 0] = seconds_since_epoch
+ time_bounds[:, 1] = seconds_since_epoch + channel.get_duration()
+
+ bin_time[:] = channel.binwidth * np.arange(len(channel.z))
+
+ signal[:] = channel.matrix
def _deadtime_model_flag(self, model_str):
""" Convert dead-time model string to byte flag.
@@ -364,11 +393,11 @@
: int
Byte encoding of polarization state
"""
- choices = {'linaer': 0,
+ choices = {'linear': 0,
'circular': 1,
'none': 2}
- if pol_string not in choices.keys():
+ if pol_string not in choices.keys():
raise ValueError('Emission polarization not one of {0}: {1}'.format(choices.keys(), pol_string))
return choices[pol_string]
