diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d1a621e4..50086b0a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,9 +22,6 @@ jobs: run: | python -m pip install --upgrade pip pip install black black[jupyter] ruff - - name: Black style check - run: | - black --check . - name: Lint with ruff run: | ruff . diff --git a/.github/workflows/upload_pypi.yml b/.github/workflows/upload_pypi.yml index 6c84dd09..767dd1a6 100644 --- a/.github/workflows/upload_pypi.yml +++ b/.github/workflows/upload_pypi.yml @@ -30,4 +30,5 @@ jobs: TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }} run: | python -m build + python setup.py build_ext --inplace twine upload dist/* diff --git a/MANIFEST.in b/MANIFEST.in index 109f8eaa..030399f8 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -8,4 +8,5 @@ recursive-include tests * recursive-exclude * __pycache__ recursive-exclude * *.py[co] +global-include *.pyx *pxd recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif diff --git a/docs/history.md b/docs/history.md index 4d37690f..3ccac141 100644 --- a/docs/history.md +++ b/docs/history.md @@ -2,6 +2,7 @@ ## Development Version (unreleased) +* ENH: Add a reader for nexrad level2 files ({pull}`147`) by [@mgrover1](https://github.com/mgrover1) * MNT: Update GitHub actions, address DeprecationWarnings ({pull}`153`) by [@kmuehlbauer](https://github.com/kmuehlbauer). * 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). diff --git a/pyproject.toml b/pyproject.toml index 1c91165a..1ece1c83 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,12 +39,14 @@ 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 = [ "setuptools>=45", "wheel", "setuptools_scm[toml]>=7.0", + "numpy" ] build-backend = "setuptools.build_meta" diff --git a/tests/conftest.py b/tests/conftest.py index b1a0ca11..e17a954e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -58,3 +58,8 @@ 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") diff --git a/tests/io/backends/test_common.py b/tests/io/backends/test_common.py new file mode 100644 index 00000000..0a0afb6b --- /dev/null +++ b/tests/io/backends/test_common.py @@ -0,0 +1,9 @@ +from xradar.io.backends import common + + +def test_lazy_dict(): + d = common.LazyLoadDict({"key1": "value1", "key2": "value2"}) + assert d["key1"] == "value1" + lazy_func = lambda: 999 + d.set_lazy("lazykey1", lazy_func) + assert d["lazykey1"] == 999 diff --git a/tests/io/test_nexrad_level2.py b/tests/io/test_nexrad_level2.py new file mode 100644 index 00000000..2a476e52 --- /dev/null +++ b/tests/io/test_nexrad_level2.py @@ -0,0 +1,30 @@ +#!/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.""" + +import xarray as xr + +from xradar.io.backends import open_nexradlevel2_datatree + + +def test_open_nexradlevel2_datatree(nexradlevel2_file): + dtree = open_nexradlevel2_datatree(nexradlevel2_file) + ds = dtree["sweep_0"] + assert ds.attrs["instrument_name"] == "KATX" + assert ds.attrs["nsweeps"] == 16 + assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters" + assert ds["DBZH"].shape == (719, 1832) + assert ds["DBZH"].dims == ("azimuth", "range") + assert int(ds.sweep_number.values) == 0 + + +def test_open_nexrad_level2_backend(nexradlevel2_file): + ds = xr.open_dataset(nexradlevel2_file, engine="nexradlevel2") + assert ds.attrs["instrument_name"] == "KATX" + assert ds.attrs["nsweeps"] == 16 + assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters" + assert ds["DBZH"].shape == (719, 1832) + assert ds["DBZH"].dims == ("azimuth", "range") + assert int(ds.sweep_number.values) == 0 diff --git a/xradar/io/__init__.py b/xradar/io/__init__.py index 560175ef..061cad19 100644 --- a/xradar/io/__init__.py +++ b/xradar/io/__init__.py @@ -13,7 +13,24 @@ .. automodule:: xradar.io.export """ -from .backends import * # noqa +from .backends import ( + CfRadial1BackendEntrypoint, # noqa + FurunoBackendEntrypoint, # noqa + GamicBackendEntrypoint, # noqa + IrisBackendEntrypoint, # noqa + NexradLevel2BackendEntrypoint, # noqa + OdimBackendEntrypoint, # noqa + RainbowBackendEntrypoint, # noqa + open_cfradial1_datatree, # noqa + open_furuno_datatree, # noqa + open_gamic_datatree, # noqa + open_iris_datatree, # noqa + open_nexradlevel2_datatree, # noqa + open_odim_datatree, # noqa + open_rainbow_datatree, # noqa +) + +# noqa from .export import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 4edf3d6d..81b02784 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,9 @@ from .iris import * # noqa from .odim import * # noqa from .rainbow import * # noqa +from .nexrad_level2 import ( + NexradLevel2BackendEntrypoint, # noqa + open_nexradlevel2_datatree, # 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..f094b11c 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -12,10 +12,15 @@ """ +import bz2 +import gzip import io +import itertools import struct from collections import OrderedDict +from collections.abc import MutableMapping +import fsspec import h5netcdf import numpy as np import xarray as xr @@ -229,3 +234,167 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): UINT1 = {"fmt": "B", "dtype": "unit8"} UINT2 = {"fmt": "H", "dtype": "uint16"} UINT4 = {"fmt": "I", "dtype": "unint32"} + + +def prepare_for_read(filename, storage_options={"anon": True}): + """ + Return a file like object read for reading. + + Open a file for reading in binary mode with transparent decompression of + Gzip and BZip2 files. The resulting file-like object should be closed. + + Parameters + ---------- + filename : str or file-like object + Filename or file-like object which will be opened. File-like objects + will not be examined for compressed data. + + storage_options : dict, optional + Parameters passed to the backend file-system such as Google Cloud Storage, + Amazon Web Service S3. + + Returns + ------- + file_like : file-like object + File like object from which data can be read. + + """ + # if a file-like object was provided, return + if hasattr(filename, "read"): # file-like object + return filename + + # look for compressed data by examining the first few bytes + fh = fsspec.open(filename, mode="rb", compression="infer", **storage_options).open() + magic = fh.read(3) + fh.close() + + # If the data is still compressed, use gunzip/bz2 to uncompress the data + if magic.startswith(b"\x1f\x8b"): + return gzip.GzipFile(filename, "rb") + + if magic.startswith(b"BZh"): + return bz2.BZ2File(filename, "rb") + + return fsspec.open( + filename, mode="rb", compression="infer", **storage_options + ).open() + + +def make_time_unit_str(dtobj): + """Return a time unit string from a datetime object.""" + return "seconds since " + dtobj.strftime("%Y-%m-%dT%H:%M:%SZ") + + +class LazyLoadDict(MutableMapping): + """ + A dictionary-like class supporting lazy loading of specified keys. + + Keys which are lazy loaded are specified using the set_lazy method. + The callable object which produces the specified key is provided as the + second argument to this method. This object gets called when the value + of the key is loaded. After this initial call the results is cached + in the traditional dictionary which is used for supplemental access to + this key. + + Testing for keys in this dictionary using the "key in d" syntax will + result in the loading of a lazy key, use "key in d.keys()" to prevent + this evaluation. + + The comparison methods, __cmp__, __ge__, __gt__, __le__, __lt__, __ne__, + nor the view methods, viewitems, viewkeys, viewvalues, are implemented. + Neither is the fromkeys method. + + This is from Py-ART. + + Parameters + ---------- + dic : dict + Dictionary containing key, value pairs which will be stored and + evaluated traditionally. This dictionary referenced not copied into + the LazyLoadDictionary and hence changed to this dictionary may change + the original. If this behavior is not desired copy dic in the + initalization. + + Examples + -------- + >>> d = LazyLoadDict({'key1': 'value1', 'key2': 'value2'}) + >>> d.keys() + ['key2', 'key1'] + >>> lazy_func = lambda : 999 + >>> d.set_lazy('lazykey1', lazy_func) + >>> d.keys() + ['key2', 'key1', 'lazykey1'] + >>> d['lazykey1'] + 999 + + """ + + def __init__(self, dic): + """initalize.""" + self._dic = dic + self._lazyload = {} + + # abstract methods + def __setitem__(self, key, value): + """Set a key which will not be stored and evaluated traditionally.""" + self._dic[key] = value + if key in self._lazyload: + del self._lazyload[key] + + def __getitem__(self, key): + """Get the value of a key, evaluating a lazy key if needed.""" + if key in self._lazyload: + value = self._lazyload[key]() + self._dic[key] = value + del self._lazyload[key] + return self._dic[key] + + def __delitem__(self, key): + """Remove a lazy or traditional key from the dictionary.""" + if key in self._lazyload: + del self._lazyload[key] + else: + del self._dic[key] + + def __iter__(self): + """Iterate over all lazy and traditional keys.""" + return itertools.chain(self._dic.copy(), self._lazyload.copy()) + + def __len__(self): + """Return the number of traditional and lazy keys.""" + return len(self._dic) + len(self._lazyload) + + # additional class to mimic dict behavior + def __str__(self): + """Return a string representation of the object.""" + if len(self._dic) == 0 or len(self._lazyload) == 0: + seperator = "" + else: + seperator = ", " + lazy_reprs = [(repr(k), repr(v)) for k, v in self._lazyload.items()] + lazy_strs = ["{}: LazyLoad({})".format(*r) for r in lazy_reprs] + lazy_str = ", ".join(lazy_strs) + "}" + return str(self._dic)[:-1] + seperator + lazy_str + + def has_key(self, key): + """True if dictionary has key, else False.""" + return key in self + + def copy(self): + """ + Return a copy of the dictionary. + + Lazy keys are not evaluated in the original or copied dictionary. + """ + dic = self.__class__(self._dic.copy()) + # load all lazy keys into the copy + for key, value_callable in self._lazyload.items(): + dic.set_lazy(key, value_callable) + return dic + + # lazy dictionary specific methods + def set_lazy(self, key, value_callable): + """Set a lazy key to load from a callable object.""" + if key in self._dic: + del self._dic[key] + self._lazyload[key] = value_callable diff --git a/xradar/io/backends/nexrad_common.py b/xradar/io/backends/nexrad_common.py new file mode 100644 index 00000000..13f3c6ac --- /dev/null +++ b/xradar/io/backends/nexrad_common.py @@ -0,0 +1,248 @@ +""" +Data and functions common to all types of NEXRAD files. + +""" + +# The functions in this module are intended to be used in other +# nexrad related modules. The functions are not and should not be +# exported into the pyart.io namespace. + + +def get_nexrad_location(station): + """ + Return the latitude, longitude and altitude of a NEXRAD station. + + Parameters + ---------- + station : str + Four letter NEXRAD station ICAO name. + + Returns + ------- + lat, lon, alt : float + Latitude (in degrees), longitude (in degrees), and altitude + (in meters above mean sea level) of the NEXRAD station. + + """ + loc = NEXRAD_LOCATIONS[station.upper()] + + # Convert from feet to meters for elevation units + loc["elev"] = loc["elev"] * 0.3048 + + return loc["lat"], loc["lon"], loc["elev"] + + +# Locations of NEXRAD locations was retrieved from NOAA's +# Historical Observing Metadata Repository (HOMR) on +# 2014-Mar-27. http://www.ncdc.noaa.gov/homr/ +# The data below was extracted with: +# cut -c 10-14,107-115,117-127,128-133 + +NEXRAD_LOCATIONS = { + "KABR": {"lat": 45.45583, "lon": -98.41306, "elev": 1302}, + "KABX": {"lat": 35.14972, "lon": -106.82333, "elev": 5870}, + "KAKQ": {"lat": 36.98389, "lon": -77.0075, "elev": 112}, + "KAMA": {"lat": 35.23333, "lon": -101.70889, "elev": 3587}, + "KAMX": {"lat": 25.61056, "lon": -80.41306, "elev": 14}, + "KAPX": {"lat": 44.90722, "lon": -84.71972, "elev": 1464}, + "KARX": {"lat": 43.82278, "lon": -91.19111, "elev": 1276}, + "KATX": {"lat": 48.19472, "lon": -122.49444, "elev": 494}, + "KBBX": {"lat": 39.49611, "lon": -121.63167, "elev": 173}, + "KBGM": {"lat": 42.19972, "lon": -75.985, "elev": 1606}, + "KBHX": {"lat": 40.49833, "lon": -124.29194, "elev": 2402}, + "KBIS": {"lat": 46.77083, "lon": -100.76028, "elev": 1658}, + "KBLX": {"lat": 45.85389, "lon": -108.60611, "elev": 3598}, + "KBMX": {"lat": 33.17194, "lon": -86.76972, "elev": 645}, + "KBOX": {"lat": 41.95583, "lon": -71.1375, "elev": 118}, + "KBRO": {"lat": 25.91556, "lon": -97.41861, "elev": 23}, + "KBUF": {"lat": 42.94861, "lon": -78.73694, "elev": 693}, + "KBYX": {"lat": 24.59694, "lon": -81.70333, "elev": 8}, + "KCAE": {"lat": 33.94861, "lon": -81.11861, "elev": 231}, + "KCBW": {"lat": 46.03917, "lon": -67.80694, "elev": 746}, + "KCBX": {"lat": 43.49083, "lon": -116.23444, "elev": 3061}, + "KCCX": {"lat": 40.92306, "lon": -78.00389, "elev": 2405}, + "KCLE": {"lat": 41.41306, "lon": -81.86, "elev": 763}, + "KCLX": {"lat": 32.65556, "lon": -81.04222, "elev": 97}, + "KCRI": {"lat": 35.2383, "lon": -97.4602, "elev": 1201}, + "KCRP": {"lat": 27.78389, "lon": -97.51083, "elev": 45}, + "KCXX": {"lat": 44.51111, "lon": -73.16639, "elev": 317}, + "KCYS": {"lat": 41.15194, "lon": -104.80611, "elev": 6128}, + "KDAX": {"lat": 38.50111, "lon": -121.67667, "elev": 30}, + "KDDC": {"lat": 37.76083, "lon": -99.96833, "elev": 2590}, + "KDFX": {"lat": 29.2725, "lon": -100.28028, "elev": 1131}, + "KDGX": {"lat": 32.28, "lon": -89.98444, "elev": -99999}, + "KDIX": {"lat": 39.94694, "lon": -74.41111, "elev": 149}, + "KDLH": {"lat": 46.83694, "lon": -92.20972, "elev": 1428}, + "KDMX": {"lat": 41.73111, "lon": -93.72278, "elev": 981}, + "KDOX": {"lat": 38.82556, "lon": -75.44, "elev": 50}, + "KDTX": {"lat": 42.69972, "lon": -83.47167, "elev": 1072}, + "KDVN": {"lat": 41.61167, "lon": -90.58083, "elev": 754}, + "KDYX": {"lat": 32.53833, "lon": -99.25417, "elev": 1517}, + "KEAX": {"lat": 38.81028, "lon": -94.26417, "elev": 995}, + "KEMX": {"lat": 31.89361, "lon": -110.63028, "elev": 5202}, + "KENX": {"lat": 42.58639, "lon": -74.06444, "elev": 1826}, + "KEOX": {"lat": 31.46028, "lon": -85.45944, "elev": 434}, + "KEPZ": {"lat": 31.87306, "lon": -106.6975, "elev": 4104}, + "KESX": {"lat": 35.70111, "lon": -114.89139, "elev": 4867}, + "KEVX": {"lat": 30.56417, "lon": -85.92139, "elev": 140}, + "KEWX": {"lat": 29.70361, "lon": -98.02806, "elev": 633}, + "KEYX": {"lat": 35.09778, "lon": -117.56, "elev": 2757}, + "KFCX": {"lat": 37.02417, "lon": -80.27417, "elev": 2868}, + "KFDR": {"lat": 34.36222, "lon": -98.97611, "elev": 1267}, + "KFDX": {"lat": 34.63528, "lon": -103.62944, "elev": 4650}, + "KFFC": {"lat": 33.36333, "lon": -84.56583, "elev": 858}, + "KFSD": {"lat": 43.58778, "lon": -96.72889, "elev": 1430}, + "KFSX": {"lat": 34.57444, "lon": -111.19833, "elev": -99999}, + "KFTG": {"lat": 39.78667, "lon": -104.54528, "elev": 5497}, + "KFWS": {"lat": 32.57278, "lon": -97.30278, "elev": 683}, + "KGGW": {"lat": 48.20639, "lon": -106.62417, "elev": 2276}, + "KGJX": {"lat": 39.06222, "lon": -108.21306, "elev": 9992}, + "KGLD": {"lat": 39.36694, "lon": -101.7, "elev": 3651}, + "KGRB": {"lat": 44.49833, "lon": -88.11111, "elev": 682}, + "KGRK": {"lat": 30.72167, "lon": -97.38278, "elev": 538}, + "KGRR": {"lat": 42.89389, "lon": -85.54472, "elev": 778}, + "KGSP": {"lat": 34.88306, "lon": -82.22028, "elev": 940}, + "KGWX": {"lat": 33.89667, "lon": -88.32889, "elev": 476}, + "KGYX": {"lat": 43.89139, "lon": -70.25694, "elev": 409}, + "KHDX": {"lat": 33.07639, "lon": -106.12222, "elev": 4222}, + "KHGX": {"lat": 29.47194, "lon": -95.07889, "elev": 18}, + "KHNX": {"lat": 36.31417, "lon": -119.63111, "elev": 243}, + "KHPX": {"lat": 36.73667, "lon": -87.285, "elev": 576}, + "KHTX": {"lat": 34.93056, "lon": -86.08361, "elev": 1760}, + "KICT": {"lat": 37.65444, "lon": -97.4425, "elev": 1335}, + "KICX": {"lat": 37.59083, "lon": -112.86222, "elev": 10600}, + "KILN": {"lat": 39.42028, "lon": -83.82167, "elev": 1056}, + "KILX": {"lat": 40.15056, "lon": -89.33667, "elev": 582}, + "KIND": {"lat": 39.7075, "lon": -86.28028, "elev": 790}, + "KINX": {"lat": 36.175, "lon": -95.56444, "elev": 668}, + "KIWA": {"lat": 33.28917, "lon": -111.66917, "elev": 1353}, + "KIWX": {"lat": 41.40861, "lon": -85.7, "elev": 960}, + "KJAX": {"lat": 30.48444, "lon": -81.70194, "elev": 33}, + "KJGX": {"lat": 32.675, "lon": -83.35111, "elev": 521}, + "KJKL": {"lat": 37.59083, "lon": -83.31306, "elev": 1364}, + "KLBB": {"lat": 33.65417, "lon": -101.81361, "elev": 3259}, + "KLCH": {"lat": 30.125, "lon": -93.21583, "elev": 13}, + "KLGX": {"lat": 47.1158, "lon": -124.1069, "elev": 252}, + "KLIX": {"lat": 30.33667, "lon": -89.82528, "elev": 24}, + "KLNX": {"lat": 41.95778, "lon": -100.57583, "elev": 2970}, + "KLOT": {"lat": 41.60444, "lon": -88.08472, "elev": 663}, + "KLRX": {"lat": 40.73972, "lon": -116.80278, "elev": 6744}, + "KLSX": {"lat": 38.69889, "lon": -90.68278, "elev": 608}, + "KLTX": {"lat": 33.98917, "lon": -78.42917, "elev": 64}, + "KLVX": {"lat": 37.97528, "lon": -85.94389, "elev": 719}, + "KLWX": {"lat": 38.97628, "lon": -77.48751, "elev": -99999}, + "KLZK": {"lat": 34.83639, "lon": -92.26194, "elev": 568}, + "KMAF": {"lat": 31.94333, "lon": -102.18889, "elev": 2868}, + "KMAX": {"lat": 42.08111, "lon": -122.71611, "elev": 7513}, + "KMBX": {"lat": 48.3925, "lon": -100.86444, "elev": 1493}, + "KMHX": {"lat": 34.77583, "lon": -76.87639, "elev": 31}, + "KMKX": {"lat": 42.96778, "lon": -88.55056, "elev": 958}, + "KMLB": {"lat": 28.11306, "lon": -80.65444, "elev": 99}, + "KMOB": {"lat": 30.67944, "lon": -88.23972, "elev": 208}, + "KMPX": {"lat": 44.84889, "lon": -93.56528, "elev": 946}, + "KMQT": {"lat": 46.53111, "lon": -87.54833, "elev": 1411}, + "KMRX": {"lat": 36.16833, "lon": -83.40194, "elev": 1337}, + "KMSX": {"lat": 47.04111, "lon": -113.98611, "elev": 7855}, + "KMTX": {"lat": 41.26278, "lon": -112.44694, "elev": 6460}, + "KMUX": {"lat": 37.15528, "lon": -121.8975, "elev": 3469}, + "KMVX": {"lat": 47.52806, "lon": -97.325, "elev": 986}, + "KMXX": {"lat": 32.53667, "lon": -85.78972, "elev": 400}, + "KNKX": {"lat": 32.91889, "lon": -117.04194, "elev": 955}, + "KNQA": {"lat": 35.34472, "lon": -89.87333, "elev": 282}, + "KOAX": {"lat": 41.32028, "lon": -96.36639, "elev": 1148}, + "KOHX": {"lat": 36.24722, "lon": -86.5625, "elev": 579}, + "KOKX": {"lat": 40.86556, "lon": -72.86444, "elev": 85}, + "KOTX": {"lat": 47.68056, "lon": -117.62583, "elev": 2384}, + "KPAH": {"lat": 37.06833, "lon": -88.77194, "elev": 392}, + "KPBZ": {"lat": 40.53167, "lon": -80.21833, "elev": 1185}, + "KPDT": {"lat": 45.69056, "lon": -118.85278, "elev": 1515}, + "KPOE": {"lat": 31.15528, "lon": -92.97583, "elev": 408}, + "KPUX": {"lat": 38.45944, "lon": -104.18139, "elev": 5249}, + "KRAX": {"lat": 35.66528, "lon": -78.49, "elev": 348}, + "KRGX": {"lat": 39.75417, "lon": -119.46111, "elev": 8299}, + "KRIW": {"lat": 43.06611, "lon": -108.47667, "elev": 5568}, + "KRLX": {"lat": 38.31194, "lon": -81.72389, "elev": 1080}, + "KRTX": {"lat": 45.715, "lon": -122.96417, "elev": -99999}, + "KSFX": {"lat": 43.10583, "lon": -112.68528, "elev": 4474}, + "KSGF": {"lat": 37.23528, "lon": -93.40028, "elev": 1278}, + "KSHV": {"lat": 32.45056, "lon": -93.84111, "elev": 273}, + "KSJT": {"lat": 31.37111, "lon": -100.49222, "elev": 1890}, + "KSOX": {"lat": 33.81778, "lon": -117.635, "elev": 3027}, + "KSRX": {"lat": 35.29056, "lon": -94.36167, "elev": -99999}, + "KTBW": {"lat": 27.70528, "lon": -82.40194, "elev": 41}, + "KTFX": {"lat": 47.45972, "lon": -111.38444, "elev": 3714}, + "KTLH": {"lat": 30.3975, "lon": -84.32889, "elev": 63}, + "KTLX": {"lat": 35.33306, "lon": -97.2775, "elev": 1213}, + "KTWX": {"lat": 38.99694, "lon": -96.2325, "elev": 1367}, + "KTYX": {"lat": 43.75583, "lon": -75.68, "elev": 1846}, + "KUDX": {"lat": 44.125, "lon": -102.82944, "elev": 3016}, + "KUEX": {"lat": 40.32083, "lon": -98.44167, "elev": 1976}, + "KVAX": {"lat": 30.89, "lon": -83.00194, "elev": 178}, + "KVBX": {"lat": 34.83806, "lon": -120.39583, "elev": 1233}, + "KVNX": {"lat": 36.74083, "lon": -98.1275, "elev": 1210}, + "KVTX": {"lat": 34.41167, "lon": -119.17861, "elev": 2726}, + "KVWX": {"lat": 38.2600, "lon": -87.7247, "elev": -99999}, + "KYUX": {"lat": 32.49528, "lon": -114.65583, "elev": 174}, + "LPLA": {"lat": 38.73028, "lon": -27.32167, "elev": 3334}, + "PABC": {"lat": 60.79278, "lon": -161.87417, "elev": 162}, + "PACG": {"lat": 56.85278, "lon": -135.52917, "elev": 270}, + "PAEC": {"lat": 64.51139, "lon": -165.295, "elev": 54}, + "PAHG": {"lat": 60.725914, "lon": -151.35146, "elev": 243}, + "PAIH": {"lat": 59.46194, "lon": -146.30111, "elev": 67}, + "PAKC": {"lat": 58.67944, "lon": -156.62944, "elev": 63}, + "PAPD": {"lat": 65.03556, "lon": -147.49917, "elev": 2593}, + "PGUA": {"lat": 13.45444, "lon": 144.80833, "elev": 264}, + "PHKI": {"lat": 21.89417, "lon": -159.55222, "elev": 179}, + "PHKM": {"lat": 20.12556, "lon": -155.77778, "elev": 3812}, + "PHMO": {"lat": 21.13278, "lon": -157.18, "elev": 1363}, + "PHWA": {"lat": 19.095, "lon": -155.56889, "elev": 1370}, + "RKJK": {"lat": 35.92417, "lon": 126.62222, "elev": 78}, + "RKSG": {"lat": 36.95972, "lon": 127.01833, "elev": 52}, + "RODN": {"lat": 26.30194, "lon": 127.90972, "elev": 218}, + "TJUA": {"lat": 18.1175, "lon": -66.07861, "elev": 2794}, + "TJFK": {"lat": 40.5668, "lon": -73.8874, "elev": 112}, + "TADW": {"lat": 38.6704, "lon": -76.8446, "elev": 346}, + "TATL": {"lat": 33.6433, "lon": -84.2524, "elev": 1075}, + "TBNA": {"lat": 35.9767, "lon": -86.6618, "elev": 817}, + "TBOS": {"lat": 42.1515, "lon": -70.9302, "elev": 264}, + "TBWI": {"lat": 39.0870, "lon": -76.6276, "elev": 297}, + "TCLT": {"lat": 35.3269, "lon": -80.8772, "elev": 871}, + "TCMH": {"lat": 39.9878, "lon": -82.71, "elev": 1148}, + "TCVG": {"lat": 38.8799, "lon": -84.5737, "elev": 1053}, + "TDAL": {"lat": 32.9076, "lon": -96.9568, "elev": 622}, + "TDAY": {"lat": 39.9875, "lon": -84.1102, "elev": 1019}, + "TDCA": {"lat": 38.7474, "lon": -76.9509, "elev": 345}, + "TDEN": {"lat": 39.7256, "lon": -104.5431, "elev": 5701}, + "TDFW": {"lat": 33.0396, "lon": -96.8974, "elev": 585}, + "TDTW": {"lat": 42.0710, "lon": -83.4704, "elev": 772}, + "TEWR": {"lat": 40.5880, "lon": -74.2503, "elev": 136}, + "TFLL": {"lat": 26.1263, "lon": -80.3478, "elev": 120}, + "THOU": {"lat": 29.5328, "lon": -95.2444, "elev": 117}, + "TIAD": {"lat": 39.0675, "lon": -77.5012, "elev": 473}, + "TIAH": {"lat": 30.0297, "lon": -95.5708, "elev": 253}, + "TICH": {"lat": 37.4069, "lon": -97.4764, "elev": 1351}, + "TIDS": {"lat": 39.5978, "lon": -86.4085, "elev": 847}, + "TLAS": {"lat": 36.1292, "lon": -115.0147, "elev": 2058}, + "TLVE": {"lat": 41.2805, "lon": -81.9659, "elev": 931}, + "TMCI": {"lat": 39.4488, "lon": -94.7396, "elev": 1090}, + "TMCO": {"lat": 28.2584, "lon": -81.3133, "elev": 169}, + "TMDW": {"lat": 41.69, "lon": -87.8034, "elev": 763}, + "TMEM": {"lat": 34.8867, "lon": -90.0007, "elev": 483}, + "TMIA": {"lat": 25.7555, "lon": -80.4932, "elev": 125}, + "TMKE": {"lat": 42.7619, "lon": -87.9994, "elev": 933}, + "TMSP": {"lat": 44.8197, "lon": -92.9392, "elev": 1121}, + "TMSY": {"lat": 29.9385, "lon": -90.3811, "elev": 99}, + "TOKC": {"lat": 35.2474, "lon": -97.5395, "elev": 1308}, + "TORD": {"lat": 41.7712, "lon": -87.8363, "elev": 744}, + "TPBI": {"lat": 26.6572, "lon": -80.2586, "elev": 133}, + "TPHL": {"lat": 39.9084, "lon": -75.0426, "elev": 153}, + "TPHX": {"lat": 33.3678, "lon": -112.1580, "elev": 1089}, + "TPIT": {"lat": 40.4641, "lon": -80.4697, "elev": 1386}, + "TRDU": {"lat": 35.9898, "lon": -78.6787, "elev": 515}, + "TSDF": {"lat": 38.0109, "lon": -85.5995, "elev": 731}, + "TSJU": {"lat": 18.4313, "lon": -66.1722, "elev": 157}, + "TSLC": {"lat": 40.9341, "lon": -111.9214, "elev": 4295}, + "TSTL": {"lat": 38.7668, "lon": -90.4698, "elev": 647}, + "TTPA": {"lat": 27.8196, "lon": -82.5179, "elev": 93}, + "TTUL": {"lat": 36.0236, "lon": -95.8175, "elev": 823}, +} diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py new file mode 100644 index 00000000..c4ae52b1 --- /dev/null +++ b/xradar/io/backends/nexrad_level2.py @@ -0,0 +1,1381 @@ +import bz2 +import struct +import warnings +from datetime import datetime, timedelta + +import numpy as np +import xarray as xr +from datatree import DataTree +from xarray.backends.common import BackendEntrypoint +from xarray.core import indexing +from xarray.core.variable import Variable + +from xradar import util +from xradar.io.backends.common import ( + LazyLoadDict, + _assign_root, + _attach_sweep_groups, + prepare_for_read, +) +from xradar.io.backends.nexrad_common import get_nexrad_location +from xradar.model import ( + get_altitude_attrs, + get_azimuth_attrs, + get_elevation_attrs, + get_latitude_attrs, + get_longitude_attrs, + get_moment_attrs, + get_range_attrs, + get_time_attrs, +) + +nexrad_mapping = { + "REF": "DBZH", + "VEL": "VRADH", + "SW": "WRADH", + "ZDR": "ZDR", + "PHI": "PHIDP", + "RHO": "RHOHV", + "CFP": "CCORH", +} + + +class _NEXRADLevel2StagedField: + """ + A class to facilitate on demand loading of field data from a Level 2 file. + """ + + def __init__(self, nfile, moment, max_ngates, scans): + """initialize.""" + self.nfile = nfile + self.moment = moment + self.max_ngates = max_ngates + self.scans = scans + + def __call__(self): + """Return the array containing the field data.""" + return self.nfile.get_data(self.moment, self.max_ngates, scans=self.scans) + + +class NEXRADLevel2File: + """ + 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. + + Attributes + ---------- + radial_records : list + Radial (1 or 31) messages in the file. + nscans : int + Number of scans in the file. + scan_msgs : list of arrays + Each element specifies the indices of the message in the + radial_records attribute which belong to a given scan. + volume_header : dict + Volume header. + vcp : dict + VCP information dictionary. + _records : list + A list of all records (message) in the file. + _fh : file-like + File like object from which data is read. + _msg_type : '31' or '1': + Type of radial messages in file. + + 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): + """initalize the object.""" + # read in the volume header and compression_record + if hasattr(filename, "read"): + fh = filename + else: + fh = open(filename, "rb") + size = _structure_size(VOLUME_HEADER) + self.volume_header = _unpack_structure(fh.read(size), VOLUME_HEADER) + compression_record = fh.read(COMPRESSION_RECORD_SIZE) + + # read the records in the file, decompressing as needed + compression_slice = slice(CONTROL_WORD_SIZE, CONTROL_WORD_SIZE + 2) + compression_or_ctm_info = compression_record[compression_slice] + if compression_or_ctm_info == b"BZ": + buf = _decompress_records(fh) + # The 12-byte compression record previously held the Channel Terminal + # Manager (CTM) information. Bytes 4 through 6 contain the size of the + # record (2432) as a big endian unsigned short, which is encoded as + # b'\t\x80' == struct.pack('>H', 2432). + # Newer files zero out this section. + elif compression_or_ctm_info in (b"\x00\x00", b"\t\x80"): + buf = fh.read() + else: + raise OSError("unknown compression record") + self._fh = fh + + # read the records from the buffer + self._records = [] + buf_length = len(buf) + pos = 0 + while pos < buf_length: + pos, dic = _get_record_from_buf(buf, pos) + self._records.append(dic) + + # pull out radial records (1 or 31) which contain the moment data. + self.radial_records = [r for r in self._records if r["header"]["type"] == 31] + self._msg_type = "31" + if len(self.radial_records) == 0: + self.radial_records = [r for r in self._records if r["header"]["type"] == 1] + self._msg_type = "1" + if len(self.radial_records) == 0: + raise ValueError("No MSG31 records found, cannot read file") + elev_nums = np.array( + [m["msg_header"]["elevation_number"] for m in self.radial_records] + ) + self.scan_msgs = [ + np.where(elev_nums == i + 1)[0] for i in range(elev_nums.max()) + ] + self.nscans = len(self.scan_msgs) + + # pull out the vcp record + msg_5 = [r for r in self._records if r["header"]["type"] == 5] + + if len(msg_5): + self.vcp = msg_5[0] + else: + # There is no VCP Data.. This is uber dodgy + warnings.warn( + "No MSG5 detected. Setting to meaningless data. " + "Rethink your life choices and be ready for errors." + "Specifically fixed angle data will be missing" + ) + + self.vcp = None + return + + def close(self): + """Close the file.""" + self._fh.close() + + def location(self): + """ + Find the location of the radar. + + Returns all zeros if location is not available. + + Returns + ------- + latitude : float + Latitude of the radar in degrees. + longitude : float + Longitude of the radar in degrees. + height : int + Height of radar and feedhorn in meters above mean sea level. + + """ + if self._msg_type == "31": + dic = self.radial_records[0]["VOL"] + height = dic["height"] + dic["feedhorn_height"] + return dic["lat"], dic["lon"], height + else: + return 0.0, 0.0, 0.0 + + def scan_info(self, scans=None): + """ + Return a list of dictionaries with scan information. + + Parameters + ---------- + scans : list ot None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + scan_info : list, optional + A list of the scan performed with a dictionary with keys + 'moments', 'ngates', 'nrays', 'first_gate' and 'gate_spacing' + for each scan. The 'moments', 'ngates', 'first_gate', and + 'gate_spacing' keys are lists of the NEXRAD moments and gate + information for that moment collected during the specific scan. + The 'nrays' key provides the number of radials collected in the + given scan. + + """ + info = [] + + if scans is None: + scans = range(self.nscans) + for scan in scans: + nrays = self.get_nrays(scan) + if nrays < 2: + self.nscans -= 1 + continue + msg31_number = self.scan_msgs[scan][0] + msg = self.radial_records[msg31_number] + nexrad_moments = ["REF", "VEL", "SW", "ZDR", "PHI", "RHO", "CFP"] + moments = [f for f in nexrad_moments if f in msg] + ngates = [msg[f]["ngates"] for f in moments] + gate_spacing = [msg[f]["gate_spacing"] for f in moments] + first_gate = [msg[f]["first_gate"] for f in moments] + info.append( + { + "nrays": nrays, + "ngates": ngates, + "gate_spacing": gate_spacing, + "first_gate": first_gate, + "moments": moments, + } + ) + return info + + def get_vcp_pattern(self): + """ + Return the numerical volume coverage pattern (VCP) or None if unknown. + """ + if self.vcp is None: + return None + else: + return self.vcp["msg5_header"]["pattern_number"] + + def get_nrays(self, scan): + """ + Return the number of rays in a given scan. + + Parameters + ---------- + scan : int + Scan of interest (0 based). + + Returns + ------- + nrays : int + Number of rays (radials) in the scan. + + """ + return len(self.scan_msgs[scan]) + + def get_range(self, scan_num, moment): + """ + Return an array of gate ranges for a given scan and moment. + + Parameters + ---------- + scan_num : int + Scan number (0 based). + moment : 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'RHO', or 'CFP' + Moment of interest. + + Returns + ------- + range : ndarray + Range in meters from the antenna to the center of gate (bin). + + """ + dic = self.radial_records[self.scan_msgs[scan_num][0]][moment] + ngates = dic["ngates"] + first_gate = dic["first_gate"] + gate_spacing = dic["gate_spacing"] + return np.arange(ngates) * gate_spacing + first_gate + + # helper functions for looping over scans + def _msg_nums(self, scans): + """Find the all message number for a list of scans.""" + return np.concatenate([self.scan_msgs[i] for i in scans]) + + def _radial_array(self, scans, key): + """ + Return an array of radial header elements for all rays in scans. + """ + msg_nums = self._msg_nums(scans) + temp = [self.radial_records[i]["msg_header"][key] for i in msg_nums] + return np.array(temp) + + def _radial_sub_array(self, scans, key): + """ + Return an array of RAD or msg_header elements for all rays in scans. + """ + msg_nums = self._msg_nums(scans) + if self._msg_type == "31": + tmp = [self.radial_records[i]["RAD"][key] for i in msg_nums] + else: + tmp = [self.radial_records[i]["msg_header"][key] for i in msg_nums] + return np.array(tmp) + + def get_times(self, scans=None): + """ + Retrieve the times at which the rays were collected. + + Parameters + ---------- + scans : list or None + Scans (0-based) to retrieve ray (radial) collection times from. + None (the default) will return the times for all scans in the + volume. + + Returns + ------- + time_start : Datetime + Initial time. + time : ndarray + Offset in seconds from the initial time at which the rays + in the requested scans were collected. + + """ + if scans is None: + scans = range(self.nscans) + days = self._radial_array(scans, "collect_date") + secs = self._radial_array(scans, "collect_ms") / 1000.0 + offset = timedelta(days=int(days[0]) - 1, seconds=int(secs[0])) + time_start = datetime(1970, 1, 1) + offset + time = secs - int(secs[0]) + (days - days[0]) * 86400 + return time_start, time + + def get_azimuth_angles(self, scans=None): + """ + Retrieve the azimuth angles of all rays in the requested scans. + + Parameters + ---------- + scans : list ot None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + angles : ndarray + Azimuth angles in degress for all rays in the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "1": + scale = 180 / (4096 * 8.0) + else: + scale = 1.0 + return self._radial_array(scans, "azimuth_angle") * scale + + def get_elevation_angles(self, scans=None): + """ + Retrieve the elevation angles of all rays in the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for + all scans in the volume. + + Returns + ------- + angles : ndarray + Elevation angles in degress for all rays in the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "1": + scale = 180 / (4096 * 8.0) + else: + scale = 1.0 + return self._radial_array(scans, "elevation_angle") * scale + + def get_target_angles(self, scans=None): + """ + Retrieve the target elevation angle of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the target elevation angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + angles : ndarray + Target elevation angles in degress for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "31": + if self.vcp is not None: + cut_parameters = self.vcp["cut_parameters"] + else: + cut_parameters = [{"elevation_angle": 0.0}] * self.nscans + scale = 360.0 / 65536.0 + return np.array( + [cut_parameters[i]["elevation_angle"] * scale for i in scans], + dtype="float32", + ) + else: + scale = 180 / (4096 * 8.0) + msgs = [self.radial_records[self.scan_msgs[i][0]] for i in scans] + return np.round( + np.array( + [m["msg_header"]["elevation_angle"] * scale for m in msgs], + dtype="float32", + ), + 1, + ) + + def get_nyquist_vel(self, scans=None): + """ + Retrieve the Nyquist velocities of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the Nyquist velocities will be + retrieved. None (the default) will return the velocities for all + scans in the volume. + + Returns + ------- + velocities : ndarray + Nyquist velocities (in m/s) for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + return self._radial_sub_array(scans, "nyquist_vel") * 0.01 + + def get_unambigous_range(self, scans=None): + """ + Retrieve the unambiguous range of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the unambiguous range will be retrieved. + None (the default) will return the range for all scans in the + volume. + + Returns + ------- + unambiguous_range : ndarray + Unambiguous range (in meters) for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + # unambiguous range is stored in tenths of km, x100 for meters + return self._radial_sub_array(scans, "unambig_range") * 100.0 + + def get_data(self, moment, max_ngates, scans=None, raw_data=False): + """ + Retrieve moment data for a given set of scans. + + Masked points indicate that the data was not collected, below + threshold or is range folded. + + Parameters + ---------- + moment : 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'RHO', or 'CFP' + Moment for which to to retrieve data. + max_ngates : int + Maximum number of gates (bins) in any ray. + requested. + raw_data : bool + True to return the raw data, False to perform masking as well as + applying the appropiate scale and offset to the data. When + raw_data is True values of 1 in the data likely indicate that + the gate was not present in the sweep, in some cases in will + indicate range folded data. + scans : list or None. + Scans to retrieve data from (0 based). None (the default) will + get the data for all scans in the volume. + + Returns + ------- + data : ndarray + + """ + if scans is None: + scans = range(self.nscans) + + # determine the number of rays + msg_nums = self._msg_nums(scans) + nrays = len(msg_nums) + # extract the data + set_datatype = False + data = np.ones((nrays, max_ngates), ">B") + for i, msg_num in enumerate(msg_nums): + msg = self.radial_records[msg_num] + if moment not in msg.keys(): + continue + if not set_datatype: + data = data.astype(">" + _bits_to_code(msg, moment)) + set_datatype = True + + ngates = min(msg[moment]["ngates"], max_ngates, len(msg[moment]["data"])) + data[i, :ngates] = msg[moment]["data"][:ngates] + # return raw data if requested + if raw_data: + return data + + # mask, scan and offset, assume that the offset and scale + # are the same in all scans/gates + for scan in scans: # find a scan which contains the moment + msg_num = self.scan_msgs[scan][0] + msg = self.radial_records[msg_num] + if moment in msg.keys(): + offset = np.float32(msg[moment]["offset"]) + scale = np.float32(msg[moment]["scale"]) + mask = data <= 1 + scaled_data = (data - offset) / scale + return np.ma.array(scaled_data, mask=mask) + + # moment is not present in any scan, mask all values + return np.ma.masked_less_equal(data, 1) + + +def _bits_to_code(msg, moment): + """ + Convert number of bits to the proper code for unpacking. + Based on the code found in MetPy: + https://github.com/Unidata/MetPy/blob/40d5c12ab341a449c9398508bd41 + d010165f9eeb/src/metpy/io/_tools.py#L313-L321 + """ + if msg["header"]["type"] == 1: + word_size = msg[moment]["data"].dtype + if word_size == "uint16": + return "H" + elif word_size == "uint8": + return "B" + else: + warnings.warn(('Unsupported bit size: %s. Returning "B"', word_size)) + return "B" + + elif msg["header"]["type"] == 31: + word_size = msg[moment]["word_size"] + if word_size == 16: + return "H" + elif word_size == 8: + return "B" + else: + warnings.warn(('Unsupported bit size: %s. Returning "B"', word_size)) + return "B" + else: + raise TypeError("Unsupported msg type %s", msg["header"]["type"]) + + +def _decompress_records(file_handler): + """ + Decompressed the records from an BZ2 compressed Archive 2 file. + """ + file_handler.seek(0) + cbuf = file_handler.read() # read all data from the file + decompressor = bz2.BZ2Decompressor() + skip = _structure_size(VOLUME_HEADER) + CONTROL_WORD_SIZE + buf = decompressor.decompress(cbuf[skip:]) + while len(decompressor.unused_data): + cbuf = decompressor.unused_data + decompressor = bz2.BZ2Decompressor() + buf += decompressor.decompress(cbuf[CONTROL_WORD_SIZE:]) + + return buf[COMPRESSION_RECORD_SIZE:] + + +def _get_record_from_buf(buf, pos): + """Retrieve and unpack a NEXRAD record from a buffer.""" + dic = {"header": _unpack_from_buf(buf, pos, MSG_HEADER)} + msg_type = dic["header"]["type"] + + if msg_type == 31: + new_pos = _get_msg31_from_buf(buf, pos, dic) + elif msg_type == 5: + # Sometimes we encounter incomplete buffers + try: + new_pos = _get_msg5_from_buf(buf, pos, dic) + except struct.error: + warnings.warn( + "Encountered incomplete MSG5. File may be corrupt.", RuntimeWarning + ) + new_pos = pos + RECORD_SIZE + elif msg_type == 29: + new_pos = _get_msg29_from_buf(pos, dic) + warnings.warn("Message 29 encountered, not parsing.", RuntimeWarning) + elif msg_type == 1: + new_pos = _get_msg1_from_buf(buf, pos, dic) + else: # not message 31 or 1, no decoding performed + new_pos = pos + RECORD_SIZE + + return new_pos, dic + + +def _get_msg29_from_buf(pos, dic): + msg_size = dic["header"]["size"] + if msg_size == 65535: + msg_size = dic["header"]["segments"] << 16 | dic["header"]["seg_num"] + msg_header_size = _structure_size(MSG_HEADER) + new_pos = pos + msg_header_size + msg_size + return new_pos + + +def _get_msg31_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG31 record from a buffer.""" + msg_size = dic["header"]["size"] * 2 - 4 + msg_header_size = _structure_size(MSG_HEADER) + new_pos = pos + msg_header_size + msg_size + mbuf = buf[pos + msg_header_size : new_pos] + msg_31_header = _unpack_from_buf(mbuf, 0, MSG_31) + block_pointers = [ + v for k, v in msg_31_header.items() if k.startswith("block_pointer") and v > 0 + ] + for block_pointer in block_pointers: + block_name, block_dic = _get_msg31_data_block(mbuf, block_pointer) + dic[block_name] = block_dic + + dic["msg_header"] = msg_31_header + return new_pos + + +def _get_msg31_data_block(buf, ptr): + """Unpack a msg_31 data block into a dictionary.""" + block_name = buf[ptr + 1 : ptr + 4].decode("ascii").strip() + + if block_name == "VOL": + dic = _unpack_from_buf(buf, ptr, VOLUME_DATA_BLOCK) + elif block_name == "ELV": + dic = _unpack_from_buf(buf, ptr, ELEVATION_DATA_BLOCK) + elif block_name == "RAD": + dic = _unpack_from_buf(buf, ptr, RADIAL_DATA_BLOCK) + elif block_name in ["REF", "VEL", "SW", "ZDR", "PHI", "RHO", "CFP"]: + dic = _unpack_from_buf(buf, ptr, GENERIC_DATA_BLOCK) + ngates = dic["ngates"] + ptr2 = ptr + _structure_size(GENERIC_DATA_BLOCK) + if dic["word_size"] == 16: + data = np.frombuffer(buf[ptr2 : ptr2 + ngates * 2], ">u2") + elif dic["word_size"] == 8: + data = np.frombuffer(buf[ptr2 : ptr2 + ngates], ">u1") + else: + warnings.warn( + 'Unsupported bit size: %s. Returning array dtype "B"', dic["word_size"] + ) + dic["data"] = data + else: + dic = {} + return block_name, dic + + +def _get_msg1_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG1 record from a buffer.""" + msg_header_size = _structure_size(MSG_HEADER) + msg1_header = _unpack_from_buf(buf, pos + msg_header_size, MSG_1) + dic["msg_header"] = msg1_header + + sur_nbins = int(msg1_header["sur_nbins"]) + doppler_nbins = int(msg1_header["doppler_nbins"]) + + sur_step = int(msg1_header["sur_range_step"]) + doppler_step = int(msg1_header["doppler_range_step"]) + + sur_first = int(msg1_header["sur_range_first"]) + doppler_first = int(msg1_header["doppler_range_first"]) + if doppler_first > 2**15: + doppler_first = doppler_first - 2**16 + + if msg1_header["sur_pointer"]: + offset = pos + msg_header_size + msg1_header["sur_pointer"] + data = np.frombuffer(buf[offset : offset + sur_nbins], ">u1") + dic["REF"] = { + "ngates": sur_nbins, + "gate_spacing": sur_step, + "first_gate": sur_first, + "data": data, + "scale": 2.0, + "offset": 66.0, + } + if msg1_header["vel_pointer"]: + offset = pos + msg_header_size + msg1_header["vel_pointer"] + data = np.frombuffer(buf[offset : offset + doppler_nbins], ">u1") + dic["VEL"] = { + "ngates": doppler_nbins, + "gate_spacing": doppler_step, + "first_gate": doppler_first, + "data": data, + "scale": 2.0, + "offset": 129.0, + } + if msg1_header["doppler_resolution"] == 4: + # 1 m/s resolution velocity, offset remains 129. + dic["VEL"]["scale"] = 1.0 + if msg1_header["width_pointer"]: + offset = pos + msg_header_size + msg1_header["width_pointer"] + data = np.frombuffer(buf[offset : offset + doppler_nbins], ">u1") + dic["SW"] = { + "ngates": doppler_nbins, + "gate_spacing": doppler_step, + "first_gate": doppler_first, + "data": data, + "scale": 2.0, + "offset": 129.0, + } + return pos + RECORD_SIZE + + +def _get_msg5_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG1 record from a buffer.""" + msg_header_size = _structure_size(MSG_HEADER) + msg5_header_size = _structure_size(MSG_5) + msg5_elev_size = _structure_size(MSG_5_ELEV) + + dic["msg5_header"] = _unpack_from_buf(buf, pos + msg_header_size, MSG_5) + dic["cut_parameters"] = [] + for i in range(dic["msg5_header"]["num_cuts"]): + pos2 = pos + msg_header_size + msg5_header_size + msg5_elev_size * i + dic["cut_parameters"].append(_unpack_from_buf(buf, pos2, MSG_5_ELEV)) + return pos + RECORD_SIZE + + +def _structure_size(structure): + """Find the size of a structure in bytes.""" + return struct.calcsize(">" + "".join([i[1] for i in structure])) + + +def _unpack_from_buf(buf, pos, structure): + """Unpack a structure from a buffer.""" + size = _structure_size(structure) + return _unpack_structure(buf[pos : pos + size], structure) + + +def _unpack_structure(string, structure): + """Unpack a structure from a string.""" + fmt = ">" + "".join([i[1] for i in structure]) # NEXRAD is big-endian + lst = struct.unpack(fmt, string) + return dict(zip([i[0] for i in structure], lst)) + + +# 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_SIZE = 2432 +COMPRESSION_RECORD_SIZE = 12 +CONTROL_WORD_SIZE = 4 + +# format of structure elements +# section 3.2.1, page 3-2 +CODE1 = "B" +CODE2 = "H" +INT1 = "B" +INT2 = "H" +INT4 = "I" +REAL4 = "f" +REAL8 = "d" +SINT1 = "b" +SINT2 = "h" +SINT4 = "i" + +# Figure 1 in Interface Control Document for the Archive II/User +# page 7-2 +VOLUME_HEADER = ( + ("tape", "9s"), + ("extension", "3s"), + ("date", "I"), + ("time", "I"), + ("icao", "4s"), +) + +# Table II Message Header Data +# page 3-7 +MSG_HEADER = ( + ("size", INT2), # size of data, no including header + ("channels", INT1), + ("type", INT1), + ("seq_id", INT2), + ("date", INT2), + ("ms", INT4), + ("segments", INT2), + ("seg_num", INT2), +) + +# Table XVII Digital Radar Generic Format Blocks (Message Type 31) +# pages 3-87 to 3-89 +MSG_31 = ( + ("id", "4s"), # 0-3 + ("collect_ms", INT4), # 4-7 + ("collect_date", INT2), # 8-9 + ("azimuth_number", INT2), # 10-11 + ("azimuth_angle", REAL4), # 12-15 + ("compress_flag", CODE1), # 16 + ("spare_0", INT1), # 17 + ("radial_length", INT2), # 18-19 + ("azimuth_resolution", CODE1), # 20 + ("radial_spacing", CODE1), # 21 + ("elevation_number", INT1), # 22 + ("cut_sector", INT1), # 23 + ("elevation_angle", REAL4), # 24-27 + ("radial_blanking", CODE1), # 28 + ("azimuth_mode", SINT1), # 29 + ("block_count", INT2), # 30-31 + ("block_pointer_1", INT4), # 32-35 Volume Data Constant XVII-E + ("block_pointer_2", INT4), # 36-39 Elevation Data Constant XVII-F + ("block_pointer_3", INT4), # 40-43 Radial Data Constant XVII-H + ("block_pointer_4", INT4), # 44-47 Moment "REF" XVII-{B/I} + ("block_pointer_5", INT4), # 48-51 Moment "VEL" + ("block_pointer_6", INT4), # 52-55 Moment "SW" + ("block_pointer_7", INT4), # 56-59 Moment "ZDR" + ("block_pointer_8", INT4), # 60-63 Moment "PHI" + ("block_pointer_9", INT4), # 64-67 Moment "RHO" + ("block_pointer_10", INT4), # Moment "CFP" +) + + +# Table III Digital Radar Data (Message Type 1) +# pages 3-7 to +MSG_1 = ( + ("collect_ms", INT4), # 0-3 + ("collect_date", INT2), # 4-5 + ("unambig_range", SINT2), # 6-7 + ("azimuth_angle", CODE2), # 8-9 + ("azimuth_number", INT2), # 10-11 + ("radial_status", CODE2), # 12-13 + ("elevation_angle", INT2), # 14-15 + ("elevation_number", INT2), # 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", INT2), # 26-27 + ("doppler_nbins", INT2), # 28-29 + ("cut_sector_num", INT2), # 30-31 + ("calib_const", REAL4), # 32-35 + ("sur_pointer", INT2), # 36-37 + ("vel_pointer", INT2), # 38-39 + ("width_pointer", INT2), # 40-41 + ("doppler_resolution", CODE2), # 42-43 + ("vcp", INT2), # 44-45 + ("spare_1", "8s"), # 46-53 + ("spare_2", "2s"), # 54-55 + ("spare_3", "2s"), # 56-57 + ("spare_4", "2s"), # 58-59 + ("nyquist_vel", SINT2), # 60-61 + ("atmos_attenuation", SINT2), # 62-63 + ("threshold", SINT2), # 64-65 + ("spot_blank_status", INT2), # 66-67 + ("spare_5", "32s"), # 68-99 + # 100+ reflectivity, velocity and/or spectral width data, CODE1 +) + +# Table XI Volume Coverage Pattern Data (Message Type 5 & 7) +# pages 3-51 to 3-54 +MSG_5 = ( + ("msg_size", INT2), + ("pattern_type", CODE2), + ("pattern_number", INT2), + ("num_cuts", INT2), + ("clutter_map_group", INT2), + ("doppler_vel_res", CODE1), # 2: 0.5 degrees, 4: 1.0 degrees + ("pulse_width", CODE1), # 2: short, 4: long + ("spare", "10s"), # halfwords 7-11 (10 bytes, 5 halfwords) +) + +MSG_5_ELEV = ( + ("elevation_angle", CODE2), # scaled by 360/65536 for value in degrees. + ("channel_config", CODE1), + ("waveform_type", CODE1), + ("super_resolution", CODE1), + ("prf_number", INT1), + ("prf_pulse_count", INT2), + ("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", INT2), + ("dop_prf_pulse_count_1", INT2), + ("spare_1", "2s"), + ("edge_angle_2", CODE2), + ("dop_prf_num_2", INT2), + ("dop_prf_pulse_count_2", INT2), + ("spare_2", "2s"), + ("edge_angle_3", CODE2), + ("dop_prf_num_3", INT2), + ("dop_prf_pulse_count_3", INT2), + ("spare_3", "2s"), +) + +# Table XVII-B Data Block (Descriptor of Generic Data Moment Type) +# pages 3-90 and 3-91 +GENERIC_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), # VEL, REF, SW, RHO, PHI, ZDR + ("reserved", INT4), + ("ngates", INT2), + ("first_gate", SINT2), + ("gate_spacing", SINT2), + ("thresh", SINT2), + ("snr_thres", SINT2), + ("flags", CODE1), + ("word_size", INT1), + ("scale", REAL4), + ("offset", REAL4), + # then data +) + +# Table XVII-E Data Block (Volume Data Constant Type) +# page 3-92 +VOLUME_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("version_major", INT1), + ("version_minor", INT1), + ("lat", REAL4), + ("lon", REAL4), + ("height", SINT2), + ("feedhorn_height", INT2), + ("refl_calib", REAL4), + ("power_h", REAL4), + ("power_v", REAL4), + ("diff_refl_calib", REAL4), + ("init_phase", REAL4), + ("vcp", INT2), + ("spare", "2s"), +) + +# Table XVII-F Data Block (Elevation Data Constant Type) +# page 3-93 +ELEVATION_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("atmos", SINT2), + ("refl_calib", REAL4), +) + +# Table XVII-H Data Block (Radial Data Constant Type) +# pages 3-93 +RADIAL_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("unambig_range", SINT2), + ("noise_h", REAL4), + ("noise_v", REAL4), + ("nyquist_vel", SINT2), + ("spare", "2s"), +) + + +def _find_range_params(scan_info): + """Return range parameters, first_gate, gate_spacing, last_gate.""" + min_first_gate = 999999 + min_gate_spacing = 999999 + max_last_gate = 0 + for scan_params in scan_info: + ngates = scan_params["ngates"][0] + for i, moment in enumerate(scan_params["moments"]): + first_gate = scan_params["first_gate"][i] + gate_spacing = scan_params["gate_spacing"][i] + last_gate = first_gate + gate_spacing * (ngates - 0.5) + + min_first_gate = min(min_first_gate, first_gate) + min_gate_spacing = min(min_gate_spacing, gate_spacing) + max_last_gate = max(max_last_gate, last_gate) + return min_first_gate, min_gate_spacing, max_last_gate + + +def format_nexrad_level2_data( + filename, + delay_field_loading=False, + station=None, + scans=None, + storage_options={"anon": True}, + first_dimension=None, + group=None, + **kwargs, +): + # Load the data file in using NEXRADLevel2File Class + nfile = NEXRADLevel2File( + prepare_for_read(filename, storage_options=storage_options) + ) + + # Access the scan info and load in the available moments + scan_info = nfile.scan_info(scans) + available_moments = {m for scan in scan_info for m in scan["moments"]} + + # Deal with time + time_start, _time = nfile.get_times(scans) + time = get_time_attrs(time_start) + time["data"] = _time + + # range + _range = get_range_attrs() + first_gate, gate_spacing, last_gate = _find_range_params(scan_info) + _range["data"] = np.arange(first_gate, last_gate, gate_spacing, "float32") + _range["meters_to_center_of_first_gate"] = float(first_gate) + _range["meters_between_gates"] = float(gate_spacing) + + # metadata + metadata = { + "Conventions": "CF/Radial instrument_parameters", + "version": "1.3", + "title": "", + "institution": "", + "references": "", + "source": "", + "history": "", + "comment": "", + "intrument_name": "", + "original_container": "NEXRAD Level II", + } + + vcp_pattern = nfile.get_vcp_pattern() + if vcp_pattern is not None: + metadata["vcp_pattern"] = vcp_pattern + if "icao" in nfile.volume_header.keys(): + metadata["instrument_name"] = nfile.volume_header["icao"].decode() + + # scan_type + + # latitude, longitude, altitude + latitude = get_latitude_attrs() + longitude = get_longitude_attrs() + altitude = get_altitude_attrs() + + if nfile._msg_type == "1" and station is not None: + lat, lon, alt = get_nexrad_location(station) + elif ( + "icao" in nfile.volume_header.keys() + and nfile.volume_header["icao"].decode()[0] == "T" + ): + lat, lon, alt = get_nexrad_location(nfile.volume_header["icao"].decode()) + else: + lat, lon, alt = nfile.location() + latitude["data"] = np.array(lat, dtype="float64") + longitude["data"] = np.array(lon, dtype="float64") + altitude["data"] = np.array(alt, dtype="float64") + + # Sweep information + sweep_number = { + "units": "count", + "standard_name": "sweep_number", + "long_name": "Sweep number", + } + + sweep_mode = { + "units": "unitless", + "standard_name": "sweep_mode", + "long_name": "Sweep mode", + "comment": 'Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi"', + } + sweep_start_ray_index = { + "long_name": "Index of first ray in sweep, 0-based", + "units": "count", + } + sweep_end_ray_index = { + "long_name": "Index of last ray in sweep, 0-based", + "units": "count", + } + + if scans is None: + nsweeps = int(nfile.nscans) + else: + nsweeps = len(scans) + sweep_number["data"] = np.arange(nsweeps, dtype="int32") + sweep_mode["data"] = np.array(nsweeps * ["azimuth_surveillance"], dtype="S") + + rays_per_scan = [s["nrays"] for s in scan_info] + sweep_end_ray_index["data"] = np.cumsum(rays_per_scan, dtype="int32") - 1 + + rays_per_scan.insert(0, 0) + sweep_start_ray_index["data"] = np.cumsum(rays_per_scan[:-1], dtype="int32") + + # azimuth, elevation, fixed_angle + azimuth = get_azimuth_attrs() + elevation = get_elevation_attrs() + fixed_angle = { + "long_name": "Target angle for sweep", + "units": "degrees", + "standard_name": "target_fixed_angle", + } + + azimuth["data"] = nfile.get_azimuth_angles(scans) + elevation["data"] = nfile.get_elevation_angles(scans).astype("float32") + fixed_agl = [] + for i in nfile.get_target_angles(scans): + if i > 180: + i = i - 360.0 + warnings.warn( + "Fixed_angle(s) greater than 180 degrees present." + " Assuming angle to be negative so subtrating 360", + UserWarning, + ) + else: + i = i + fixed_agl.append(i) + fixed_angles = np.array(fixed_agl, dtype="float32") + fixed_angle["data"] = fixed_angles + + # fields + max_ngates = len(_range["data"]) + + fields = {} + for moment in available_moments: + dic = get_moment_attrs(nexrad_mapping[moment]) + dic["_FillValue"] = -9999 + if delay_field_loading: + dic = LazyLoadDict(dic) + data_call = _NEXRADLevel2StagedField(nfile, moment, max_ngates, scans) + dic.set_lazy("data", data_call) + else: + mdata = nfile.get_data(moment, max_ngates, scans=scans) + dic["data"] = mdata + fields[nexrad_mapping[moment]] = dic + instrument_parameters = {} + + return create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters, + first_dimension=first_dimension, + group=group, + ) + + +def create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters={}, + first_dimension=None, + group=None, +): + """ + Creates an xarray dataset given a selection of fields defined as a collection of dictionaries + """ + if first_dimension is None: + first_dimension = "azimuth" + + if group is not None: + encoding = {"group": group} + group = int(group[6:]) + else: + encoding = {} + group = 0 + + group_slice = slice( + sweep_start_ray_index["data"][group], sweep_end_ray_index["data"][group] + ) + + # Track the number of sweeps + nsweeps = len(sweep_number["data"]) + + # Handle the coordinates + azimuth["data"] = indexing.LazilyOuterIndexedArray(azimuth["data"][group_slice]) + elevation["data"] = indexing.LazilyOuterIndexedArray(elevation["data"][group_slice]) + time["data"] = indexing.LazilyOuterIndexedArray(time["data"][group_slice]) + sweep_mode["data"] = sweep_mode["data"][group] + sweep_number["data"] = sweep_number["data"][group] + fixed_angle["data"] = fixed_angle["data"][group] + + coords = { + "azimuth": dictionary_to_xarray_variable(azimuth, encoding, first_dimension), + "elevation": dictionary_to_xarray_variable( + elevation, encoding, first_dimension + ), + "time": dictionary_to_xarray_variable(time, encoding, first_dimension), + "range": dictionary_to_xarray_variable(_range, encoding, "range"), + "longitude": dictionary_to_xarray_variable(longitude), + "latitude": dictionary_to_xarray_variable(latitude), + "altitude": dictionary_to_xarray_variable(altitude), + "sweep_mode": dictionary_to_xarray_variable(sweep_mode), + "sweep_number": dictionary_to_xarray_variable(sweep_number), + "sweep_fixed_angle": dictionary_to_xarray_variable(fixed_angle), + } + data_vars = { + field: dictionary_to_xarray_variable( + data, dims=(first_dimension, "range"), subset=group_slice + ) + for (field, data) in fields.items() + } + + ds = xr.Dataset(data_vars=data_vars, coords=coords, attrs=metadata) + + ds.attrs["nsweeps"] = nsweeps + + return ds + + +def dictionary_to_xarray_variable(variable_dict, encoding=None, dims=None, subset=None): + attrs = variable_dict.copy() + data = variable_dict.copy() + attrs.pop("data") + if dims is None: + dims = () + if subset is not None: + data["data"] = data["data"][subset] + return Variable(dims, data["data"], encoding, attrs) + + +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, + *, + decode_coords=True, + mask_and_scale=True, + decode_times=True, + concat_characters=True, + drop_variables=None, + use_cftime=None, + decode_timedelta=None, + group=None, + first_dim="auto", + reindex_angle=False, + ): + ds = format_nexrad_level2_data(filename_or_obj, group=group) + + # 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) + + 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) + + return ds + + +def _get_nexradlevel2_number_of_sweeps(filename_or_obj, **kwargs): + number_of_sweeps = xr.open_dataset( + filename_or_obj, engine="nexradlevel2", **kwargs + ).nsweeps + return [f"sweep_{sweep_number}" for sweep_number in np.arange(0, number_of_sweeps)] + + +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: + sweeps = _get_nexradlevel2_number_of_sweeps(filename_or_obj) + + ds = [ + xr.open_dataset(filename_or_obj, group=swp, engine="nexradlevel2", **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:])