From 56a9ca1f42ca23074dc0b1d2a86d51e1bd4eafa2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kai=20M=C3=BChlbauer?= Date: Wed, 13 Mar 2024 17:25:42 +0100 Subject: [PATCH] NEXRAD Level2 structured reader (#158) * NEXRAD structured reader * WIP nexrad * WIP lint * WIP lint * WIP fixture * WIP lint * WIP entrypoint * WIP fix tests * WIP nexrad bz2 compression * MNT CI/open-radar-data * ADD io tests * fix nexradlevel2 tests * add NexradLevel2 notebook, docstring, fix tests/moments * add history.md entry --- ci/notebooktests.yml | 3 +- ci/unittests.yml | 3 +- docs/history.md | 1 + docs/importers.md | 19 + docs/usage.md | 1 + environment.yml | 3 +- examples/notebooks/NexradLevel2.ipynb | 271 +++++ pyproject.toml | 2 + tests/conftest.py | 28 + tests/io/test_io.py | 113 ++ tests/io/test_nexrad_level2.py | 434 +++++++ xradar/io/backends/__init__.py | 2 + xradar/io/backends/common.py | 2 +- xradar/io/backends/iris.py | 12 +- xradar/io/backends/nexrad_level2.py | 1546 +++++++++++++++++++++++++ xradar/util.py | 8 + 16 files changed, 2436 insertions(+), 12 deletions(-) create mode 100644 examples/notebooks/NexradLevel2.ipynb create mode 100644 tests/io/test_nexrad_level2.py create mode 100644 xradar/io/backends/nexrad_level2.py diff --git a/ci/notebooktests.yml b/ci/notebooktests.yml index a0587a09..012e0a14 100644 --- a/ci/notebooktests.yml +++ b/ci/notebooktests.yml @@ -16,6 +16,7 @@ dependencies: - netCDF4 - notebook - numpy + - open-radar-data>=0.1.0 - pip - pyproj - pytest @@ -32,5 +33,3 @@ dependencies: - arm_pyart - ffmpeg - zarr - - pip: - - git+https://github.com/openradar/open-radar-data.git diff --git a/ci/unittests.yml b/ci/unittests.yml index 4ffc7a79..eb9f3f82 100644 --- a/ci/unittests.yml +++ b/ci/unittests.yml @@ -13,6 +13,7 @@ dependencies: - lat_lon_parser - netCDF4 - numpy + - open-radar-data>=0.1.0 - pip - pyproj - pytest @@ -25,5 +26,3 @@ dependencies: - xarray - xarray-datatree>=0.0.10 - xmltodict - - pip: - - git+https://github.com/openradar/open-radar-data.git diff --git a/docs/history.md b/docs/history.md index 4d37690f..c45f2725 100644 --- a/docs/history.md +++ b/docs/history.md @@ -6,6 +6,7 @@ * MNT: restructure odim.py/gamic.py, add test_odim.py/test_gamic.py ({pull}`154`) by [@kmuehlbauer](https://github.com/kmuehlbauer). * MNT: use CODECOV token ({pull}`155`) by [@kmuehlbauer](https://github.com/kmuehlbauer). * MNT: fix path for notebook coverage ({pull}`157`) by [@kmuehlbauer](https://github.com/kmuehlbauer). +* ADD: NEXRAD Level2 structured reader ({pull}`158`) by [@kmuehlbauer](https://github.com/kmuehlbauer) and [@mgrover1](https://github.com/mgrover1). ## 0.4.3 (2024-02-24) diff --git a/docs/importers.md b/docs/importers.md index f87c0eae..9ef93fca 100644 --- a/docs/importers.md +++ b/docs/importers.md @@ -10,6 +10,7 @@ Currently xradar can import: - [](#furuno-scn-and-scnx) - [](#rainbow) - [](#iris-sigmet) +- [](#nexradlevel2) ## CfRadial1 @@ -115,3 +116,21 @@ more functions are applied on that {py:class}`xarray:xarray.Dataset`. With {func}`~xradar.io.backends.iris.open_iris_datatree` all groups (eg. ``1``) are extracted. From that the ``root`` group is processed. Everything is finally added as ParentNodes and ChildNodes to a {py:class}`datatree:datatree.DataTree`. + + +## NexradLevel2 + +### NexradLevel2BackendEntryPoint + +The xarray backend {class}`~xradar.io.backends.nexrad_level2.NexradLevel2BackendEntryPoint` +opens the file with {class}`~xradar.io.backends.nexrad_level2.NexradLevel2Store`. Several +private helper functions are used to conveniently access data and +metadata. Finally, the xarray machinery returns a {py:class}`xarray:xarray.Dataset` +with wanted group (eg. ``0``). Depending on the used backend kwargs several +more functions are applied on that {py:class}`xarray:xarray.Dataset`. + +### open_nexradlevel2_datatree + +With {func}`~xradar.io.backends.nexrad_level2.open_nexradlevel2_datatree` +all groups (eg. ``1``) are extracted. From that the ``root`` group is processed. +Everything is finally added as ParentNodes and ChildNodes to a {py:class}`datatree:datatree.DataTree`. diff --git a/docs/usage.md b/docs/usage.md index e554fe51..540ad1ff 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -32,6 +32,7 @@ notebooks/GAMIC notebooks/Furuno notebooks/Rainbow notebooks/Iris +notebooks/NexradLevel2 notebooks/Read-plot-Sigmet-data-from-AWS notebooks/plot-ppi notebooks/angle_reindexing diff --git a/environment.yml b/environment.yml index c13f203a..b812dac4 100644 --- a/environment.yml +++ b/environment.yml @@ -16,6 +16,7 @@ dependencies: - h5py - lat_lon_parser - netCDF4 + - open-radar-data>=0.1.0 - pyproj - xarray - xarray-datatree>=0.0.10 @@ -29,5 +30,3 @@ dependencies: - arm_pyart - ffmpeg - zarr - - pip: - - git+https://github.com/openradar/open-radar-data.git diff --git a/examples/notebooks/NexradLevel2.ipynb b/examples/notebooks/NexradLevel2.ipynb new file mode 100644 index 00000000..a4cd1cb6 --- /dev/null +++ b/examples/notebooks/NexradLevel2.ipynb @@ -0,0 +1,271 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# NEXRAD Level 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import xradar as xd\n", + "import cmweather\n", + "from open_radar_data import DATASETS" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Download\n", + "\n", + "Fetching NEXRAD Level2 radar data file from [open-radar-data](https://github.com/openradar/open-radar-data) repository." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "filename = DATASETS.fetch(\"KLBB20160601_150025_V06\")" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## xr.open_dataset\n", + "\n", + "Making use of the xarray `nexradlevel2` backend. We also need to provide the group. Note, that we are using CfRadial2 group access pattern." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(filename, group=\"sweep_0\", engine=\"nexradlevel2\")\n", + "display(ds)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Plot Time vs. Azimuth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "ds.azimuth.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "### Plot Range vs. Time\n", + "\n", + "We need to sort by time and specify the y-coordinate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "ds.DBZH.sortby(\"time\").plot(y=\"time\", cmap=\"HomeyerRainbow\")" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "### Plot Range vs. Azimuth\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "ds.DBZH.plot(cmap=\"HomeyerRainbow\")" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## backend_kwargs\n", + "\n", + "Beside `first_dim` there are several additional backend_kwargs for the nexradlevel2 backend, which handle different aspects of angle alignment. This comes into play, when azimuth and/or elevation arrays are not evenly spacend and other issues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "help(xd.io.NexradLevel2BackendEntrypoint)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(filename, group=\"sweep_0\", engine=\"nexradlevel2\", first_dim=\"time\")\n", + "display(ds)" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## open_nexradlevel2_datatree\n", + "\n", + "The same works analoguous with the datatree loader. But additionally we can provide a sweep string, number or list." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "help(xd.io.open_nexradlevel2_datatree)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "dtree = xd.io.open_nexradlevel2_datatree(filename, sweep=4)\n", + "display(dtree)" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### Plot Sweep Range vs. Time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "dtree[\"sweep_0\"].ds.DBZH.sortby(\"time\").plot(y=\"time\", cmap=\"HomeyerRainbow\")" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "### Plot Sweep Range vs. Azimuth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "dtree[\"sweep_0\"].ds.DBZH.plot(cmap=\"HomeyerRainbow\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "dtree = xd.io.open_nexradlevel2_datatree(filename, sweep=\"sweep_8\")\n", + "display(dtree)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "dtree = xd.io.open_nexradlevel2_datatree(filename, sweep=[0, 1, 8])\n", + "display(dtree)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "dtree = xd.io.open_nexradlevel2_datatree(\n", + " filename, sweep=[\"sweep_1\", \"sweep_2\", \"sweep_8\"]\n", + ")\n", + "display(dtree)" + ] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 1c91165a..5f022ead 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,8 @@ gamic = "xradar.io.backends:GamicBackendEntrypoint" iris = "xradar.io.backends:IrisBackendEntrypoint" odim = "xradar.io.backends:OdimBackendEntrypoint" rainbow = "xradar.io.backends:RainbowBackendEntrypoint" +nexradlevel2 = "xradar.io.backends:NexradLevel2BackendEntrypoint" + [build-system] requires = [ diff --git a/tests/conftest.py b/tests/conftest.py index b1a0ca11..0a251e33 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -58,3 +58,31 @@ def iris0_file(): @pytest.fixture(scope="session") def iris1_file(): return DATASETS.fetch("SUR210819000227.RAWKPJV") + + +@pytest.fixture(scope="session") +def nexradlevel2_file(): + return DATASETS.fetch("KATX20130717_195021_V06") + + +@pytest.fixture(scope="session") +def nexradlevel2_gzfile(): + fnamei = DATASETS.fetch("KLBB20160601_150025_V06.gz") + fnameo = f"{fnamei[:-3]}_gz" + import gzip + import shutil + + with gzip.open(fnamei) as fin: + with open(fnameo, "wb") as fout: + shutil.copyfileobj(fin, fout) + return fnameo + + +@pytest.fixture(scope="session") +def nexradlevel2_bzfile(): + return DATASETS.fetch("KLBB20160601_150025_V06") + + +@pytest.fixture +def nexradlevel2_files(request): + return request.getfixturevalue(request.param) diff --git a/tests/io/test_io.py b/tests/io/test_io.py index 98b9bb47..1d41b547 100644 --- a/tests/io/test_io.py +++ b/tests/io/test_io.py @@ -19,6 +19,7 @@ open_cfradial1_datatree, open_gamic_datatree, open_iris_datatree, + open_nexradlevel2_datatree, open_odim_datatree, open_rainbow_datatree, ) @@ -811,3 +812,115 @@ def test_cfradfial2_roundtrip(cfradial1_file, first_dim): assert "time" in dtree1[d1].dims assert "time" in dtree2[d2].dims xr.testing.assert_equal(dtree1[d1].ds, dtree2[d2].ds) + + +@pytest.mark.parametrize("sweep", ["sweep_0", 0, [0, 1], ["sweep_0", "sweep_1"]]) +@pytest.mark.parametrize( + "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True +) +def test_open_nexradlevel2_datatree_sweep(nexradlevel2_files, sweep): + dtree = open_nexradlevel2_datatree(nexradlevel2_files, sweep=sweep) + if isinstance(sweep, (str, int)): + lswp = len([sweep]) + else: + lswp = len(sweep) + assert len(dtree.groups[1:]) == lswp + + +@pytest.mark.parametrize( + "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True +) +def test_open_nexradlevel2_datatree(nexradlevel2_files): + dtree = open_nexradlevel2_datatree(nexradlevel2_files) + + # root_attrs + attrs = dtree.attrs + assert attrs["Conventions"] == "None" + assert attrs["instrument_name"] == "KLBB" + + # root vars + rvars = dtree.data_vars + 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" + 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)) + + # iterate over subgroups and check some values + moments = [ + ["DBZH", "PHIDP", "RHOHV", "ZDR"], + ["DBZH", "WRADH", "VRADH"], + ["DBZH", "PHIDP", "RHOHV", "ZDR"], + ["DBZH", "WRADH", "VRADH"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ["DBZH", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"], + ] + elevations = [ + 0.5, + 0.5, + 1.5, + 1.5, + 2.4, + 3.4, + 4.3, + 6.0, + 9.9, + 14.6, + 19.5, + ] + azimuths = [ + 720, + 720, + 720, + 720, + 360, + 360, + 360, + 360, + 360, + 360, + 360, + ] + ranges = [ + 1832, + 1192, + 1632, + 1192, + 1312, + 1076, + 908, + 696, + 448, + 308, + 232, + ] + assert len(dtree.groups[1:]) == 11 + for i, grp in enumerate(dtree.groups[1:]): + print(i) + ds = dtree[grp].ds + assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]} + assert set(ds.data_vars) & ( + sweep_dataset_vars | non_standard_sweep_dataset_vars + ) == set(moments[i]) + assert set(ds.data_vars) & (required_sweep_metadata_vars) == set( + required_sweep_metadata_vars ^ {"azimuth", "elevation"} + ) + assert set(ds.coords) == { + "azimuth", + "elevation", + "time", + "latitude", + "longitude", + "altitude", + "range", + } + assert np.round(ds.elevation.mean().values.item(), 1) == elevations[i] + assert ds.sweep_number.values == int(grp[7:]) diff --git a/tests/io/test_nexrad_level2.py b/tests/io/test_nexrad_level2.py new file mode 100644 index 00000000..ce3f8ee0 --- /dev/null +++ b/tests/io/test_nexrad_level2.py @@ -0,0 +1,434 @@ +#!/usr/bin/env python +# Copyright (c) 2024, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +"""Tests for `xradar.io.nexrad_archive` module.""" + +from collections import OrderedDict + +import pytest +import xarray + +from xradar.io.backends.nexrad_level2 import ( + NexradLevel2BackendEntrypoint, + NEXRADLevel2File, + open_nexradlevel2_datatree, +) + + +@pytest.mark.parametrize( + "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True +) +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["DBZH"].shape == (720, 1832) + assert ds["DBZH"].dims == ("azimuth", "range") + assert int(ds.sweep_number.values) == 0 + + +@pytest.mark.parametrize( + "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True +) +def test_open_nexrad_level2_backend(nexradlevel2_files): + with NEXRADLevel2File(nexradlevel2_files, loaddata=False) as nex: + nsweeps = nex.msg_5["number_elevation_cuts"] + sweeps = [f"sweep_{i}" for i in range(nsweeps)] + assert nsweeps == 11 + for i, group in enumerate(sweeps): + ds = xarray.open_dataset( + nexradlevel2_files, engine=NexradLevel2BackendEntrypoint, group=group + ) + assert ds.attrs["instrument_name"] == "KLBB" + assert int(ds.sweep_number.values) == i + + +@pytest.mark.parametrize( + "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True +) +def test_open_nexradlevel2_file(nexradlevel2_files): + fh = NEXRADLevel2File(nexradlevel2_files) + + # volume header + assert fh.volume_header["tape"] == b"AR2V0006." + assert fh.volume_header["extension"] == b"736" + assert fh.volume_header["date"] == 977403904 + assert fh.volume_header["time"] == 274675715 + assert fh.volume_header["icao"] == b"KLBB" + + # meta_header 15 + assert len(fh.meta_header["msg_15"]) == 5 + assert fh.meta_header["msg_15"][0]["size"] == 1208 + assert fh.meta_header["msg_15"][0]["channels"] == 8 + assert fh.meta_header["msg_15"][0]["type"] == 15 + assert fh.meta_header["msg_15"][0]["seq_id"] == 5 + assert fh.meta_header["msg_15"][0]["date"] == 16940 + assert fh.meta_header["msg_15"][0]["ms"] == 65175863 + assert fh.meta_header["msg_15"][0]["segments"] == 5 + assert fh.meta_header["msg_15"][0]["seg_num"] == 1 + assert fh.meta_header["msg_15"][0]["record_number"] == 0 + # meta_header 13 + assert len(fh.meta_header["msg_13"]) == 49 + assert fh.meta_header["msg_13"][0]["size"] == 1208 + assert fh.meta_header["msg_13"][0]["channels"] == 8 + assert fh.meta_header["msg_13"][0]["type"] == 13 + assert fh.meta_header["msg_13"][0]["seq_id"] == 15690 + assert fh.meta_header["msg_13"][0]["date"] == 16954 + assert fh.meta_header["msg_13"][0]["ms"] == 54016980 + assert fh.meta_header["msg_13"][0]["segments"] == 49 + assert fh.meta_header["msg_13"][0]["seg_num"] == 1 + assert fh.meta_header["msg_13"][0]["record_number"] == 77 + # meta header 18 + assert len(fh.meta_header["msg_18"]) == 4 + assert fh.meta_header["msg_18"][0]["size"] == 1208 + assert fh.meta_header["msg_18"][0]["channels"] == 8 + assert fh.meta_header["msg_18"][0]["type"] == 18 + assert fh.meta_header["msg_18"][0]["seq_id"] == 6 + assert fh.meta_header["msg_18"][0]["date"] == 16940 + assert fh.meta_header["msg_18"][0]["ms"] == 65175864 + assert fh.meta_header["msg_18"][0]["segments"] == 4 + assert fh.meta_header["msg_18"][0]["seg_num"] == 1 + assert fh.meta_header["msg_18"][0]["record_number"] == 126 + # meta header 3 + assert len(fh.meta_header["msg_3"]) == 1 + assert fh.meta_header["msg_3"][0]["size"] == 488 + assert fh.meta_header["msg_3"][0]["channels"] == 8 + assert fh.meta_header["msg_3"][0]["type"] == 3 + assert fh.meta_header["msg_3"][0]["seq_id"] == 15694 + assert fh.meta_header["msg_3"][0]["date"] == 16954 + assert fh.meta_header["msg_3"][0]["ms"] == 54025336 + assert fh.meta_header["msg_3"][0]["segments"] == 1 + assert fh.meta_header["msg_3"][0]["seg_num"] == 1 + assert fh.meta_header["msg_3"][0]["record_number"] == 131 + # meta header 5 + assert len(fh.meta_header["msg_5"]) == 1 + assert fh.meta_header["msg_5"][0]["size"] == 272 + assert fh.meta_header["msg_5"][0]["channels"] == 8 + assert fh.meta_header["msg_5"][0]["type"] == 5 + assert fh.meta_header["msg_5"][0]["seq_id"] == 15695 + assert fh.meta_header["msg_5"][0]["date"] == 16954 + assert fh.meta_header["msg_5"][0]["ms"] == 54025336 + assert fh.meta_header["msg_5"][0]["segments"] == 1 + assert fh.meta_header["msg_5"][0]["seg_num"] == 1 + assert fh.meta_header["msg_5"][0]["record_number"] == 132 + assert fh.msg_5 == OrderedDict( + [ + ("message_size", 264), + ("pattern_type", 2), + ("pattern_number", 21), + ("number_elevation_cuts", 11), + ("clutter_map_group_number", 1), + ("doppler_velocity_resolution", 2), + ("pulse_width", 2), + ( + "elevation_data", + [ + OrderedDict( + [ + ("elevation_angle", 0.4833984375), + ("channel_config", 0), + ("waveform_type", 1), + ("super_resolution", 11), + ("prf_number", 1), + ("prf_pulse_count", 28), + ("azimuth_rate", 8256), + ("ref_thresh", 16), + ("vel_thresh", 16), + ("sw_thresh", 16), + ("zdr_thres", 16), + ("phi_thres", 16), + ("rho_thres", 16), + ("edge_angle_1", 0), + ("dop_prf_num_1", 0), + ("dop_prf_pulse_count_1", 0), + ("edge_angle_2", 0), + ("dop_prf_num_2", 0), + ("dop_prf_pulse_count_2", 0), + ("edge_angle_3", 0), + ("dop_prf_num_3", 0), + ("dop_prf_pulse_count_3", 0), + ] + ), + OrderedDict( + [ + ("elevation_angle", 0.4833984375), + ("channel_config", 0), + ("waveform_type", 2), + ("super_resolution", 7), + ("prf_number", 0), + ("prf_pulse_count", 0), + ("azimuth_rate", 8272), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 75), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 75), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 75), + ] + ), + OrderedDict( + [ + ("elevation_angle", 1.4501953125), + ("channel_config", 0), + ("waveform_type", 1), + ("super_resolution", 11), + ("prf_number", 1), + ("prf_pulse_count", 28), + ("azimuth_rate", 8256), + ("ref_thresh", 16), + ("vel_thresh", 16), + ("sw_thresh", 16), + ("zdr_thres", 16), + ("phi_thres", 16), + ("rho_thres", 16), + ("edge_angle_1", 0), + ("dop_prf_num_1", 0), + ("dop_prf_pulse_count_1", 0), + ("edge_angle_2", 0), + ("dop_prf_num_2", 0), + ("dop_prf_pulse_count_2", 0), + ("edge_angle_3", 0), + ("dop_prf_num_3", 0), + ("dop_prf_pulse_count_3", 0), + ] + ), + OrderedDict( + [ + ("elevation_angle", 1.4501953125), + ("channel_config", 0), + ("waveform_type", 2), + ("super_resolution", 7), + ("prf_number", 0), + ("prf_pulse_count", 0), + ("azimuth_rate", 8272), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 75), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 75), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 75), + ] + ), + OrderedDict( + [ + ("elevation_angle", 2.4169921875), + ("channel_config", 0), + ("waveform_type", 4), + ("super_resolution", 14), + ("prf_number", 2), + ("prf_pulse_count", 8), + ("azimuth_rate", 8144), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 59), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 59), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 59), + ] + ), + OrderedDict( + [ + ("elevation_angle", 3.3837890625), + ("channel_config", 0), + ("waveform_type", 4), + ("super_resolution", 14), + ("prf_number", 2), + ("prf_pulse_count", 8), + ("azimuth_rate", 8144), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 59), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 59), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 59), + ] + ), + OrderedDict( + [ + ("elevation_angle", 4.306640625), + ("channel_config", 0), + ("waveform_type", 4), + ("super_resolution", 14), + ("prf_number", 2), + ("prf_pulse_count", 8), + ("azimuth_rate", 8144), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 59), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 59), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 59), + ] + ), + OrderedDict( + [ + ("elevation_angle", 6.0205078125), + ("channel_config", 0), + ("waveform_type", 4), + ("super_resolution", 14), + ("prf_number", 3), + ("prf_pulse_count", 12), + ("azimuth_rate", 8144), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 4), + ("dop_prf_pulse_count_1", 59), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 4), + ("dop_prf_pulse_count_2", 59), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 4), + ("dop_prf_pulse_count_3", 59), + ] + ), + OrderedDict( + [ + ("elevation_angle", 9.8876953125), + ("channel_config", 0), + ("waveform_type", 3), + ("super_resolution", 10), + ("prf_number", 0), + ("prf_pulse_count", 0), + ("azimuth_rate", 10384), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 7), + ("dop_prf_pulse_count_1", 82), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 7), + ("dop_prf_pulse_count_2", 82), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 7), + ("dop_prf_pulse_count_3", 82), + ] + ), + OrderedDict( + [ + ("elevation_angle", 14.58984375), + ("channel_config", 0), + ("waveform_type", 3), + ("super_resolution", 10), + ("prf_number", 0), + ("prf_pulse_count", 0), + ("azimuth_rate", 10432), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 7), + ("dop_prf_pulse_count_1", 82), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 7), + ("dop_prf_pulse_count_2", 82), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 7), + ("dop_prf_pulse_count_3", 82), + ] + ), + OrderedDict( + [ + ("elevation_angle", 19.51171875), + ("channel_config", 0), + ("waveform_type", 3), + ("super_resolution", 10), + ("prf_number", 0), + ("prf_pulse_count", 0), + ("azimuth_rate", 10496), + ("ref_thresh", 28), + ("vel_thresh", 28), + ("sw_thresh", 28), + ("zdr_thres", 28), + ("phi_thres", 28), + ("rho_thres", 28), + ("edge_angle_1", 5464), + ("dop_prf_num_1", 7), + ("dop_prf_pulse_count_1", 82), + ("edge_angle_2", 38232), + ("dop_prf_num_2", 7), + ("dop_prf_pulse_count_2", 82), + ("edge_angle_3", 60984), + ("dop_prf_num_3", 7), + ("dop_prf_pulse_count_3", 82), + ] + ), + ], + ), + ] + ) + + # meta header 2 + assert len(fh.meta_header["msg_2"]) == 1 + assert fh.meta_header["msg_2"][0]["size"] == 48 + assert fh.meta_header["msg_2"][0]["channels"] == 8 + assert fh.meta_header["msg_2"][0]["type"] == 2 + assert fh.meta_header["msg_2"][0]["seq_id"] == 15693 + assert fh.meta_header["msg_2"][0]["date"] == 16954 + assert fh.meta_header["msg_2"][0]["ms"] == 54025336 + assert fh.meta_header["msg_2"][0]["segments"] == 1 + assert fh.meta_header["msg_2"][0]["seg_num"] == 1 + assert fh.meta_header["msg_2"][0]["record_number"] == 133 + + # data header + assert len(fh.data_header) == 5400 + msg_31_header_length = [720, 720, 720, 720, 360, 360, 360, 360, 360, 360, 360] + for i, head in enumerate(fh.msg_31_header): + assert len(head) == msg_31_header_length[i] diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 4edf3d6d..8c1b2e19 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -15,6 +15,7 @@ .. automodule:: xradar.io.backends.furuno .. automodule:: xradar.io.backends.rainbow .. automodule:: xradar.io.backends.iris +.. automodule:: xradar.io.backends.nexrad_level2 """ @@ -24,5 +25,6 @@ from .iris import * # noqa from .odim import * # noqa from .rainbow import * # noqa +from .nexrad_level2 import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/common.py b/xradar/io/backends/common.py index d954bce5..25385e6f 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -118,6 +118,7 @@ def _assign_root(sweeps): # assign root attributes attrs = {} attrs["Conventions"] = sweeps[0].attrs.get("Conventions", "None") + attrs["instrument_name"] = sweeps[0].attrs.get("instrument_name", "None") attrs.update( { "version": "None", @@ -127,7 +128,6 @@ def _assign_root(sweeps): "source": "None", "history": "None", "comment": "im/exported using xradar", - "instrument_name": "None", } ) root = root.assign_attrs(attrs) diff --git a/xradar/io/backends/iris.py b/xradar/io/backends/iris.py index 18701416..07a5fdbe 100644 --- a/xradar/io/backends/iris.py +++ b/xradar/io/backends/iris.py @@ -271,8 +271,8 @@ def decode_string(data): UINT1 = {"fmt": "B", "dtype": "unit8"} UINT2 = {"fmt": "H", "dtype": "uint16"} UINT4 = {"fmt": "I", "dtype": "unint32"} -FLT4 = {"fmt": "f"} -FLT8 = {"fmt": "d"} +FLT4 = {"fmt": "f", "dtype": "float32"} +FLT8 = {"fmt": "d", "dtype": "float64"} BIN1 = { "name": "BIN1", "dtype": "uint8", @@ -358,7 +358,7 @@ def _get_struct_dtype(dictionary): return np.dtype(dtypes) -def _unpack_dictionary(buffer, dictionary, rawdata=False): +def _unpack_dictionary(buffer, dictionary, rawdata=False, byte_order="<"): """Unpacks binary data using the given dictionary structure. Parameters @@ -373,7 +373,7 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): Ordered Dictionary with unpacked data """ # get format and substructures of dictionary - fmt, sub = _get_fmt_string(dictionary, retsub=True) + fmt, sub = _get_fmt_string(dictionary, retsub=True, byte_order=byte_order) # unpack into OrderedDict data = OrderedDict(zip(dictionary, struct.unpack(fmt, buffer))) @@ -398,7 +398,9 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): pass # unpack sub dictionary try: - data[k] = _unpack_dictionary(data[k], v, rawdata=rawdata) + data[k] = _unpack_dictionary( + data[k], v, rawdata=rawdata, byte_order=byte_order + ) except TypeError: pass diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py new file mode 100644 index 00000000..a2a0e4d1 --- /dev/null +++ b/xradar/io/backends/nexrad_level2.py @@ -0,0 +1,1546 @@ +#!/usr/bin/env python +# Copyright (c) 2024, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +""" +NEXRAD Level2 Data I/O +^^^^^^^^^^^^^^^^^^^^^^ + +Reads data from NEXRAD Level2 data format + +See https://www.roc.noaa.gov/WSR88D/BuildInfo/Files.aspx + +Documents: + - ICD 2620002 M ICD FOR RDA/RPG - Build RDA 11.5/RPG 13.0 (PDF) + - ICD 2620010 E ICD FOR ARCHIVE II/USER - Build 12.0 (PDF) + +To read from NEXRAD Level2 files :class:`numpy:numpy.memmap` is used for +uncompressed files (pre 2016-06-01) and :class:`bz2:BZ2Decompressor` for BZ2 +compressed data. The NEXRAD header (`VOLUME_HEADER`, `MSG_HEADER`) are read in +any case into dedicated OrderedDict's. Reading sweep data can be skipped by +setting `loaddata=False`. By default, the data is decoded on the fly. +Using `rawdata=True` the data will be kept undecoded. + +Code adapted from Py-ART. + +.. autosummary:: + :nosignatures: + :toctree: generated/ + + {} +""" +__all__ = [ + "NexradLevel2BackendEntrypoint", + "open_nexradlevel2_datatree", +] +__doc__ = __doc__.format("\n ".join(__all__)) + +import bz2 +import struct +from collections import OrderedDict + +import numpy as np +import xarray as xr +from datatree import DataTree +from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint +from xarray.backends.file_manager import CachingFileManager +from xarray.backends.store import StoreBackendEntrypoint +from xarray.core import indexing +from xarray.core.utils import FrozenDict +from xarray.core.variable import Variable + +from xradar import util +from xradar.io.backends.common import ( + _assign_root, + _attach_sweep_groups, +) +from xradar.model import ( + get_altitude_attrs, + get_azimuth_attrs, + get_elevation_attrs, + get_latitude_attrs, + get_longitude_attrs, + get_range_attrs, + get_time_attrs, + moment_attrs, + sweep_vars_mapping, +) + +from .iris import ( + BIN2, + IrisRecord, + _get_fmt_string, + _unpack_dictionary, + string_dict, +) + +#: mapping from NEXRAD names to CfRadial2/ODIM +nexrad_mapping = { + "REF": "DBZH", + "VEL": "VRADH", + "SW ": "WRADH", + "ZDR": "ZDR", + "PHI": "PHIDP", + "RHO": "RHOHV", + "CFP": "CCORH", +} + +# NEXRAD Level II file structures and sizes +# The deails on these structures are documented in: +# "Interface Control Document for the Achive II/User" RPG Build 12.0 +# Document Number 2620010E +# and +# "Interface Control Document for the RDA/RPG" Open Build 13.0 +# Document Number 2620002M +# Tables and page number refer to those in the second document unless +# otherwise noted. +RECORD_BYTES = 2432 +COMPRESSION_RECORD_SIZE = 12 +CONTROL_WORD_SIZE = 4 + +# format of structure elements +# section 3.2.1, page 3-2 +UINT1 = {"fmt": "B", "dtype": "unit8"} +UINT2 = {"fmt": "H", "dtype": "uint16"} +UINT4 = {"fmt": "I", "dtype": "uint32"} +SINT1 = {"fmt": "b", "dtype": "int8"} +SINT2 = {"fmt": "h", "dtype": "int16"} +SINT4 = {"fmt": "i", "dtype": "int32"} +FLT4 = {"fmt": "f", "dtype": "float32"} +CODE1 = {"fmt": "B"} +CODE2 = {"fmt": "H"} + + +class NEXRADFile: + """ + Class for accessing data in a NEXRAD (WSR-88D) Level II file. + + NEXRAD Level II files [1]_, also know as NEXRAD Archive Level II or + WSR-88D Archive level 2, are available from the NOAA National Climate Data + Center [2]_ as well as on the UCAR THREDDS Data Server [3]_. Files with + uncompressed messages and compressed messages are supported. This class + supports reading both "message 31" and "message 1" type files. + + Parameters + ---------- + filename : str + Filename of Archive II file to read. + + References + ---------- + .. [1] http://www.roc.noaa.gov/WSR88D/Level_II/Level2Info.aspx + .. [2] http://www.ncdc.noaa.gov/ + .. [3] http://thredds.ucar.edu/thredds/catalog.html + + """ + + def __init__(self, filename, mode="r", loaddata=False): + """initalize the object.""" + self._fp = None + self._filename = filename + # read in the volume header and compression_record + if hasattr(filename, "read"): + self._fh = filename + else: + self._fp = open(filename, "rb") + self._fh = np.memmap(self._fp, mode=mode) + self._filepos = 0 + self._rawdata = False + self._loaddata = loaddata + self._bz2_indices = None + self.volume_header = self.get_header(VOLUME_HEADER) + return + + @property + def filename(self): + return self._filename + + @property + def fh(self): + return self._fh + + @property + def filepos(self): + return self._filepos + + @filepos.setter + def filepos(self, pos): + self._filepos = pos + + def read_from_file(self, size): + """Read from file. + + Parameters + ---------- + size : int + Number of data words to read. + + Returns + ------- + data : array-like + numpy array of data + """ + start = self._filepos + self._filepos += size + return self._fh[start : self._filepos] + + def get_header(self, header): + len = struct.calcsize(_get_fmt_string(header)) + head = _unpack_dictionary(self.read_from_file(len), header, self._rawdata) + return head + + @property + def bz2_record_indices(self): + """Get file offsets of bz2 records.""" + if self._bz2_indices is None: + # magic number inside BZ2 + seq = np.array([49, 65, 89, 38, 83, 89], dtype=np.uint8) + rd = util.rolling_dim(self._fh, 6) + self._bz2_indices = np.nonzero((rd == seq).all(1))[0] - 8 + return self._bz2_indices + + @property + def is_compressed(self): + """File contains bz2 compressed data.""" + size = self._fh[24:28].view(dtype=">u4")[0] + return size > 0 + + def close(self): + if self._fp is not None: + self._fp.close() + + __del__ = close + + def __enter__(self): + return self + + def __exit__(self, type, value, traceback): + self.close() + + +class NEXRADRecordFile(NEXRADFile): + """NEXRADRecordFile Record File class""" + + def __init__(self, filename, **kwargs): + super().__init__(filename=filename, **kwargs) + self._rh = None + self._rc = None + self._ldm = dict() + self._record_number = None + + @property + def rh(self): + """Returns current record object.""" + return self._rh + + @rh.setter + def rh(self, value): + """Sets current record object.""" + self._rh = value + + @property + def record_number(self): + """Returns current record number.""" + return self._record_number + + @record_number.setter + def record_number(self, value): + """Sets current record number.""" + self._record_number = value + + @property + def record_size(self): + """Returns current record number.""" + return self._record_size + + @record_size.setter + def record_size(self, value): + """Sets current record number.""" + self._record_size = value + + def _check_record(self): + """Checks record for correct size. + + Need to be implemented in the derived classes + """ + return True + + def get_end(self, buf): + if len(buf) < 16: + return 0 + header = _unpack_dictionary(buf, MSG_HEADER, self._rawdata, byte_order=">") + size = header["size"] * 2 + 12 + if header["type"] != 31: + size = size if size >= RECORD_BYTES else RECORD_BYTES + return size + + def init_record(self, recnum): + """Initialize record using given number.""" + + # map record numbers to ldm compressed records + def get_ldm(recnum): + if recnum < 134: + return 0 + mod = ((recnum - 134) // 120) + 1 + return mod + + if self.is_compressed: + ldm = get_ldm(recnum) + # get uncompressed ldm record + if self._ldm.get(ldm, None) is None: + # otherwise extract wanted ldm compressed record + if ldm >= len(self.bz2_record_indices): + return False + start = self.bz2_record_indices[ldm] + size = self._fh[start : start + 4].view(dtype=">u4")[0] + self._fp.seek(start + 4) + dec = bz2.BZ2Decompressor() + self._ldm[ldm] = np.frombuffer( + dec.decompress(self._fp.read(size)), dtype=np.uint8 + ) + + # rectrieve wanted record and put into self.rh + if recnum < 134: + start = recnum * RECORD_BYTES + if not self.is_compressed: + start += 24 + stop = start + RECORD_BYTES + else: + if self.is_compressed: + # get index into current compressed ldm record + rnum = (recnum - 134) % 120 + start = self.record_size + self.filepos if rnum else 0 + buf = self._ldm[ldm][start + 12 : start + 12 + LEN_MSG_HEADER] + size = self.get_end(buf) + if not size: + return False + stop = start + size + else: + start = self.record_size + self.filepos + buf = self.fh[start + 12 : start + 12 + LEN_MSG_HEADER] + size = self.get_end(buf) + if not size: + return False + stop = start + size + self.record_number = recnum + self.record_size = stop - start + if self.is_compressed: + self.rh = IrisRecord(self._ldm[ldm][start:stop], recnum) + else: + self.rh = IrisRecord(self.fh[start:stop], recnum) + self.filepos = start + return self._check_record() + + def init_record_by_filepos(self, recnum, filepos): + """Initialize record using given record number and position.""" + start = filepos + buf = self.fh[start + 12 : start + 12 + LEN_MSG_HEADER] + size = self.get_end(buf) + if not size: + return False + stop = start + size + self.record_number = recnum + self.record_size = stop - start + self.rh = IrisRecord(self.fh[start:stop], recnum) + self.filepos = start + return self._check_record() + + def init_next_record(self): + """Get next record from file. + + This increases record_number count and initialises a new IrisRecord + with the calculated start and stop file offsets. + + Returns + ------- + chk : bool + True, if record is truncated. + """ + return self.init_record(self.record_number + 1) + + def array_from_record(self, words, width, dtype): + """Retrieve array from current record. + + Parameters + ---------- + words : int + Number of data words to read. + width : int + Size of the data word to read in bytes. + dtype : str + dtype string specifying data format. + + Returns + ------- + data : array-like + numpy array of data + """ + return self.rh.read(words, width=width).view(dtype=dtype) + + def bytes_from_record(self, words, width): + """Retrieve bytes from current record. + + Parameters + ---------- + words : int + Number of data words to read. + width : int + Size of the data word to read in bytes. + + Returns + ------- + data : array-like + numpy array of data + """ + return self.rh.read(words, width=width) + + +class NEXRADLevel2File(NEXRADRecordFile): + def __init__(self, filename, **kwargs): + super().__init__(filename, **kwargs) + + # check compression + + self._data_header = None + # get all metadata headers + # message 15, 13, 18, 3, 5 and 2 + self._meta_header = None + + # message 15 + # RDA Clutter Map Data + # message 13 + # RDA Clutter Filter Bypass Map + # message 18 + # RDA Adaptable Parameters + # message 3 + # RDA Performance/Maintenance Data + # message 5 + # RDA Volume Coverage Data + self._msg_5_data = None + # message 2 + # RDA Status Data + + # message 31 headers + # Digital Radar Data Generic Format + self._data_header = None + self._msg_31_header = None + self._msg_31_data_header = None + self._data = OrderedDict() + + @property + def data_header(self): + """Retrieve data header.""" + if self._data_header is None: + ( + self._data_header, + self._msg_31_header, + self._msg_31_data_header, + ) = self.get_data_header() + return self._data_header + + @property + def msg_31_header(self): + """Retrieve MSG31 Header.""" + if self._msg_31_header is None: + ( + self._data_header, + self._msg_31_header, + self._msg_31_data_header, + ) = self.get_data_header() + return self._msg_31_header + + @property + def msg_31_data_header(self): + """Retrieve MSG31 data header.""" + if self._msg_31_data_header is None: + ( + self._data_header, + self._msg_31_header, + self._msg_31_data_header, + ) = self.get_data_header() + return self._msg_31_data_header + + @property + def meta_header(self): + """Retrieve metadat header.""" + if self._meta_header is None: + self._meta_header = self.get_metadata_header() + return self._meta_header + + @property + def msg_5(self): + """Retrieve MSG5 data.""" + if self._msg_5_data is None: + self._msg_5_data = self.get_msg_5_data() + return self._msg_5_data + + @property + def data(self): + """Retrieve data.""" + return self._data + + def get_sweep(self, sweep_number, moments=None): + """Retrieve sweep according sweep_number.""" + if moments is None: + moments = self.msg_31_data_header[sweep_number]["msg_31_data_header"].keys() + # get selected coordinate type names + selected = [k for k in DATA_BLOCK_CONSTANT_IDENTIFIERS.keys() if k in moments] + for moment in selected: + self.get_moment(sweep_number, moment, "sweep_constant_data") + + # get selected data type names + selected = [k for k in DATA_BLOCK_VARIABLE_IDENTIFIERS.keys() if k in moments] + for moment in selected: + self.get_moment(sweep_number, moment, "sweep_data") + + def get_moment(self, sweep_number, moment, mtype): + """Retrieve Moment according sweep_number and moment.""" + sweep = self.data[sweep_number] + # create new sweep_data OrderedDict + if not sweep.get(mtype, False): + sweep[mtype] = OrderedDict() + # check if moment exist and return early + data = sweep[mtype].get(moment, False) + if data is not False: + return + + data = OrderedDict() + data[moment] = self.msg_31_data_header[sweep_number]["msg_31_data_header"].pop( + moment + ) + sweep.update(self.msg_31_data_header[sweep_number]) + sweep[mtype].update(data) + + def get_metadata_header(self): + """Get metadaata header""" + # data offsets + # ICD 2620010E + # 7.3.5 Metadata Record + meta_offsets = dict( + msg_15=(0, 77), + msg_13=(77, 126), + msg_18=(126, 131), + msg_3=(131, 132), + msg_5=(132, 133), + msg_2=(133, 134), + ) + meta_headers = {} + for msg, (ms, me) in meta_offsets.items(): + mheader = [] + for rec in np.arange(ms, me): + self.init_record(rec) + filepos = self.filepos + message_header = self.get_message_header() + # do not read zero blocks of data + if message_header["size"] == 0: + break + message_header["record_number"] = self.record_number + message_header["filepos"] = filepos + mheader.append(message_header) + meta_headers[msg] = mheader + return meta_headers + + def get_msg_5_data(self): + """Get MSG5 data.""" + self.init_record(132) + self.rh.pos += LEN_MSG_HEADER + 12 + # unpack msg header + msg_5 = _unpack_dictionary( + self._rh.read(LEN_MSG_5, width=1), + MSG_5, + self._rawdata, + byte_order=">", + ) + msg_5["elevation_data"] = [] + # unpack elevation cuts + for i in range(msg_5["number_elevation_cuts"]): + msg_5_elev = _unpack_dictionary( + self._rh.read(LEN_MSG_5_ELEV, width=1), + MSG_5_ELEV, + self._rawdata, + byte_order=">", + ) + msg_5["elevation_data"].append(msg_5_elev) + return msg_5 + + def get_message_header(self): + """Read and unpack message header.""" + self._rh.pos += 12 + return _unpack_dictionary( + self._rh.read(LEN_MSG_HEADER, width=1), + MSG_HEADER, + self._rawdata, + byte_order=">", + ) + + def get_data(self, sweep_number, moment=None): + """Load sweep data from file.""" + sweep = self.data[sweep_number] + start = sweep["record_number"] + stop = sweep["record_end"] + intermediate_records = [ + rec["record_number"] for rec in sweep["intermediate_records"] + ] + filepos = sweep["filepos"] + moments = sweep["sweep_data"] + if moment is None: + moment = moments + elif isinstance(moment, str): + moment = [moment] + for name in moment: + if self.is_compressed: + self.init_record(start) + else: + self.init_record_by_filepos(start, filepos) + ngates = moments[name]["ngates"] + word_size = moments[name]["word_size"] + data_offset = moments[name]["data_offset"] + ws = {8: 1, 16: 2} + width = ws[word_size] + data = [] + self.rh.pos += data_offset + data.append(self._rh.read(ngates, width=width).view(f"uint{word_size}")) + while self.init_next_record() and self.record_number <= stop: + if self.record_number in intermediate_records: + continue + self.rh.pos += data_offset + data.append(self._rh.read(ngates, width=width).view(f"uint{word_size}")) + moments[name].update(data=data) + + def get_data_header(self): + """Load all data header from file.""" + self.init_record(133) + current_sweep = -1 + current_header = -1 + sweep_msg_31_header = [] + sweep_intermediate_records = [] + sweep = OrderedDict() + + data_header = [] + _msg_31_header = [] + _msg_31_data_header = [] + + while self.init_next_record(): + current_header += 1 + # get message headers + msg_header = self.get_message_header() + # append to data headers list + msg_header["record_number"] = self.record_number + msg_header["filepos"] = self.filepos + + # keep all data headers + data_header.append(msg_header) + # only check type 31 for now + if msg_header["type"] == 31: + # get msg_31_header + msg_31_header = _unpack_dictionary( + self._rh.read(LEN_MSG_31, width=1), + MSG_31, + self._rawdata, + byte_order=">", + ) + # add record_number/filepos + msg_31_header["record_number"] = self.record_number + 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 = [] + current_sweep += 1 + sweep["record_number"] = self.rh.recnum + sweep["filepos"] = self.filepos + sweep["msg_31_header"] = msg_31_header + block_pointers = [ + v + for k, v in msg_31_header.items() + if k.startswith("block_pointer") and v > 0 + ] + block_header = {} + for i, block_pointer in enumerate( + block_pointers[: msg_31_header["block_count"]] + ): + self.rh.pos = block_pointer + 12 + LEN_MSG_HEADER + buf = self._rh.read(4, width=1) + dheader = _unpack_dictionary( + buf, + DATA_BLOCK_HEADER, + self._rawdata, + byte_order=">", + ) + block = DATA_BLOCK_TYPE_IDENTIFIER[dheader["block_type"]][ + dheader["data_name"] + ] + LEN_BLOCK = struct.calcsize( + _get_fmt_string(block, byte_order=">") + ) + block_header[dheader["data_name"]] = _unpack_dictionary( + self._rh.read(LEN_BLOCK, width=1), + block, + self._rawdata, + byte_order=">", + ) + if dheader["block_type"] == "D": + block_header[dheader["data_name"]][ + "data_offset" + ] = self.rh.pos + + sweep["msg_31_data_header"] = block_header + _msg_31_data_header.append(sweep) + sweep_msg_31_header.append(msg_31_header) + 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): + """Checks record for correct size. + + Returns + ------- + chk : bool + False, if record is truncated. + """ + chk = self._rh.record.shape[0] == self.record_size + if not chk: + raise EOFError(f"Unexpected file end detected at record {self.rh.recnum}") + return chk + + +# Figure 1 in Interface Control Document for the Archive II/User +# page 7-2 +VOLUME_HEADER = OrderedDict( + [ + ("tape", {"fmt": "9s"}), + ("extension", {"fmt": "3s"}), + ("date", UINT4), + ("time", UINT4), + ("icao", {"fmt": "4s"}), + ] +) + +# Table II Message Header Data +# page 3-7 +MSG_HEADER = OrderedDict( + [ + ("size", UINT2), # size of data, no including header + ("channels", UINT1), + ("type", UINT1), + ("seq_id", UINT2), + ("date", UINT2), + ("ms", UINT4), + ("segments", UINT2), + ("seg_num", UINT2), + ] +) +LEN_MSG_HEADER = struct.calcsize(_get_fmt_string(MSG_HEADER, byte_order=">")) + +# Table XVII Digital Radar Generic Format Blocks (Message Type 31) +# pages 3-87 to 3-89 +MSG_31 = OrderedDict( + [ + ("id", {"fmt": "4s"}), # 0-3 + ("collect_ms", UINT4), # 4-7 + ("collect_date", UINT2), # 8-9 + ("azimuth_number", UINT2), # 10-11 + ("azimuth_angle", FLT4), # 12-15 + ("compress_flag", CODE1), # 16 + ("spare_0", UINT1), # 17 + ("radial_length", UINT2), # 18-19 + ("azimuth_resolution", CODE1), # 20 + ("radial_status", CODE1), # 21 + ("elevation_number", UINT1), # 22 + ("cut_sector", UINT1), # 23 + ("elevation_angle", FLT4), # 24-27 + ("radial_blanking", CODE1), # 28 + ("azimuth_mode", SINT1), # 29 + ("block_count", UINT2), # 30-31 + ("block_pointer_1", UINT4), # 32-35 Volume Data Constant XVII-E + ("block_pointer_2", UINT4), # 36-39 Elevation Data Constant XVII-F + ("block_pointer_3", UINT4), # 40-43 Radial Data Constant XVII-H + ("block_pointer_4", UINT4), # 44-47 Moment "REF" XVII-{B/I} + ("block_pointer_5", UINT4), # 48-51 Moment "VEL" + ("block_pointer_6", UINT4), # 52-55 Moment "SW" + ("block_pointer_7", UINT4), # 56-59 Moment "ZDR" + ("block_pointer_8", UINT4), # 60-63 Moment "PHI" + ("block_pointer_9", UINT4), # 64-67 Moment "RHO" + ("block_pointer_10", UINT4), # Moment "CFP" + ] +) +LEN_MSG_31 = struct.calcsize(_get_fmt_string(MSG_31, byte_order=">")) + + +# Table III Digital Radar Data (Message Type 1) +# pages 3-7 to +MSG_1 = OrderedDict( + [ + ("collect_ms", UINT4), # 0-3 + ("collect_date", UINT2), # 4-5 + ("unambig_range", SINT2), # 6-7 + ("azimuth_angle", CODE2), # 8-9 + ("azimuth_number", UINT2), # 10-11 + ("radial_status", CODE2), # 12-13 + ("elevation_angle", UINT2), # 14-15 + ("elevation_number", UINT2), # 16-17 + ("sur_range_first", CODE2), # 18-19 + ("doppler_range_first", CODE2), # 20-21 + ("sur_range_step", CODE2), # 22-23 + ("doppler_range_step", CODE2), # 24-25 + ("sur_nbins", UINT2), # 26-27 + ("doppler_nbins", UINT2), # 28-29 + ("cut_sector_num", UINT2), # 30-31 + ("calib_const", FLT4), # 32-35 + ("sur_pointer", UINT2), # 36-37 + ("vel_pointer", UINT2), # 38-39 + ("width_pointer", UINT2), # 40-41 + ("doppler_resolution", CODE2), # 42-43 + ("vcp", UINT2), # 44-45 + ("spare_1", {"fmt": "8s"}), # 46-53 + ("spare_2", {"fmt": "2s"}), # 54-55 + ("spare_3", {"fmt": "2s"}), # 56-57 + ("spare_4", {"fmt": "2s"}), # 58-59 + ("nyquist_vel", SINT2), # 60-61 + ("atmos_attenuation", SINT2), # 62-63 + ("threshold", SINT2), # 64-65 + ("spot_blank_status", UINT2), # 66-67 + ("spare_5", {"fmt": "32s"}), # 68-99 + # 100+ reflectivity, velocity and/or spectral width data, CODE1 + ] +) +LEN_MSG_1 = struct.calcsize(_get_fmt_string(MSG_1, byte_order=">")) + +# Table XI Volume Coverage Pattern Data (Message Type 5 & 7) +# pages 3-51 to 3-54 +MSG_5 = OrderedDict( + [ + ("message_size", UINT2), + ("pattern_type", CODE2), + ("pattern_number", UINT2), + ("number_elevation_cuts", UINT2), + ("clutter_map_group_number", UINT2), + ("doppler_velocity_resolution", CODE1), # 2: 0.5 degrees, 4: 1.0 degrees + ("pulse_width", CODE1), # 2: short, 4: long + ("spare", {"fmt": "10s"}), # halfwords 7-11 (10 bytes, 5 halfwords) + ] +) +LEN_MSG_5 = struct.calcsize(_get_fmt_string(MSG_5, byte_order=">")) + +MSG_5_ELEV = OrderedDict( + [ + ("elevation_angle", BIN2), # scaled by 360/65536 for value in degrees. + ("channel_config", CODE1), + ("waveform_type", CODE1), + ("super_resolution", CODE1), + ("prf_number", UINT1), + ("prf_pulse_count", UINT2), + ("azimuth_rate", CODE2), + ("ref_thresh", SINT2), + ("vel_thresh", SINT2), + ("sw_thresh", SINT2), + ("zdr_thres", SINT2), + ("phi_thres", SINT2), + ("rho_thres", SINT2), + ("edge_angle_1", CODE2), + ("dop_prf_num_1", UINT2), + ("dop_prf_pulse_count_1", UINT2), + ("spare_1", {"fmt": "2s"}), + ("edge_angle_2", CODE2), + ("dop_prf_num_2", UINT2), + ("dop_prf_pulse_count_2", UINT2), + ("spare_2", {"fmt": "2s"}), + ("edge_angle_3", CODE2), + ("dop_prf_num_3", UINT2), + ("dop_prf_pulse_count_3", UINT2), + ("spare_3", {"fmt": "2s"}), + ] +) +LEN_MSG_5_ELEV = struct.calcsize(_get_fmt_string(MSG_5_ELEV, byte_order=">")) + +MSG_18 = OrderedDict( + [ + ("adap_file_name", {"fmt": "12s"}), + ("adap_format", {"fmt": "4s"}), + ("adap_revision", {"fmt": "4s"}), + ("adap_date", {"fmt": "12s"}), + ("adap_time", {"fmt": "12s"}), + ("k1", FLT4), + ("az_lat", FLT4), + ("k3", FLT4), + ("el_lat", FLT4), + ("park_az", FLT4), + ("park_el", FLT4), + ("a_fuel_conv0", FLT4), + ("a_fuel_conv1", FLT4), + ("a_fuel_conv2", FLT4), + ("a_fuel_conv3", FLT4), + ("a_fuel_conv4", FLT4), + ("a_fuel_conv5", FLT4), + ("a_fuel_conv6", FLT4), + ("a_fuel_conv7", FLT4), + ("a_fuel_conv8", FLT4), + ("a_fuel_conv9", FLT4), + ("a_fuel_conv10", FLT4), + ("a_min_shelter_temp", FLT4), + ("a_max_shelter_temp", FLT4), + ("a_min_shelter_ac_temp_diff", FLT4), + ("a_max_xmtr_air_temp", FLT4), + ("a_max_rad_temp", FLT4), + ("a_max_rad_temp_rise", FLT4), + ("ped_28v_reg_lim", FLT4), + ("ped_5v_reg_lim", FLT4), + ("ped_15v_reg_lim", FLT4), + ("a_min_gen_room_temp", FLT4), + ("a_max_gen_room_temp", FLT4), + ("dau_5v_reg_lim", FLT4), + ("dau_15v_reg_lim", FLT4), + ("dau_28v_reg_lim", FLT4), + ("en_5v_reg_lim", FLT4), + ("en_5v_nom_volts", FLT4), + ("rpg_co_located", {"fmt": "4s"}), + ("spec_filter_installed", {"fmt": "4s"}), + ("tps_installed", {"fmt": "4s"}), + ("rms_installed", {"fmt": "4s"}), + ("a_hvdl_tst_int", UINT4), + ("a_rpg_lt_int", UINT4), + ("a_min_stab_util_pwr_time", UINT4), + ("a_gen_auto_exer_interval", UINT4), + ("a_util_pwr_sw_req_interval", UINT4), + ("a_low_fuel_level", FLT4), + ("config_chan_number", UINT4), + ("a_rpg_link_type", UINT4), + ("redundant_chan_config", UINT4), + ] +) +for i in np.arange(104): + MSG_18[f"atten_table_{i:03d}"] = FLT4 +MSG_18.update( + [ + ("spare_01", {"fmt": "24s"}), + ("path_losses_07", FLT4), + ("spare_02", {"fmt": "12s"}), + ("path_losses_11", FLT4), + ("spare_03", {"fmt": "4s"}), + ("path_losses_13", FLT4), + ("path_losses_14", FLT4), + ("spare_04", {"fmt": "16s"}), + ("path_losses_19", FLT4), + ("spare_05", {"fmt": "32s"}), + ] +) +for i in np.arange(28, 48): + MSG_18[f"path_losses_{i:02d}"] = FLT4 +MSG_18.update( + [ + ("h_coupler_cw_loss", FLT4), + ("v_coupler_xmt_loss", FLT4), + ("path_losses_50", FLT4), + ("spare_06", {"fmt": "4s"}), + ("path_losses_52", FLT4), + ("v_coupler_cw_loss", FLT4), + ("path_losses_54", FLT4), + ("spare_07", {"fmt": "12s"}), + ("path_losses_58", FLT4), + ("path_losses_59", FLT4), + ("path_losses_60", FLT4), + ("path_losses_61", FLT4), + ("spare_08", {"fmt": "4s"}), + ("path_losses_63", FLT4), + ("path_losses_64", FLT4), + ("path_losses_65", FLT4), + ("path_losses_66", FLT4), + ("path_losses_67", FLT4), + ("path_losses_68", FLT4), + ("path_losses_69", FLT4), + ("chan_cal_diff", FLT4), + ("spare_09", {"fmt": "4s"}), + ("log_amp_factor_1", FLT4), + ("log_amp_factor_2", FLT4), + ("v_ts_cw", FLT4), + ] +) +for i in np.arange(13): + MSG_18[f"rnscale_{i:02d}"] = FLT4 +for i in np.arange(13): + MSG_18[f"atmos_{i:02d}"] = FLT4 +for i in np.arange(12): + MSG_18[f"el_index_{i:02d}"] = FLT4 +MSG_18.update( + [ + ("tfreq_mhz", UINT4), + ("base_data_tcn", FLT4), + ("refl_data_tover", FLT4), + ("tar_h_dbz0_lp", FLT4), + ("tar_v_dbz0_lp", FLT4), + ("init_phi_dp", UINT4), + ("norm_init_phi_dp", UINT4), + ("lx_lp", FLT4), + ("lx_sp", FLT4), + ("meteor_param", FLT4), + ("beamwidth", FLT4), + ("antenna_gain", FLT4), + ("spare_10", {"fmt": "4s"}), + ("vel_maint_limit", FLT4), + ("wth_maint_limit", FLT4), + ("vel_degrad_limit", FLT4), + ("wth_degrad_limit", FLT4), + ("h_noisetemp_degrad_limit", FLT4), + ("h_noisetemp_maint_limit", FLT4), + ("v_noisetemp_degrad_limit", FLT4), + ("v_noisetemp_maint_limit", FLT4), + ("kly_degrade_limit", FLT4), + ("ts_coho", FLT4), + ("h_ts_cw", FLT4), + ("ts_rf_sp", FLT4), + ("ts_rf_lp", FLT4), + ("ts_stalo", FLT4), + ("ame_h_noise_enr", FLT4), + ("xmtr_peak_pwr_high_limit", FLT4), + ("xmtr_peak_pwr_low_limit", FLT4), + ("h_dbz0_delta_limit", FLT4), + ("threshold1", FLT4), + ("threshold2", FLT4), + ("clut_supp_dgrad_lim", FLT4), + ("clut_supp_maint_lim", FLT4), + ("range0_value", FLT4), + ("xmtr_pwr_mtr_scale", FLT4), + ("v_dbz0_delta_limit", FLT4), + ("tar_h_dbz0_sp", FLT4), + ("tar_v_dbz0_sp", FLT4), + ("deltaprf", UINT4), + ("spare_11", {"fmt": "8s"}), + ("tau_sp", UINT4), + ("tau_lp", UINT4), + ("nc_dead_value", UINT4), + ("tau_rf_sp", UINT4), + ("tau_rf_lp", UINT4), + ("seg1lim", FLT4), + ("slatsec", FLT4), + ("slonsec", FLT4), + ("spare_12", {"fmt": "4s"}), + ("slatdeg", UINT4), + ("slatmin", UINT4), + ("slondeg", UINT4), + ("slonmin", UINT4), + ("slatdir", {"fmt": "4s"}), + ("slondir", {"fmt": "4s"}), + ("spare_13", {"fmt": "4s"}), + ("vcpat_11", {"fmt": "1172s"}), + ("vcpat_21", {"fmt": "1172s"}), + ("vcpat_31", {"fmt": "1172s"}), + ("vcpat_32", {"fmt": "1172s"}), + ("vcpat_300", {"fmt": "1172s"}), + ("vcpat_301", {"fmt": "1172s"}), + ("az_correction_factor", FLT4), + ("el_correction_factor", FLT4), + ("site_name", {"fmt": "4s"}), + ("ant_manual_setup_ielmin", SINT4), + ("ant_manual_setup_ielmax", SINT4), + ("ant_manual_setup_fazvelmax", SINT4), + ("ant_manual_setup_felvelmax", SINT4), + ("ant_manual_setup_ignd_hgt", SINT4), + ("ant_manual_setup_irad_hgt", SINT4), + ("spare_14", {"fmt": "300s"}), + ("rvp8nv_iwaveguide_lenght", UINT4), + ("spare_15", {"fmt": "44s"}), + ("vel_data_tover", FLT4), + ("width_data_tover", FLT4), + ("spare_16", {"fmt": "12s"}), + ("doppler_range_start", FLT4), + ("max_el_index", UINT4), + ("seg2lim", FLT4), + ("seg3lim", FLT4), + ("seg4lim", FLT4), + ("nbr_el_segments", UINT4), + ("h_noise_long", FLT4), + ("ant_noise_temp", FLT4), + ("h_noise_short", FLT4), + ("h_noise_tolerance", FLT4), + ("min_h_dyn_range", FLT4), + ("gen_installed", {"fmt": "4s"}), + ("gen_exercise", {"fmt": "4s"}), + ("v_noise_tolerance", FLT4), + ("min_v_dyn_range", FLT4), + ("zdr_bias_dgrad_lim", FLT4), + ("spare_17", {"fmt": "16s"}), + ("v_noise_long", FLT4), + ("v_noise_short", FLT4), + ("zdr_data_tover", FLT4), + ("phi_data_tover", FLT4), + ("rho_data_tover", FLT4), + ("stalo_power_dgrad_limit", FLT4), + ("stalo_power_maint_limit", FLT4), + ("min_h_pwr_sense", FLT4), + ("min_v_pwr_sense", FLT4), + ("h_pwr_sense_offset", FLT4), + ("v_pwr_sense_offset", FLT4), + ("ps_gain_ref", FLT4), + ("rf_pallet_broad_loss", FLT4), + ("zdr_check_threshold", FLT4), + ("phi_check_threshold", FLT4), + ("rho_check_threshold", FLT4), + ("spare_18", {"fmt": "52s"}), + ("ame_ps_tolerance", FLT4), + ("ame_max_temp", FLT4), + ("ame_min_temp", FLT4), + ("rcvr_mod_max_temp", FLT4), + ("rcvr_mod_min_temp", FLT4), + ("bite_mod_max_temp", FLT4), + ("bite_mod_min_temp", FLT4), + ("default_polarization", UINT4), + ("tr_limit_dgrad_limit", FLT4), + ("tr_limit_fail_limit", FLT4), + ("spare_19", {"fmt": "8s"}), + ("ame_current_tolerance", FLT4), + ("h_only_polarization", UINT4), + ("v_only_polarization", UINT4), + ("spare_20", {"fmt": "8s"}), + ("reflector_bias", FLT4), + ("a_min_shelter_temp_warn", FLT4), + ("spare_21", {"fmt": "432s"}), + ] +) +LEN_MSG_18 = struct.calcsize(_get_fmt_string(MSG_18, byte_order=">")) + +DATA_BLOCK_HEADER = OrderedDict( + [("block_type", string_dict(1)), ("data_name", string_dict(3))] +) +LEN_DATA_BLOCK_HEADER = struct.calcsize( + _get_fmt_string(DATA_BLOCK_HEADER, byte_order=">") +) + + +# Table XVII-B Data Block (Descriptor of Generic Data Moment Type) +# pages 3-90 and 3-91 +GENERIC_DATA_BLOCK = OrderedDict( + [ + ("reserved", UINT4), + ("ngates", UINT2), + ("first_gate", SINT2), + ("gate_spacing", SINT2), + ("thresh", SINT2), + ("snr_thres", SINT2), + ("flags", CODE1), + ("word_size", UINT1), + ("scale", FLT4), + ("offset", FLT4), + # then data + ] +) +LEN_GENERIC_DATA_BLOCK = struct.calcsize( + _get_fmt_string(GENERIC_DATA_BLOCK, byte_order=">") +) + +# Table XVII-E Data Block (Volume Data Constant Type) +# page 3-92 +VOLUME_DATA_BLOCK = OrderedDict( + [ + ("lrtup", UINT2), + ("version_major", UINT1), + ("version_minor", UINT1), + ("lat", FLT4), + ("lon", FLT4), + ("height", SINT2), + ("feedhorn_height", UINT2), + ("refl_calib", FLT4), + ("power_h", FLT4), + ("power_v", FLT4), + ("diff_refl_calib", FLT4), + ("init_phase", FLT4), + ("vcp", UINT2), + ("spare", {"fmt": "2s"}), + ] +) +LEN_VOLUME_DATA_BLOCK = struct.calcsize( + _get_fmt_string(VOLUME_DATA_BLOCK, byte_order=">") +) + +# Table XVII-F Data Block (Elevation Data Constant Type) +# page 3-93 +ELEVATION_DATA_BLOCK = OrderedDict( + [ + ("lrtup", UINT2), + ("atmos", SINT2), + ("refl_calib", FLT4), + ] +) +LEN_ELEVATION_DATA_BLOCK = struct.calcsize( + _get_fmt_string(ELEVATION_DATA_BLOCK, byte_order=">") +) + +# Table XVII-H Data Block (Radial Data Constant Type) +# pages 3-93 +RADIAL_DATA_BLOCK = OrderedDict( + [ + ("lrtup", UINT2), + ("unambig_range", SINT2), + ("noise_h", FLT4), + ("noise_v", FLT4), + ("nyquist_vel", SINT2), + ("spare", {"fmt": "2s"}), + ] +) +LEN_RADIAL_DATA_BLOCK = struct.calcsize( + _get_fmt_string(RADIAL_DATA_BLOCK, byte_order=">") +) + +DATA_BLOCK_CONSTANT_IDENTIFIERS = OrderedDict( + [ + ("VOL", VOLUME_DATA_BLOCK), + ("ELV", ELEVATION_DATA_BLOCK), + ("RAD", RADIAL_DATA_BLOCK), + ] +) + +DATA_BLOCK_VARIABLE_IDENTIFIERS = OrderedDict( + [ + ("REF", GENERIC_DATA_BLOCK), + ("VEL", GENERIC_DATA_BLOCK), + ("SW ", GENERIC_DATA_BLOCK), + ("ZDR", GENERIC_DATA_BLOCK), + ("PHI", GENERIC_DATA_BLOCK), + ("RHO", GENERIC_DATA_BLOCK), + ("CFP", GENERIC_DATA_BLOCK), + ] +) + +DATA_BLOCK_TYPE_IDENTIFIER = OrderedDict( + [ + ("R", DATA_BLOCK_CONSTANT_IDENTIFIERS), + ("D", DATA_BLOCK_VARIABLE_IDENTIFIERS), + ] +) + + +class NexradLevel2ArrayWrapper(BackendArray): + """Wraps array of NexradLevel2 data.""" + + def __init__(self, datastore, name, var): + self.datastore = datastore + self.group = datastore._group + self.name = name + # get rays and bins + nrays = ( + datastore.ds["record_end"] + - datastore.ds["record_number"] + + 1 + - len(datastore.ds["intermediate_records"]) + ) + nbins = max([v["ngates"] for k, v in datastore.ds["sweep_data"].items()]) + self.dtype = np.dtype("uint8") + if name == "PHI": + self.dtype = np.dtype("uint16") + self.shape = (nrays, nbins) + + def _getitem(self, key): + # read the data if not available + try: + data = self.datastore.ds["sweep_data"][self.name]["data"] + except KeyError: + self.datastore.root.get_data(self.group, self.name) + data = self.datastore.ds["sweep_data"][self.name]["data"] + if self.name == "PHI": + x = np.uint16(0x3FF) + else: + x = np.uint8(0xFF) + if len(data[0]) < self.shape[1]: + return np.pad( + np.vstack(data) & x, + ((0, 0), (0, self.shape[1] - len(data[0]))), + mode="constant", + constant_values=0, + )[key] + else: + return (np.vstack(data) & x)[key] + + def __getitem__(self, key): + return indexing.explicit_indexing_adapter( + key, + self.shape, + indexing.IndexingSupport.BASIC, + self._getitem, + ) + + +class NexradLevel2Store(AbstractDataStore): + def __init__(self, manager, group=None): + self._manager = manager + self._group = int(group[6:]) + self._filename = self.filename + + @classmethod + def open(cls, filename, mode="r", group=None, **kwargs): + manager = CachingFileManager( + NEXRADLevel2File, filename, mode=mode, kwargs=kwargs + ) + return cls(manager, group=group) + + @property + def filename(self): + with self._manager.acquire_context(False) as root: + return root.filename + + @property + def root(self): + with self._manager.acquire_context(False) as root: + return root + + def _acquire(self, needs_lock=True): + with self._manager.acquire_context(needs_lock) as root: + root.get_sweep(self._group) + ds = root.data[self._group] + return ds + + @property + def ds(self): + return self._acquire() + + def open_store_variable(self, name, var): + dim = "azimuth" + + data = indexing.LazilyOuterIndexedArray( + NexradLevel2ArrayWrapper(self, name, var) + ) + encoding = {"group": self._group, "source": self._filename} + + mname = nexrad_mapping.get(name, name) + mapping = sweep_vars_mapping.get(mname, {}) + attrs = {key: mapping[key] for key in moment_attrs if key in mapping} + attrs["scale_factor"] = 1.0 / var["scale"] + attrs["add_offset"] = -var["offset"] / var["scale"] + attrs["coordinates"] = ( + "elevation azimuth range latitude longitude altitude time" + ) + return mname, Variable((dim, "range"), data, attrs, encoding) + + def open_store_coordinates(self): + msg_31_header = self.root.msg_31_header[self._group] + + # azimuth/elevation + azimuth = np.array([ms["azimuth_angle"] for ms in msg_31_header]) + 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 + milliseconds = np.array([ms["collect_ms"] for ms in msg_31_header]) + rtime = date + milliseconds + time_prefix = "milli" + rtime_attrs = get_time_attrs(date_unit=f"{time_prefix}seconds") + + # site coords + vol = self.ds["sweep_constant_data"]["VOL"] + lon, lat, alt = vol["lon"], vol["lat"], vol["height"] + vol["feedhorn_height"] + + # range + sweep_data = list(self.ds["sweep_data"].values())[0] + first_gate = sweep_data["first_gate"] + gate_spacing = sweep_data["gate_spacing"] + ngates = max([v["ngates"] for k, v in self.ds["sweep_data"].items()]) + last_gate = first_gate + gate_spacing * (ngates - 0.5) + rng = np.arange(first_gate, last_gate, gate_spacing, "float32") + range_attrs = get_range_attrs(rng) + + encoding = {"group": self._group} + dim = "azimuth" + sweep_mode = "azimuth_surveillance" + sweep_number = self._group + prt_mode = "not_set" + follow_mode = "not_set" + + fixed_angle = self.root.msg_5["elevation_data"][self._group] + + coords = { + "azimuth": Variable((dim,), azimuth, get_azimuth_attrs(), encoding), + "elevation": Variable((dim,), elevation, get_elevation_attrs(), encoding), + "time": Variable((dim,), rtime, rtime_attrs, encoding), + "range": Variable(("range",), rng, range_attrs), + "longitude": Variable((), lon, get_longitude_attrs()), + "latitude": Variable((), lat, get_latitude_attrs()), + "altitude": Variable((), alt, get_altitude_attrs()), + "sweep_mode": Variable((), sweep_mode), + "sweep_number": Variable((), sweep_number), + "prt_mode": Variable((), prt_mode), + "follow_mode": Variable((), follow_mode), + "sweep_fixed_angle": Variable((), fixed_angle), + } + + return coords + + def get_variables(self): + return FrozenDict( + (k1, v1) + for k1, v1 in { + **dict( + self.open_store_variable(k, v) + for k, v in self.ds["sweep_data"].items() + ), + **self.open_store_coordinates(), + }.items() + ) + + def get_attrs(self): + attributes = [("instrument_name", self.root.volume_header["icao"].decode())] + + return FrozenDict(attributes) + + +class NexradLevel2BackendEntrypoint(BackendEntrypoint): + """ + Xarray BackendEntrypoint for NEXRAD Level2 Data + """ + + description = "Open NEXRAD Level2 files in Xarray" + url = "tbd" + + def open_dataset( + self, + filename_or_obj, + *, + mask_and_scale=True, + decode_times=True, + concat_characters=True, + decode_coords=True, + drop_variables=None, + use_cftime=None, + decode_timedelta=None, + group=None, + first_dim="auto", + reindex_angle=False, + fix_second_angle=False, + site_coords=True, + optional=True, + ): + store = NexradLevel2Store.open( + filename_or_obj, + group=group, + loaddata=False, + ) + + store_entrypoint = StoreBackendEntrypoint() + + ds = store_entrypoint.open_dataset( + store, + mask_and_scale=mask_and_scale, + decode_times=decode_times, + concat_characters=concat_characters, + decode_coords=decode_coords, + drop_variables=drop_variables, + use_cftime=use_cftime, + decode_timedelta=decode_timedelta, + ) + + # reassign azimuth/elevation/time coordinates + ds = ds.assign_coords({"azimuth": ds.azimuth}) + ds = ds.assign_coords({"elevation": ds.elevation}) + ds = ds.assign_coords({"time": ds.time}) + + ds.encoding["engine"] = "nexradlevel2" + + # handle duplicates and reindex + if decode_coords and reindex_angle is not False: + ds = ds.pipe(util.remove_duplicate_rays) + ds = ds.pipe(util.reindex_angle, **reindex_angle) + ds = ds.pipe(util.ipol_time, **reindex_angle) + + # handling first dimension + dim0 = "elevation" if ds.sweep_mode.load() == "rhi" else "azimuth" + + # todo: could be optimized + if first_dim == "time": + ds = ds.swap_dims({dim0: "time"}) + ds = ds.sortby("time") + else: + ds = ds.sortby(dim0) + + # assign geo-coords + if site_coords: + ds = ds.assign_coords( + { + "latitude": ds.latitude, + "longitude": ds.longitude, + "altitude": ds.altitude, + } + ) + + return ds + + +def open_nexradlevel2_datatree(filename_or_obj, **kwargs): + """Open NEXRAD Level2 dataset as :py:class:`datatree.DataTree`. + + Parameters + ---------- + filename_or_obj : str, Path, file-like or DataStore + Strings and Path objects are interpreted as a path to a local or remote + radar file + + Keyword Arguments + ----------------- + sweep : int, list of int, optional + Sweep number(s) to extract, default to first sweep. If None, all sweeps are + extracted into a list. + first_dim : str + Can be ``time`` or ``auto`` first dimension. If set to ``auto``, + first dimension will be either ``azimuth`` or ``elevation`` depending on + type of sweep. Defaults to ``auto``. + reindex_angle : bool or dict + Defaults to False, no reindexing. Given dict should contain the kwargs to + reindex_angle. Only invoked if `decode_coord=True`. + fix_second_angle : bool + If True, fixes erroneous second angle data. Defaults to ``False``. + site_coords : bool + Attach radar site-coordinates to Dataset, defaults to ``True``. + kwargs : dict + Additional kwargs are fed to :py:func:`xarray.open_dataset`. + + Returns + ------- + dtree: datatree.DataTree + DataTree + """ + # handle kwargs, extract first_dim + backend_kwargs = kwargs.pop("backend_kwargs", {}) + # first_dim = backend_kwargs.pop("first_dim", None) + sweep = kwargs.pop("sweep", None) + sweeps = [] + kwargs["backend_kwargs"] = backend_kwargs + + if isinstance(sweep, str): + sweeps = [sweep] + elif isinstance(sweep, int): + sweeps = [f"sweep_{sweep}"] + elif isinstance(sweep, list): + if isinstance(sweep[0], int): + sweeps = [f"sweep_{i}" for i in sweep] + else: + sweeps.extend(sweep) + else: + with NEXRADLevel2File(filename_or_obj, loaddata=False) as nex: + nsweeps = nex.msg_5["number_elevation_cuts"] + sweeps = [f"sweep_{i}" for i in range(nsweeps)] + + engine = NexradLevel2BackendEntrypoint + + ds = [ + xr.open_dataset(filename_or_obj, group=swp, engine=engine, **kwargs) + for swp in sweeps + ] + + ds.insert(0, xr.Dataset()) + + # create datatree root node with required data + dtree = DataTree(data=_assign_root(ds), name="root") + # return datatree with attached sweep child nodes + return _attach_sweep_groups(dtree, ds[1:]) diff --git a/xradar/util.py b/xradar/util.py index c4441a17..83666528 100644 --- a/xradar/util.py +++ b/xradar/util.py @@ -21,6 +21,7 @@ "reindex_angle", "extract_angle_parameters", "ipol_time", + "rolling_dim", ] __doc__ = __doc__.format("\n ".join(__all__)) @@ -479,3 +480,10 @@ def _get_data_file(file, file_or_filelike): yield io.BytesIO(f.read()) else: yield file + + +def rolling_dim(data, window): + """Return array with rolling dimension of window-length added at the end.""" + shape = data.shape[:-1] + (data.shape[-1] - window + 1, window) + strides = data.strides + (data.strides[-1],) + return np.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)