From 1fca580ed56a2c11d796fae6fd180bb0b837f62a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kai=20M=C3=BChlbauer?= Date: Thu, 1 Aug 2024 15:43:02 +0200 Subject: [PATCH] FIX: nexrad level2 time and missing elevations (#181) * FIX: time issue due to 1-based day count * FIX: restructure retrieval of msg_31, do not try to read missing elevations, even when listed in msg_5 * add history.md entry * add comment on radial_status 5 --- docs/history.md | 2 ++ tests/io/test_io.py | 4 +-- tests/io/test_nexrad_level2.py | 3 ++ xradar/io/backends/nexrad_level2.py | 51 +++++++++++++++++------------ 4 files changed, 37 insertions(+), 23 deletions(-) diff --git a/docs/history.md b/docs/history.md index 4d2de549..cce02e2c 100644 --- a/docs/history.md +++ b/docs/history.md @@ -6,6 +6,8 @@ * ADD: Add Alfonso to citation doc ({pull}`169`) by [@mgrover1](https://github.com/mgrover1). * ENH: Adding global variables and attributes to iris datatree ({pull}`166`) by [@aladinor](https://github.com/aladinor). * FIX: Set fillvalue before applying scale/offset when exporting to odim ({issue}`122`) by [@pavlikp](https://github.com/pavlikp), ({pull}`173`) by [@kmuehlbauer](https://github.com/kmuehlbauer). +* FIX: Nexrad level2 time offset of 1 day, skip reading missing elevations, introduce new radial_status of 5 + ({issue}`180`) by [@ghiggi](https://github.com/ghiggi), ({pull}`180`) by [@kmuehlbauer](https://github.com/kmuehlbauer). * ADD: Reader for Halo Photonics Doppler lidar data by [@rcjackson](https://github.com/rcjackson) ## 0.5.0 (2024-03-28) diff --git a/tests/io/test_io.py b/tests/io/test_io.py index 291919d7..3417c577 100644 --- a/tests/io/test_io.py +++ b/tests/io/test_io.py @@ -940,8 +940,8 @@ def test_open_nexradlevel2_datatree(nexradlevel2_files): assert rvars["volume_number"] == 0 assert rvars["platform_type"] == "fixed" assert rvars["instrument_type"] == "radar" - assert rvars["time_coverage_start"] == "2016-06-02T15:00:25Z" - assert rvars["time_coverage_end"] == "2016-06-02T15:06:06Z" + assert rvars["time_coverage_start"] == "2016-06-01T15:00:25Z" + assert rvars["time_coverage_end"] == "2016-06-01T15:06:06Z" np.testing.assert_almost_equal(rvars["latitude"].values, np.array(33.65414047)) np.testing.assert_almost_equal(rvars["longitude"].values, np.array(-101.81416321)) np.testing.assert_almost_equal(rvars["altitude"].values, np.array(1029)) diff --git a/tests/io/test_nexrad_level2.py b/tests/io/test_nexrad_level2.py index 54307bea..f46842a7 100644 --- a/tests/io/test_nexrad_level2.py +++ b/tests/io/test_nexrad_level2.py @@ -24,6 +24,9 @@ def test_open_nexradlevel2_datatree(nexradlevel2_files): dtree = open_nexradlevel2_datatree(nexradlevel2_files) ds = dtree["sweep_0"].ds assert ds.attrs["instrument_name"] == "KLBB" + assert ds["time"].min() == np.array( + "2016-06-01T15:00:25.232000000", dtype="datetime64[ns]" + ) assert ds["DBZH"].shape == (720, 1832) assert ds["DBZH"].dims == ("azimuth", "range") assert int(ds.sweep_number.values) == 0 diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 74bb8a0b..5a951eba 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -644,18 +644,26 @@ def get_data_header(self): msg_31_header["filepos"] = self.filepos # retrieve data/const headers from msg 31 # check if this is a new sweep - if msg_31_header["radial_status"] in [0, 3]: - # do not run for the first time - if current_sweep > -1: - sweep["record_end"] = self.rh.recnum - 1 - sweep["intermediate_records"] = sweep_intermediate_records - self._data[current_sweep] = sweep - # create new sweep object - sweep = OrderedDict() - _msg_31_header.append(sweep_msg_31_header) - sweep_msg_31_header = [] - sweep_intermediate_records = [] + status = msg_31_header["radial_status"] + if status == 1: + # 1 - intermediate radial + pass + if status in [2, 4]: + # 2 - end of elevation + # 4 - end of volume + sweep["record_end"] = self.rh.recnum + sweep["intermediate_records"] = sweep_intermediate_records + self._data[current_sweep] = sweep + _msg_31_header.append(sweep_msg_31_header) + if status in [0, 3, 5]: + # 0 - start of new elevation + # 3 - start of new volume + # 5 - start of new elevation, last elevation in VCP current_sweep += 1 + # create new sweep object + sweep = OrderedDict() + sweep_msg_31_header = [] + sweep_intermediate_records = [] sweep["record_number"] = self.rh.recnum sweep["filepos"] = self.filepos sweep["msg_31_header"] = msg_31_header @@ -699,11 +707,6 @@ def get_data_header(self): else: sweep_intermediate_records.append(msg_header) - # finalize last sweep - sweep["record_end"] = self.rh.recnum - sweep["intermediate_records"] = sweep_intermediate_records - self._data[current_sweep] = sweep - _msg_31_header.append(sweep_msg_31_header) return data_header, _msg_31_header, _msg_31_data_header def _check_record(self): @@ -1333,7 +1336,8 @@ def open_store_coordinates(self): elevation = np.array([ms["elevation_angle"] for ms in msg_31_header]) # time - date = np.array([ms["collect_date"] for ms in msg_31_header]) * 86400e3 + # date is 1-based (1 -> 1970-01-01T00:00:00), so we need to subtract 1 + date = (np.array([ms["collect_date"] for ms in msg_31_header]) - 1) * 86400e3 milliseconds = np.array([ms["collect_ms"] for ms in msg_31_header]) rtime = date + milliseconds time_prefix = "milli" @@ -1533,10 +1537,15 @@ def open_nexradlevel2_datatree(filename_or_obj, **kwargs): engine = NexradLevel2BackendEntrypoint - ds = [ - xr.open_dataset(filename_or_obj, group=swp, engine=engine, **kwargs) - for swp in sweeps - ] + # todo: only open once! Needs new xarray built in datatree! + ds = [] + for swp in sweeps: + try: + dsx = xr.open_dataset(filename_or_obj, group=swp, engine=engine, **kwargs) + except IndexError: + break + else: + ds.append(dsx) ds.insert(0, xr.Dataset())