diff --git a/.gitignore b/.gitignore index 405ba97dad..bc28853000 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ *.org .project .pydevproject +.spyproject *.rej .settings/ .spyproject diff --git a/pyart/default_config.py b/pyart/default_config.py index 4a11addd5c..931f071a3a 100644 --- a/pyart/default_config.py +++ b/pyart/default_config.py @@ -1016,6 +1016,7 @@ 173: radar_estimated_rain_rate, # Digital User-Selectable Accum. 174: radar_estimated_rain_rate, # Digital 1 hr Diff. Accum. 175: radar_estimated_rain_rate, # Digital Storm Total Diff. Accum. + 176: radar_estimated_rain_rate, # Digital Inst. Precipitation Rate 177: radar_echo_classification, # Hybrid Hydrometeor Classification 181: reflectivity, # Base Reflectivity 182: velocity, # Base Velocity diff --git a/pyart/io/nexrad_level3.py b/pyart/io/nexrad_level3.py index 904c1a6471..c4501155ef 100644 --- a/pyart/io/nexrad_level3.py +++ b/pyart/io/nexrad_level3.py @@ -66,8 +66,10 @@ """ import bz2 +from collections import namedtuple from datetime import datetime, timedelta import struct +from xdrlib import Unpacker import warnings import numpy as np @@ -139,26 +141,30 @@ def __init__(self, filename): supports max version of %d. Most recent product version has not \ yet been implemented/tested.' % (ver,supp_ver), UserWarning) - # uncompressed symbology block if necessary + # Uncompress symbology block if necessary if buf[bpos:bpos+2] == b'BZ': buf2 = bz2.decompress(buf[bpos:]) else: buf2 = buf[bpos:] - self._read_symbology_block(buf2) - - def close(self): - """ Close the file. """ - self._fh.close() - - def _read_symbology_block(self, buf2): - """ Read symbology block. """ # Read and decode symbology header self.symbology_header = _unpack_from_buf(buf2, 0, SYMBOLOGY_HEADER) - # Read radial packets packet_code = struct.unpack('>h', buf2[16:18])[0] assert packet_code in SUPPORTED_PACKET_CODES + + bpos = 16 + + if packet_code == 28: + self._read_symbology_block_28(buf2, bpos, packet_code) + else: + self._read_symbology_block(buf2, bpos, packet_code) + + def close(self): + """ Close the file. """ + self._fh.close() + + def _read_symbology_block(self, buf2, pos, packet_code): self.packet_header = _unpack_from_buf(buf2, 16, RADIAL_PACKET_HEADER) self.radial_headers = [] nbins = self.packet_header['nbins'] @@ -186,6 +192,34 @@ def _read_symbology_block(self, buf2): pos += rle_size self.radial_headers.append(radial_header) + def _read_symbology_block_28(self, buf2, bpos, packet_code): + """ Read symbology block for Packet Code 28 (Product 176). """ + self.packet_header = _unpack_from_buf(buf2, bpos, GEN_DATA_PACK_HEADER) + bpos += 8 + + # Read number of bytes (2 HW) and return + num_bytes = self.packet_header['num_bytes'] + hunk = buf2[bpos : bpos+num_bytes] + xdrparser = Level3XDRParser(hunk) + self.gen_data_pack = xdrparser(packet_code) + + # Rearrange some of the info so it matches the format of packet codes + # 16 and AF1F so method calls can be done properly + self.packet_header['nradials'] = len(self.gen_data_pack['components'].radials) + nradials = self.packet_header['nradials'] + self.packet_header['nbins'] = self.gen_data_pack['components'].radials[0].num_bins + nbins = self.packet_header['nbins'] + self.packet_header['first_bin'] = self.gen_data_pack['components'].first_gate + self.packet_header['range_scale'] = 1000 # 1000m in 1 km + + # Read azimuths + self.azimuths = [rad.azimuth for rad in self.gen_data_pack['components'].radials] + + # Pull each radial's data into an array + self.raw_data = np.empty((nradials, nbins), dtype='uint8') + for i in range(0,nradials): + self.raw_data[i,:] = self.gen_data_pack['components'].radials[i].data + def get_location(self): """ Return the latitude, longitude and height of the radar. """ latitude = self.prod_descr['latitude'] * 0.001 @@ -195,8 +229,11 @@ def get_location(self): def get_azimuth(self): """ Return an array of starting azimuth angles in degrees. """ - azimuths = [d['angle_start'] for d in self.radial_headers] - return np.array(azimuths, dtype='float32') * 0.1 + if self.packet_header['packet_code'] == 28: + azimuths = self.azimuths + else: + azimuths = [d['angle_start'] * 0.1 for d in self.radial_headers] + return np.array(azimuths, dtype='float32') def get_range(self): """ Return an array of gate range spacing in meters. """ @@ -221,7 +258,7 @@ def get_volume_start_datetime(self): self.prod_descr['vol_scan_time']) def get_data(self): - """ Return an masked array containing the field data. """ + """ Return a masked array containing the field data. """ msg_code = self.msg_header['code'] threshold_data = self.prod_descr['threshold_data'] @@ -256,6 +293,11 @@ def get_data(self): scale, offset = np.frombuffer(threshold_data[:8], '>f4') data = (self.raw_data - offset) / (scale) * 0.01 mdata = np.ma.array(data, mask=self.raw_data < 1) + + elif msg_code in [176]: + scale, offset = np.frombuffer(threshold_data[:8], '>f4') + data = (self.raw_data - offset) / (scale) + mdata = np.ma.array(data, mask=self.raw_data < 1) elif msg_code in [165, 177]: # Corresponds to classifications in table on page 3-37 @@ -354,6 +396,171 @@ def nexrad_level3_message_code(filename): # Tables and page number refer to those in this document. +class Level3XDRParser(Unpacker): + """Handle XDR-formatted Level 3 NEXRAD products. + + This class is virtually identical to the Metpy implementation. It has been + pulled into this module to avoid future changes to the Metpy package from + breaking something. The class may be imported from Metpy as a dependency + if someday the project has matured so that features breaking are unlikely. + + This class has been modified from MetPy + Copyright (c) 2009,2015,2016,2017 MetPy Developers. + Distributed under the terms of the BSD 3-Clause License. + SPDX-License-Identifier: BSD-3-Clause + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + POSSIBILITY OF SUCH DAMAGE. + """ + + def __call__(self, packet_code): + """Perform the actual unpacking.""" + xdr = {} + + if packet_code == 28: + xdr.update(self._unpack_prod_desc()) + else: + raise NotImplementedError( + 'Unknown XDR Component: %d' % (packet_code)) + + # Check that we got it all + self.done() + return xdr + + def unpack_string(self): + """Unpack the internal data as a string.""" + return Unpacker.unpack_string(self).decode('ascii') + + def _unpack_prod_desc(self): + xdr = {} + + # NOTE: The ICD incorrectly lists op-mode, vcp, el_num, and + # spare as int*2. Changing to int*4 makes things parse correctly. + xdr['name'] = self.unpack_string() + xdr['description'] = self.unpack_string() + xdr['code'] = self.unpack_int() + xdr['type'] = self.unpack_int() + xdr['prod_time'] = self.unpack_uint() + xdr['radar_name'] = self.unpack_string() + xdr['latitude'] = self.unpack_float() + xdr['longitude'] = self.unpack_float() + xdr['height'] = self.unpack_float() + xdr['vol_time'] = self.unpack_uint() + xdr['el_time'] = self.unpack_uint() + xdr['el_angle'] = self.unpack_float() + xdr['vol_num'] = self.unpack_int() + xdr['op_mode'] = self.unpack_int() + xdr['vcp_num'] = self.unpack_int() + xdr['el_num'] = self.unpack_int() + xdr['compression'] = self.unpack_int() + xdr['uncompressed_size'] = self.unpack_int() + xdr['parameters'] = self._unpack_parameters() + xdr['components'] = self._unpack_components() + + return xdr + + def _unpack_parameters(self): + num = self.unpack_int() + + # ICD documents a "pointer" here, that seems to be garbage. Just read + # and use the number, starting the list immediately. + self.unpack_int() + + if num == 0: + return None + + ret = [] + for i in range(num): + ret.append((self.unpack_string(), self.unpack_string())) + if i < num - 1: + self.unpack_int() # Another pointer for the 'list' ? + + if num == 1: + ret = ret[0] + + return ret + + def _unpack_components(self): + num = self.unpack_int() + + # ICD documents a "pointer" here, that seems to be garbage. Just read + # and use the number, starting the list immediately. + self.unpack_int() + + ret = [] + for i in range(num): + try: + code = self.unpack_int() + ret.append(self._component_lookup[code](self)) + if i < num - 1: + self.unpack_int() # Another pointer for the 'list' ? + except KeyError: + raise NotImplementedError( + 'Unknown XDR Component: %d' % (code)) + break + + if num == 1: + ret = ret[0] + + return ret + + radial_fmt = namedtuple('RadialComponent', ['description', 'gate_width', + 'first_gate', 'parameters', + 'radials']) + radial_data_fmt = namedtuple('RadialData', ['azimuth', 'elevation', 'width', + 'num_bins', 'attributes', + 'data']) + + def _unpack_radial(self): + ret = self.radial_fmt(description=self.unpack_string(), + gate_width=self.unpack_float(), + first_gate=self.unpack_float(), + parameters=self._unpack_parameters(), + radials=None) + num_rads = self.unpack_int() + rads = [] + for _ in range(num_rads): + # ICD is wrong, says num_bins is float, should be int + rads.append(self.radial_data_fmt(azimuth=self.unpack_float(), + elevation=self.unpack_float(), + width=self.unpack_float(), + num_bins=self.unpack_int(), + attributes=self.unpack_string(), + data=self.unpack_array(self.unpack_int))) + return ret._replace(radials=rads) + + text_fmt = namedtuple('TextComponent', ['parameters', 'text']) + + def _unpack_text(self): + return self.text_fmt(parameters=self._unpack_parameters(), + text=self.unpack_string()) + + _component_lookup = {1: _unpack_radial, 4: _unpack_text} + + def _int16_to_float16(val): """ Convert a 16 bit interger into a 16 bit float. """ # NEXRAD Level III float16 format defined on page 3-33. @@ -402,6 +609,7 @@ def _int16_to_float16(val): 173: 1., 174: 1., 175: 1., + 176: 0.25, 177: 0.25, 181: 150., 182: 150., @@ -452,6 +660,7 @@ def _int16_to_float16(val): INT4 = 'i' UINT4 = 'I' REAL4 = 'f' +LONG = 'l' # 3.3.1 Graphic Product Messages @@ -513,13 +722,14 @@ def _int16_to_float16(val): # Display data packets ) +AF1F = -20705 # struct.unpack('>h', 'AF1F'.decode('hex')) +SUPPORTED_PACKET_CODES = [16, AF1F, 28] + # Digital Radial Data Array Packet - Packet Code 16 (Sheet 2) -# Figure 3-11c (Sheet 1 and 2), page 3-120. +# Figure 3-11c (Sheet 1 and 2), page 3-120 # and # Radial Data Packet - Packet Code AF1F -# Figure 3-10 (Sheet 1 and 2), page 3-113. -AF1F = -20705 # struct.unpack('>h', 'AF1F'.decode('hex')) -SUPPORTED_PACKET_CODES = [16, AF1F] # elsewhere +# Figure 3-10 (Sheet 1 and 2), page 3-113 RADIAL_PACKET_HEADER = ( ('packet_code', INT2), # Packet Code, Type 16 ('first_bin', INT2), # Location of first range bin. @@ -536,6 +746,14 @@ def _int16_to_float16(val): ('angle_delta', INT2) # Delta angle from previous radial. ) +# Generic Data Packet - Packet Code 28 +# Figure 3-15c (Sheet 1), page 3-132 +GEN_DATA_PACK_HEADER = ( + ('packet_code', INT2), # Packet Code, Type 28 + ('reserved', INT2), # Reserved for future use. Should be set to 0. + ('num_bytes', LONG), # Number of bytes to follow in this packet +) + # A list of the NEXRAD Level 3 Product supported by this module taken # from the "Message Code for Products" Table III pages 3-15 to 3-22 # All the supported products have a Radial Image Message format. @@ -587,6 +805,8 @@ def _int16_to_float16(val): # Difference Accumulation 175, # Digital Storm Total # Difference Accumulation + 176, # Digital Instantaneous + # Precipitation Rate 177, # Hybrid Hydrometeor # Classification 181, # Base Reflectivity @@ -747,8 +967,6 @@ def _int16_to_float16(val): # # Confidence # 166, # Melting Layer Linked Contour Vectors/ # # Set Color Level -# 176, # Digital Instantaneous Generic Radial Product Format -# # Precipitation Rate # 178-193,# Reserved for Future Products # 196-198,# Reserved for Future Products # 200-210,# Reserved for Future Products diff --git a/pyart/io/tests/test_nexrad_level3.py b/pyart/io/tests/test_nexrad_level3.py index 9c3c1ef06d..416e1700cb 100644 --- a/pyart/io/tests/test_nexrad_level3.py +++ b/pyart/io/tests/test_nexrad_level3.py @@ -65,7 +65,7 @@ def test_nexrad_level3_msg19(): assert 'reflectivity' in radar.fields.keys() assert radar.fields['reflectivity']['data'].shape == (360, 230) assert type(radar.fields['reflectivity']['data']) is MaskedArray - assert round(radar.fields['reflectivity']['data'][10, 10]) == 25. + assert round(radar.fields['reflectivity']['data'][10, 10]) == np.float32(25) def test_nexrad_level3_msg161(): @@ -127,7 +127,7 @@ def test_nexrad_level3_msg161(): assert field_name in radar.fields.keys() assert radar.fields[field_name]['data'].shape == (360, 1200) assert type(radar.fields[field_name]['data']) is MaskedArray - assert round(radar.fields[field_name]['data'][103, 170]) == 2. + assert round(radar.fields[field_name]['data'][103, 170]) == np.float32(2) def test_nexrad_level3_msg161_fileobj(): @@ -138,4 +138,60 @@ def test_nexrad_level3_msg161_fileobj(): assert field_name in radar.fields.keys() assert radar.fields[field_name]['data'].shape == (360, 1200) assert type(radar.fields[field_name]['data']) is MaskedArray - assert round(radar.fields[field_name]['data'][103, 170]) == 2. + assert round(radar.fields[field_name]['data'][103, 170]) == np.float32(2) + +def test_nexrad_level3_msg176(): + radar = pyart.io.read_nexrad_level3(pyart.testing.NEXRAD_LEVEL3_MSG176) + + assert radar.time['units'] == 'seconds since 2020-03-19T18:01:43Z' + assert radar.time['data'].shape == (360, ) + assert round(radar.time['data'][0]) == 0. + + assert radar.range['data'].shape == (920, ) + assert round(radar.range['data'][100]) == 25125. + + assert radar.scan_type == 'ppi' + + assert radar.latitude['data'].shape == (1, ) + assert round(radar.latitude['data'][0]) == 42.0 + + assert radar.longitude['data'].shape == (1, ) + assert round(radar.longitude['data'][0]) == -88.0 + + assert radar.altitude['data'].shape == (1, ) + assert round(radar.altitude['data'][0]) == 232. + + assert radar.altitude_agl is None + + assert radar.sweep_number['data'].shape == (1, ) + assert radar.sweep_number['data'][0] == 0 + + assert radar.sweep_mode['data'].shape == (1, ) + assert np.all(radar.sweep_mode['data'] == [b'azimuth_surveillance']) + + assert radar.sweep_start_ray_index['data'].shape == (1, ) + assert round(radar.sweep_start_ray_index['data'][0]) == 0.0 + + assert radar.sweep_end_ray_index['data'].shape == (1, ) + assert round(radar.sweep_end_ray_index['data'][0]) == 359.0 + + assert radar.target_scan_rate is None + + assert round(radar.azimuth['data'][0]) == 0.0 + assert round(radar.azimuth['data'][10]) == 10.0 + + assert radar.scan_rate is None + assert radar.antenna_transition is None + assert radar.instrument_parameters is None + assert radar.radar_calibration is None + + assert radar.ngates == 920 + assert radar.nrays == 360 + assert radar.nsweeps == 1 + + field_name = 'radar_estimated_rain_rate' + assert field_name in radar.fields.keys() + assert radar.fields[field_name]['data'].shape == (360, 920) + assert type(radar.fields[field_name]['data']) is MaskedArray + assert round(radar.fields[field_name]['data'][14, 14],3) == np.float32(0.244) + diff --git a/pyart/testing/__init__.py b/pyart/testing/__init__.py index 890d493e24..73236cf455 100644 --- a/pyart/testing/__init__.py +++ b/pyart/testing/__init__.py @@ -13,6 +13,7 @@ from .sample_files import NEXRAD_CDM_FILE from .sample_files import NEXRAD_ARCHIVE_MSG31_COMPRESSED_FILE from .sample_files import NEXRAD_LEVEL3_MSG19, NEXRAD_LEVEL3_MSG163 +from .sample_files import NEXRAD_LEVEL3_MSG176 from .sample_objects import make_empty_ppi_radar, make_target_radar from .sample_objects import make_single_ray_radar, make_velocity_aliased_radar from .sample_objects import make_empty_grid diff --git a/pyart/testing/data/example_nexrad_level3_msg176 b/pyart/testing/data/example_nexrad_level3_msg176 new file mode 100644 index 0000000000..819477c5ff Binary files /dev/null and b/pyart/testing/data/example_nexrad_level3_msg176 differ diff --git a/pyart/testing/sample_files.py b/pyart/testing/sample_files.py index 29bcb2fcf3..f97a859fd7 100644 --- a/pyart/testing/sample_files.py +++ b/pyart/testing/sample_files.py @@ -45,6 +45,7 @@ # KBMX_SDUS84_N0KBMX_201501020205 NEXRAD_LEVEL3_MSG19 = os.path.join(DATA_PATH, 'example_nexrad_level3_msg19') NEXRAD_LEVEL3_MSG163 = os.path.join(DATA_PATH, 'example_nexrad_level3_msg163') +NEXRAD_LEVEL3_MSG176 = os.path.join(DATA_PATH, 'example_nexrad_level3_msg176') NEXRAD_CDM_FILE = os.path.join(DATA_PATH, 'example_nexrad_cdm.bz2') UF_FILE = os.path.join(DATA_PATH, 'example_uf_ppi.uf') INTERP_SOUNDE_FILE = os.path.join(DATA_PATH, 'example_interpolatedsonde.cdf')