Skip to content
Snippets Groups Projects
Commit 7417146f authored by David Hoese's avatar David Hoese
Browse files

Merge branch 'marco' into 'master' 

Rewrite of NCDF and Quicklooks generation

See merge request metobs/AossCeilo!1
parents 5380ca0e b775046a
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
MANIFEST.in
+ 1
− 1
View file @ 7417146f
include aossceilo/ceilo.ncml include aossceilo/level_b1/ceilo.ncml
aossceilo/level_b1/ceilo.ncml 0 → 100644
+ 404
− 0
View file @ 7417146f
<?xml version="1.0" encoding="UTF-8"?>
<netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2">
<dimension name="time" isUnlimited="true" />
<dimension name="range" length="256" />
<dimension name="string_len" length="8" />
<attribute name="Conventions" type="String" value="CF-1.3" />
<attribute name="title" type="String" value="Rooftop Instrument Group(RIG) Ceilometer obervations" />
<attribute name="institution" type="String" value="University of Wisconsin
Space Science and Engineering Center
1225 W. Dayton St.
Madison, WI 53706" />
<attribute name="source" type="String" value="surface observation" />
<attribute name="history" type="String" value="" />
<attribute name="references" type="String" value="http://metobs.ssec.wisc.edu" />
<attribute name="comment" type="String" value="" />
<attribute name="site_id" type="String" value="" />
<attribute name="facility_id" type="String" value="" />
<attribute name="sample_int" type="String" value="15 seconds" />
<attribute name="serial_number" type="String" value="" />
<attribute name="Mentor_QC_Field_Information" type="String" value="For each qc field interpret the values as follows:
Basic mentor QC checks (bit values):
==========================================
0x0 = value is within the specified range
0x1 = value is equal to 'missing_value'
0x2 = value is less than the 'valid_min'
0x4 = value is greater than the 'valid_max'
0x8 = value failed the 'valid_delta' check
If the value is a 'missing_value' no min, max, or delta checks are performed.
The delta checks are done by comparing each data value to the one just
prior to it in time. If a previous data value does not exist or is a
'missing_value' the delta check will not be performed.
Note that the delta computation for multi-dimensioned data compares the
absolute value between points in the same spatial location to the previous
point in time.
If the associated non-QC field does not contain any mentor-specified minimum,
maximum, or delta information a qc_field is not generated.
" />
<variable name="base_time" type="int">
<attribute name="long_name" type="String" value="Time of first record in Epoch" />
<attribute name="units" type="String" value="seconds since 1970-1-1 0:00:00 0:00" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="time_offset" shape="time" type="int">
<attribute name="long_name" type="String" value="Time offset from base_time" />
<attribute name="units" type="String" value="" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="time" shape="time" type="int">
<attribute name="long_name" type="String" value="Time offset from midnight" />
<attribute name="units" type="String" value="" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="range" shape="range" type="double">
<attribute name="long_name" type="String" value="Distance to the center of the corresponding range bin." />
<attribute name="units" type="String" value="m" />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="accuracy" type="String" value="30 meters" />
</variable>
<variable name="first_cbh" shape="time" type="double">
<attribute name="long_name" type="String" value="Lowest cloud base height detected." />
<attribute name="standard_name" type="String" value="cloud_base_altitude" />
<attribute name="units" type="String" value="m" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="valid_min" type="double" value="0" />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="description" type="String" value=" If detection_status = 1, 2 or 3,
lowest cloud base height measured.
If detection_status = 4 this value
will be stored as -9999 and a value
for vertical_visibility will be stored.
If detection_status = 0 or 5, -9999 will
be stored for this value and vertical_visibility.
(5 slashes were read from field)." />
</variable>
<variable name="vertical_visibility" shape="time" type="double">
<attribute name="long_name" type="String" value="Vertical visibility" />
<attribute name="standard_name" type="String" value="visibility_in_air" />
<attribute name="units" type="String" value="m" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="description" type="String" value="This field will only be stored if detection_status
is = 4. In all other cases, a value of -9999 will
be recorded." />
</variable>
<variable name="second_cbh" shape="time" type="double">
<attribute name="long_name" type="String" value="Second lowest cloud base height" />
<attribute name="standard_name" type="String" value="cloud_base_altitude" />
<attribute name="units" type="String" value="m" />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="valid_min" type="double" value="0" />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="description" type="String" value=" If detection_status = 2 or 3, a second
cloud base was detected. This value is
the distance from the ground to the second
cloud base.
If detection_status = 4, -9999 will be recorded
here and a value for alt_highest_signal
will be recorded.
If detection_status = 5, -9999 will be recorded
for this value and for alt_highest_signal" />
</variable>
<variable name="alt_highest_signal" shape="time" type="double">
<attribute name="long_name" type="String" value="Altitude of highest signal" />
<attribute name="units" type="String" value="m" />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="valid_min" type="double" value="0" />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="description" type="String" value=" If detection_status = 4, this value is the
distance from the ground to the highest signal
detected.
If detection_status = 0, 1, 2, 3 or 5, -9999
will be recorded." />
</variable>
<variable name="third_cbh" shape="time" type="double">
<attribute name="long_name" type="String" value="Third cloud base height" />
<attribute name="standard_name" type="String" value="cloud_base_altitude" />
<attribute name="units" type="String" value="m" />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="valid_min" type="double" value="0" />
<attribute name="valid_max" type="double" value="7620." />
<attribute name="description" type="String" value=" If detection_status = 3, this is the distance
from the ground to the base of the highest
cloud.
If detection = 0, 1, 2, 4 or 5, -9999 will be
stored here. In this case, ///// was read from field." />
</variable>
<variable name="sum_backscatter" shape="time" type="double">
<attribute name="long_name" type="String" value="SUM of detected and normalized backscatter" />
<attribute name="units" type="String" value="1.0 sr-1" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="999." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="description" type="String" value="Multiplied by scaling factor
times 10^4. At scaling factor 100 the sum range,
which is 0 to 999, this corresponds to integrated
backscatter, which ranges from 0 to 0.0999 srad^-1" />
</variable>
<variable name="backscatter" shape="time range" type="double">
<attribute name="long_name" type="String" value="Backscatter" />
<attribute name="units" type="String" value="1.0E-7 m-1.sr-1" />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="description" type="String" value="Data is range and sensitivity normalized backscatter, units, (10000 - srad - km)^-1 unless otherwise scaled by
scaling parameter.
The length of each range bin corresponds to
(c/2)*(200ns +/14 ns) = (approximately 30 m)." />
</variable>
<variable name="detection_status" shape="time" type="int">
<attribute name="long_name" type="String" value="Detection status." />
<attribute name="_FillValue" type="int" value="-9999" />
<attribute name="valid_min" type="int" value="0" />
<attribute name="valid_max" type="int" value="5" />
<attribute name="description" type="String" value="0=No significant backscatter
1=One cloud base detected
2=Two cloud bases detected
3=Three cloud bases detected
4=Full obscuration determined but no cloud
base detected.
5=Some obscuration detected but determined
to be transparent." />
</variable>
<variable name="status_flag" shape="time" type="char">
<attribute name="long_name" type="String" value="Ceilometer status indicator." />
<attribute name="description" type="String" value="Warning and alarm information as follows
0 = Self-check OK
W = At least one Warning active, no Alarms
A = At least one Alarm active" />
</variable>
<variable name="laser_pulse_energy" shape="time" type="double">
<attribute name="long_name" type="String" value="Laser pulse energy" />
<attribute name="units" type="String" value="%" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="999." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="comment" type="String" value="Percentage of nominal factory setting." />
</variable>
<variable name="laser_temperature" shape="time" type="double">
<attribute name="long_name" type="String" value="Laser temperature" />
<attribute name="units" type="String" value="degC" />
<attribute name="valid_min" type="double" value="-50." />
<attribute name="valid_max" type="double" value="99." />
<attribute name="_FillValue" type="double" value="-9999." />
</variable>
<variable name="receiver_sensitivity" shape="time" type="double">
<attribute name="long_name" type="String" value="Receiver sensitivity" />
<attribute name="units" type="String" value="%" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="999." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="comment" type="String" value="Percentage of nominal factory setting." />
</variable>
<variable name="window_contamination" shape="time" type="double">
<attribute name="long_name" type="String" value="Window contamination" />
<attribute name="units" type="String" value="0.0010 V" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="2500." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="comment" type="String" value="Millivolts measured at the internal ADC input" />
</variable>
<variable name="tilt_angle" shape="time" type="double">
<attribute name="long_name" type="String" value="Tilt angle" />
<attribute name="units" type="String" value="degrees" />
<attribute name="valid_min" type="double" value="-15." />
<attribute name="valid_max" type="double" value="90." />
<attribute name="valid_delta" type="double" value="1." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="description" type="String" value="If Angle Correction is Off detection values are distances not altitudes" />
</variable>
<variable name="background_light" shape="time" type="double">
<attribute name="long_name" type="String" value="Background light" />
<attribute name="units" type="String" value="0.0010 V" />
<attribute name="valid_min" type="double" value="0." />
<attribute name="valid_max" type="double" value="2500." />
<attribute name="_FillValue" type="double" value="-9999." />
<attribute name="comment" type="String" value="Millivolts measured at the internal ADC input" />
</variable>
<variable name="measurement_parameters" shape="time string_len" type="char">
<attribute name="long_name" type="String" value="6 character string describing instrument measurement parameters" />
<attribute name="comment" type="String" value="(L)ong/(S)hort - (F)req - pulse qty r ^(7) + 1
- gain (H)igh/(L)ow - bandwidth (N)arrow/(W)ide -
sampling (5) / (1)0 / (2)0 Mhz" />
</variable>
<variable name="status_string" shape="time string_len" type="char">
<attribute name="long_name" type="String" value="Warning and alarm status bits" />
<attribute name="comment1" type="String" value="Has Form FEDCBA98. Each character is a hexadecimal representation of
four bits, altogether 32 bits. See additional comments for interpretation." />
<attribute name="commentF" type="String" value="(8000 0000) = Laser temperature shut-off(A)
(4000 0000)= Laser failure(A)
(2000 0000)= Receiver failure(A)
(1000 0000)= Voltage failure(A)" />
<attribute name="commentE" type="String" value="(0800 0000) = Profile error(A)" />
<attribute name="commentD" type="String" value="(0080 0000) = Window contaminated(W)
(0040 0000)= Battery low(W)
(0020 0000)= Laser power low(W)
(0010 0000)= Laser temperature high or low(W)" />
<attribute name="commentC" type="String" value="(0008 0000) = Internal temperature high or low(W)
(0004 0000)= Voltage high or low(W)
(0002 0000)= Relative humidity is > 85 % (W)
(0001 0000)= Receiver optical cross-talk compensation poor(W)" />
<attribute name="commentB" type="String" value="(0000 8000) = Fan suspect(W)
(0000 4000)= Profile warning(W)" />
<attribute name="commentA" type="String" value="(0000 0800) = Blower is ON
(0000 0400)= Blower heater is ON
(0000 0200)= Internal heater is ON
(0000 0100)= Units are METERS if ON, else FEET." />
<attribute name="comment9" type="String" value="(0000 0080) = Polling mode is ON
(0000 0040)= Working from battery
(0000 0020)= Single sequence mode is ON
(0000 0010)= Manual Settings are effective." />
<attribute name="comment8" type="String" value="(0000 0008) = Tilt angle is > 45 degrees." />
</variable>
<variable name="lat" type="double">
<attribute name="standard_name" type="String" value="latitude" />
<attribute name="long_name" type="String" value="north latitude" />
<attribute name="units" type="String" value="degrees_north" />
<attribute name="valid_min" type="double" value="-90." />
<attribute name="valid_max" type="double" value="90." />
<attribute name="_FillValue" type="double" value="-9999." />
</variable>
<variable name="lon" type="double">
<attribute name="standard_name" type="String" value="longitude" />
<attribute name="long_name" type="String" value="east longitude" />
<attribute name="units" type="String" value="degrees_east" />
<attribute name="valid_min" type="double" value="-180." />
<attribute name="valid_max" type="double" value="180." />
<attribute name="_FillValue" type="double" value="-9999." />
</variable>
<variable name="alt" type="double">
<attribute name="standard_name" type="String" value="altitude" />
<attribute name="long_name" type="String" value="altitude" />
<attribute name="positive" type="String" value="up" />
<attribute name="units" type="String" value="m" />
<attribute name="_FillValue" type="double" value="-9999." />
</variable>
<variable name="qc_time" shape="time" type="int">
<attribute name="long_name" type="String" value="Results of quality checks on sample time" />
<attribute name="description" type="String" value="The qc_time values are calculated by comparing each sample
time with the previous time (i.e. delta_t = t[n] - t[n-1]).
If the 'qc_check_prior' flag is set the first sample time
from a new raw file will be compared against the time just
previous to it in the stored data. If the 'qc_check_prior'
flag is not set the qc_time value for the first sample time
will be set to 0
The qc_time bit values are as follows:
=========================================================
0x0 = delta time is within the specified range
0x1 = delta time is equal to 0, duplicate sample times
0x2 = delta time is less than the 'delta_t_lower_limit'
0x4 = delta time is greater than the 'delta_t_upper_limit'
" />
<attribute name="delta_t_lower_limit" type="int" value="12" />
<attribute name="delta_t_upper_limit" type="int" value="18" />
<attribute name="prior_sample_flag" type="int" value="1" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_first_cbh" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Lowest cloud base height detected." />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_vertical_visibility" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Vertical visibility" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_second_cbh" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Second lowest cloud base height" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_alt_highest_signal" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Altitude of highest signal" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_third_cbh" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Third cloud base height" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_laser_pulse_energy" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Laser pulse energy" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_laser_temperature" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Laser temperature" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
<variable name="qc_tilt_angle" shape="time" type="int">
<attribute name="long_name" type="String" value="Quality check results on field: Tilt angle" />
<attribute name="_FillValue" type="int" value="-9999" />
</variable>
</netcdf>
aossceilo/level_b1/message.py 0 → 100644
+ 583
− 0
View file @ 7417146f
"""
For processing Ceilometer CT25K Messages
"""
import logging
from calendar import timegm
from datetime import datetime
from numpy import array
LOG = logging.getLogger(__name__)
#: Value to use for NetCDF missing_value attribute
missing_value = -9999
#: Function to determine if a value is a missing_value or not
def is_missing(val):
return val == missing_value
#: number of allowable missing values before backscatter is considered invalid
BK_THRESH = 256 * 0.25
START_OF_HEADING = "\x01"
START_OF_TEXT = "\x02"
END_OF_TEXT = "\x03"
MANUAL_BLOWER_CONTROL = 2
HIGH_BACKGROUND_RADIANCE = 4
TILT_ANGLE_45_PLUS = 8
MANUAL_SETTINGS_EFF = 16
SINGLE_SEQ_ON = 32
WORKING_FROM_BATT = 64
POLLING_ON = 128
UNITS_IN_METERS = 256
INTERNAL_HEATER_ON = 512
BLOWER_HEATER_ON = 1024
BLOWER_ON = 2048
BLOWER_SUSPECT = 32768
RCV_OPTICAL_XTLK_POOR = 65536
REL_HUMIDITY_HIGH = 131072
VOLTAGE_ABNORMAL = 262144
INT_TEMP_HIGH = 524288
LASER_TEMP_HIGH = 1048576
LASER_PWR_LOW = 2097152
BATTERY_LOW = 4194304
WINDOW_CONTAMINATED = 8388608
VOLTAGE_FAILURE = 268435456
RECEIVER_FAILURE = 536870912
LASER_FAILURE = 1073741824
LASER_TEMP_SHUTOFF = 2147483648
WARNING_FLAGS = [
BLOWER_SUSPECT,
RCV_OPTICAL_XTLK_POOR,
REL_HUMIDITY_HIGH,
VOLTAGE_ABNORMAL,
INT_TEMP_HIGH,
LASER_TEMP_HIGH,
LASER_PWR_LOW,
BATTERY_LOW,
WINDOW_CONTAMINATED,
]
WARNING_MSGS = [
"Blower suspect",
"Receiver optical cross-talk poor",
"Relative humidity is greater than 85%",
"Voltage is low or high",
"Internal temperature is high or low",
"Laser temperature is high",
"Laser power is low",
"Battery is low",
"Window is contaminated",
]
INFO_FLAGS = [
MANUAL_BLOWER_CONTROL,
HIGH_BACKGROUND_RADIANCE,
TILT_ANGLE_45_PLUS,
MANUAL_SETTINGS_EFF,
SINGLE_SEQ_ON,
WORKING_FROM_BATT,
POLLING_ON,
UNITS_IN_METERS,
INTERNAL_HEATER_ON,
BLOWER_HEATER_ON,
BLOWER_ON,
]
INFO_MSGS = [
"Using manual blower control",
"Background radiance is high",
"Tilt angle is greater than 45 degrees",
"Manual settings are effective",
"Single sequence mode is ON",
"Working from the battery",
"Polling mode is ON",
"Units are in Meters",
"Internal heater is ON",
"Blower heater is ON",
"Blower is ON",
]
ALARM_FLAGS = [
VOLTAGE_FAILURE,
RECEIVER_FAILURE,
LASER_FAILURE,
LASER_TEMP_SHUTOFF,
]
ALARM_MSGS = [
"VOLTAGE FAILURE",
"RECEIVER FAILURE",
"LASER FAILURE",
"LASER TEMPERATURE SHUT_OFF",
]
class MessageError(RuntimeError):
"""General message error."""
class Message2(object):
NUM_LINES = 20
def __init__(
self, lines, stamp,
):
"""Initialize this message, disecting the various message components
according to the CT25K Message 2 format.
The first line in a message is the header and it must start with a
SOH and end with SOT, and the entire message must end with EOT
ASCII chars to be valid. Assertion errors are raised if these
conditions are not met.
The backscatter is initially filled with L{missing_value} so if a
line is the incorrect length all of it's values will be missing_value.
Likewise, if a value cannot be parsed as a HEX string it's value will
be missing_value.
@param lines: Exactly the 20 lines that comprise a valid message
@type lines: list
@param stamp: The message time. If the 'stamp' is a naieve datetime the
tzinfo will be set to metobs.util.time.UTC, otherwise time operations
will proceed with provided tzinfo.
@type stamp: datetime
@raises MessageError: If this instance cannot be created due to an error
parsing.
"""
assert len(lines) == self.NUM_LINES, (
"A Message2 must contain %s lines" % self.NUM_LINES
)
self._epoch = timegm(stamp.utctimetuple())
self._lines = lines
self._stamp = stamp
self._header = lines[0]
# strip non-printables
self._header = self._header.replace(START_OF_HEADING, "",)
self._header = self._header.replace(START_OF_TEXT, "",)
self._msg_num = int(self.header[5])
if self._msg_num != 2:
raise MessageError(
"Invalid message number", self.header,
)
self._status_string = lines[1][21:29]
self._status_flag = lines[1][1]
self._detection_status = _int(lines[1][0])
self._first_cbh = missing_value
self._second_cbh = missing_value
self._third_cbh = missing_value
if self.detection_status == 3:
self._third_cbh = _int(lines[1][15:20])
if self.detection_status == 2:
self._second_cbh = _int(lines[1][9:14])
if self.detection_status == 1:
self._first_cbh = _int(lines[1][3:8])
if self.detection_status == 4:
self._vertical_visibility = _int(lines[1][3:8])
self._alt_highest_signal = _int(lines[1][9:14])
else:
self._vertical_visibility = missing_value
self._alt_highest_signal = missing_value
meas_params = lines[2].split()
if len(meas_params) < 10:
LOG.warn(
"Invalid measurement parameters for message with time %s",
self.epoch,
)
self._scale = _int(meas_params[0])
self._measurement_mode = _str(meas_params[1])
self._laser_pulse_energy = _float(meas_params[2])
self._laser_temperature = _float(meas_params[3])
self._receiver_sensitivity = _float(meas_params[4])
self._window_contamination = _float(meas_params[5])
self._tilt_angle = _float(meas_params[6])
self._background_light = _float(meas_params[7])
self._measurement_parameters = _str(meas_params[8])
self._sum_backscatter = _float(meas_params[9])
self._backscatter = parse_backscatter(lines[3:])
# backscatter can contain no more than BK_THRESH missing values
missing = [b for b in self._backscatter if is_missing(b)]
if len(missing) >= BK_THRESH:
raise MessageError("Backscatter errors exceeds threshold")
@property
def epoch(self):
return self._epoch
@property
def lines(self):
return self._lines
@property
def stamp(self):
return self._stamp
@property
def header(self,):
return self._header
@property
def msg_num(self,):
return self._msg_num
@property
def status_string(self,):
return self._status_string
@property
def status_flag(self,):
return self._status_flag
@property
def detection_status(self,):
return self._detection_status
@property
def cbh(self):
return (
self._first_cbh,
self._second_cbh,
self._third_cbh,
)
@property
def first_cbh(self,):
return self._first_cbh
@property
def second_cbh(self,):
return self._second_cbh
@property
def third_cbh(self,):
return self._third_cbh
@property
def vertical_visibility(self,):
return self._vertical_visibility
@property
def alt_highest_signal(self,):
return self._alt_highest_signal
@property
def scale(self):
return self._scale
@property
def measurement_mode(self,):
return self._measurement_mode
@property
def laser_pulse_energy(self,):
return self._laser_pulse_energy
@property
def laser_temperature(self,):
return self._laser_temperature
@property
def receiver_sensitivity(self,):
return self._receiver_sensitivity
@property
def window_contamination(self,):
return self._window_contamination
@property
def tilt_angle(self,):
return self._tilt_angle
@property
def background_light(self,):
return self._background_light
@property
def measurement_parameters(self,):
return self._measurement_parameters
@property
def sum_backscatter(self,):
return self._sum_backscatter
@property
def backscatter(self,):
return self._backscatter
def warning(self,):
"""
@return: true if there is a warning indicated in the detection status string.
"""
return str(self.detection_status) in "wW"
def alarm(self):
"""
@return: true if there is a alarm indicated in the detection status string.
"""
return str(self.detection_status) in "aA"
def check_status(self, flag):
"""
Check the status of a particular flag.
@param flag: One of the flag constants
@type flag: number
@return: 0 if the flag is not set, otherwise the value of the flag.
"""
return int(self.status_string, 16,) & flag
def get_status_messages(self,):
"""
Get all status messages indicated in the device status string.
@return: 'alarms', 'warnings', and 'informational' messages.
"""
messages = {
"alarms": [],
"warnings": [],
"informational": [],
}
for idx in range(len(WARNING_FLAGS)):
if int(self.status_string, 16,) & WARNING_FLAGS[idx]:
messages["warnings"].append(WARNING_MSGS[idx])
for idx in range(len(ALARM_FLAGS)):
if int(self.status_string, 16,) & ALARM_FLAGS[idx]:
messages["alarms"].append(ALARM_MSGS[idx])
for idx in range(len(INFO_FLAGS)):
if int(self.status_string, 16,) & INFO_FLAGS[idx]:
messages["informational"].append(INFO_MSGS[idx])
return messages
def __eq__(self, that):
"""
Comparisons are done on message timestamps.
"""
return self.stamp == that.stamp
def __gt__(self, that):
"""
Comparisons are done on message timestamps.
"""
return self.stamp > that.stamp
def __lt__(self, that):
"""
Comparisons are done on message timestamps.
"""
return self.stamp < that.stamp
def __str__(self,):
return "\n".join([str(self.epoch)] + self.lines)
class Message7(Message2):
"""Same as message 2 with additional sky condition line.
"""
NUM_LINES = 21
def __init__(
self, lines, stamp,
):
super(Message7, self,).__init__(
lines, stamp,
)
self._sky_cond = parse_sky_condition(self.lines[-1])
@property
def sky_condition(self,):
return self._sky_cond
def hex_to_int(hval,):
# convert from 2's complement
val = int(hval, 16)
if val > 32767:
val = -65536 + val
return val
def _get(fcn, v, d):
try:
return fcn(v)
except:
pass
return d
def _float(v):
return _get(float, v, missing_value,)
def _int(v):
return _get(int, v, missing_value,)
def _str(v):
return _get(str, v, "")
def parse_sky_condition(line, feet=False):
parts = line.strip().split()
sky_cond = []
for idx in range(0, 8, 2):
amt = int(parts[idx])
if amt > 0 and amt < 9:
hght = int(parts[idx + 1])
if feet: # feet
hght *= 100
else: # meters
hght *= 10
sky_cond.append((amt, hght,))
elif amt == 0: # no clouds
sky_cond.append((0, missing_value,))
elif amt == 9: # vertical visibility
# FIXME: ignoring this because we already have it
# from line 2, right?
pass
else: # amt in (-1, 99): # no data
sky_cond.append((missing_value, missing_value,))
return sky_cond
def parse_backscatter(data,):
if isinstance(data, str):
lines = data.split("\n")
else:
lines = data
if len(lines) < 16:
raise MessageError("Expected 16 lines for backscatter")
backscatter = [missing_value] * 256
bidx = 0
for line in lines[:16]:
if len(line) < 67: # 3 chars for height, 64 for data
bidx += 16
continue
for idx in range(3, 67, 4):
hval = line[idx : idx + 4]
try:
backscatter[bidx] = hex_to_int(hval)
except:
pass
bidx += 1
return backscatter
def get_message_times(messages,):
times = array([missing_value] * len(messages), "l",)
for i in range(len(messages)):
times[i] = timegm(messages[i].stamp.utctimetuple())
return times
def get_message_met_data(messages,):
"""Compile message data into a single data dictionary.
@param messages: messages
@type messages: list
@return: Ceilometer data where keys are variable names from the netcdf file.
"""
data = dict(
vertical_visibility=array([missing_value] * len(messages)),
alt_highest_signal=array([missing_value] * len(messages)),
first_cbh=array([missing_value] * len(messages)),
second_cbh=array([missing_value] * len(messages)),
third_cbh=array([missing_value] * len(messages)),
backscatter=array(
[[missing_value] * 256 for idx in range(len(messages))]
),
)
for i in range(len(messages)):
data["alt_highest_signal"][i] = messages[i].highest_signal
data["vertical_visibility"][i] = messages[i].vertical_vis
data["first_cbh"][i] = messages[i].bases[0]
data["second_cbh"][i] = messages[i].bases[1]
data["third_cbh"][i] = messages[i].bases[2]
data["backscatter"][i][:] = array(messages[i].backscatter)
return data
def get_message_hk_data(messages,):
"""
"""
return {}
def load_messages(
url, on_error=None,
):
"""Load messages from a message data file. All lines in the data file
are skipped until a proper start of message is encounterd.
@param url: URL to a message data file as generated by the Ceilometer
ingest software
@type url: str
@returns: Sequence of Messages sorted from oldest to newest.
"""
fp = open(url)
line = fp.readline()
# readline returned '', so there are no lines
if not line:
return tuple()
# message must start with epoch time, so disregard partial
# messages until we reach one
while line and not line.strip().isdigit():
line = fp.readline()
# start of first complete message
messages = []
done = not line # if we have a line, we're not done
while not done:
lines = [line] + [fp.readline() for i in range(Message2.NUM_LINES)]
if "" in lines:
done = True
continue
stamp = datetime.fromtimestamp(int(lines[0]))
try:
messages.append(
Message2([line.strip() for line in lines[1:]], stamp,)
)
except Exception as e:
if hasattr(on_error, "__call__",):
on_error(
lines, e,
)
# disregard possible garbage between messages
line = fp.readline()
while line and not line.strip().isdigit():
line = fp.readline()
messages.sort()
return tuple(messages)
aossceilo/level_b1/nc.py 0 → 100644
+ 256
− 0
View file @ 7417146f
import numpy
from calendar import timegm
from metobscommon.util import mytime
from netCDF4 import Dataset
import pkg_resources
import sys
import os
import argparse
from aossceilo.level_b1 import message
try:
from xml.etree.ElementTree import XML
except ImportError:
from elementtree.ElementTree import XML
get_value = (
lambda var, value: var._FillValue
if numpy.isnan(value) and hasattr(var, "_FillValue",)
else value
)
def create_nc(input, out=None):
input = os.path.realpath(input)
if out is None:
out = f"{os.path.splitext(input)[0]}.nc"
# SETUP
ncml = XML(pkg_resources.resource_string(__name__, "ceilo.ncml",))
nc = Dataset(out, "w")
def tag(s):
return f"{{http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2}}{s}"
# Dimensions
for e in ncml.findall(tag("dimension")):
name = e.attrib["name"]
if "isUnlimited" in e.attrib:
nc.createDimension(
name, None,
) # None for unlimited length
else:
nc.createDimension(
name, int(e.attrib["length"]),
)
# Attributes
for e in ncml.findall(tag("attribute")):
name = e.attrib["name"]
value = e.attrib["value"]
setattr(
nc, name, value,
)
# Variables
var_map = {
"int": ("i4", numpy.int32,),
"char": ("S1", str,),
"float": ("f4", numpy.float32,),
"short": ("i2", numpy.int16,),
"double": ("f8", numpy.float64,),
"String": ("S1", str,),
}
attr_map = {
"int": ("i4", int,),
"char": ("S1", str,),
"float": ("f4", float,),
"short": ("i2", int,),
"double": ("f8", float,),
"String": ("S1", str,),
}
for v in ncml.findall(tag("variable")):
v_name = v.attrib["name"]
v_type = var_map[v.attrib["type"]][0]
shape = v.attrib.get("shape")
attributes = v.findall(tag("attribute"))
try:
fill_value_attr = v.find(tag('attribute[@name="_FillValue"]'))
fill_value = fill_value_attr.attrib["value"]
attributes.remove(fill_value_attr)
except AttributeError:
fill_value = None
if shape:
var = nc.createVariable(
v_name,
v_type,
dimensions=tuple(shape.split(" ")),
fill_value=fill_value,
)
else:
var = nc.createVariable(v_name, v_type, fill_value=fill_value,)
for a in attributes:
a_name = a.attrib["name"]
a_type = attr_map[a.attrib["type"]][1]
value = a.attrib["value"]
setattr(
var, a_name, a_type(value),
)
nc.sync()
ssec_loc = (
43.08543,
-89.271632,
)
nc.variables["lat"][:] = ssec_loc[0]
nc.variables["lon"][:] = ssec_loc[1]
nc.sync()
messages = message.load_messages(input)
base = mytime.datetime_to_epoch(messages[0].stamp)
var = nc.variables["base_time"]
setattr(
var, "units", "seconds since 1970-01-01 00:00:00 0:00",
)
var[:] = base
times = numpy.array(
[numpy.int64(timegm(m.stamp.utctimetuple())) for m in messages]
)
met_data = get_message_met_data(nc, messages)
hk_data = get_message_hk_data(nc, messages)
midnight = mytime.datetime_to_epoch(
mytime.day_begin(messages[len(messages) // 2].stamp)
)
offsets = times - midnight # use scalar array subtraction
if offsets.dtype != numpy.int32:
offsets = numpy.array(offsets, numpy.int32,)
var = nc.variables["time"]
var[: len(offsets)] = offsets
var.units = messages[len(messages) // 2].stamp.strftime(
"seconds since %Y-%m-%d 00:00:00 0:00"
)
# time offsets from base_time
base_secs = mytime.datetime_to_epoch(mytime.day_begin(messages[0].stamp))
offsets = times - base_secs
var = nc.variables["time_offset"]
var.units = messages[0].stamp.strftime(
"seconds since %Y-%m-%d %H:%M:%S 0:00"
)
if offsets.dtype != numpy.int32:
offsets = numpy.array(offsets, numpy.int32,)
var[: len(offsets)] = offsets
for val_dict in [
met_data,
hk_data,
]:
for var_name in val_dict.keys():
arr = val_dict[var_name]
var = nc.variables[var_name]
if len(arr.shape) == 1:
var[: arr.shape[0]] = arr
else:
var[: arr.shape[0], : arr.shape[1],] = arr
nc.sync()
return nc
def get_message_hk_data(nc, messages):
var_names = (
"laser_pulse_energy",
"laser_temperature",
"receiver_sensitivity",
"window_contamination",
"tilt_angle",
"background_light",
"sum_backscatter",
"tilt_angle",
)
str_names = [
"status_string",
"measurement_parameters",
]
data = {}
str_len = getattr(nc.dimensions["string_len"], "size",)
data["detection_status"] = numpy.array(
[m.detection_status for m in messages], numpy.int32,
)
data["range"] = numpy.array(
range(15, 256 * 30, 30,)
) # center of range buckets. Buckets are 30m
# so start at 15 and count 256 buckets.
for var_name in str_names:
data[var_name] = numpy.zeros((len(messages), str_len,))
for var_name in var_names:
data[var_name] = numpy.zeros(len(messages))
for i in range(len(messages)):
for var_name in var_names:
data[var_name][i] = get_value(
nc.variables[var_name], getattr(messages[i], var_name,),
)
for var_name in str_names:
s = getattr(messages[i], var_name,)
data[var_name][i][:] = list(map(ord, list(s.ljust(str_len)),))
return data
def get_message_met_data(nc, messages):
data = dict(
vertical_visibility=numpy.zeros(len(messages)),
alt_highest_signal=numpy.zeros(len(messages)),
first_cbh=numpy.zeros(len(messages)),
second_cbh=numpy.zeros(len(messages)),
third_cbh=numpy.zeros(len(messages)),
backscatter=numpy.array(
[numpy.zeros(256) for idx in range(len(messages))]
),
)
for i in range(len(messages)):
data["alt_highest_signal"][i] = get_value(
nc.variables["alt_highest_signal"], messages[i].alt_highest_signal,
)
data["vertical_visibility"][i] = get_value(
nc.variables["vertical_visibility"],
messages[i].vertical_visibility,
)
data["first_cbh"][i] = get_value(
nc.variables["first_cbh"], messages[i].first_cbh,
)
data["second_cbh"][i] = get_value(
nc.variables["second_cbh"], messages[i].second_cbh,
)
data["third_cbh"][i] = get_value(
nc.variables["third_cbh"], messages[i].third_cbh,
)
data["backscatter"][i][:] = numpy.array(messages[i].backscatter)
return data
def main():
parser = argparse.ArgumentParser(
description="Library for creating and manipulating Viasala CT25K Ceilometer NetCDF files."
)
parser.add_argument(
"input", help="ASCII filepath",
)
parser.add_argument(
"-o", "--output-file", help="output file name/location",
)
args = parser.parse_args()
create_nc(
args.input, out=args.output_file,
)
if __name__ == "__main__":
sys.exit(main())
aossceilo/tests/rig_ceilo.2021-10-01.ascii 0 → 100644
+ 120960
− 0
View file @ 7417146f
This diff is collapsed.
aossceilo/tests/test_nc.py 0 → 100644
+ 81
− 0
View file @ 7417146f
from aossceilo.level_b1 import message, nc
import os
import numpy
ascii_test_file = os.path.join(
os.path.dirname(__file__), "rig_ceilo.2021-10-01.ascii"
)
def test_load_messages():
# Normal Operation
messages = message.load_messages(ascii_test_file)
for m in messages:
assert type(m) == message.Message2
assert hasattr(m, "stamp")
def test_empty_load_messages(tmp_path):
# test with empty file
filename = tmp_path / "temp.ascii"
with open(filename, "w"):
pass
messages = message.load_messages(filename)
os.remove(filename)
assert messages == tuple()
def test_create_nc(tmp_path):
# create nc file then check attributes
filename = tmp_path / "temp.nc"
nc_file = nc.create_nc(ascii_test_file, out=filename)
# check global attributes
attrs = (
"Conventions",
"title",
"institution",
"source",
"history",
"references",
"comment",
"site_id",
"facility_id",
"sample_int",
"serial_number",
"Mentor_QC_Field_Information",
)
for attr in attrs:
assert hasattr(nc_file, attr)
# check dimensions
for dim in nc_file.dimensions.values():
assert hasattr(dim, "name")
assert hasattr(dim, "size")
# check variables
expect_fill = {"second_cbh", "third_cbh", "alt"}
for var in nc_file.variables.values():
assert hasattr(var, "long_name")
if var.name == "backscatter":
assert hasattr(var, "description")
assert hasattr(var, "_FillValue")
elif var.name == "sum_backscatter":
assert hasattr(var, "description")
assert hasattr(var, "_FillValue")
assert hasattr(var, "units")
assert hasattr(var, "valid_min")
assert hasattr(var, "valid_max")
elif not var.shape:
assert hasattr(var, "units")
assert hasattr(var, "_FillValue")
if (
var[:].dtype != "S1"
and var.name not in expect_fill
and var.name[:2] != "qc"
):
assert var[:].all() is not numpy.ma.masked
os.remove(filename)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment