From dd7418c23b69c4d3713f2a0b58b348bd1bb12b36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kai=20M=C3=BChlbauer?= Date: Tue, 27 Sep 2022 10:09:04 +0200 Subject: [PATCH] WIP: furuno --- docs/history.md | 4 + docs/importers.md | 18 +++ docs/usage.md | 1 + examples/notebooks/Furuno.ipynb | 209 ++++++++++++++++++++++++++++++++ tests/conftest.py | 22 ++++ tests/test_io.py | 34 ++++++ xradar/io/backends/__init__.py | 4 +- xradar/io/backends/common.py | 13 +- xradar/io/backends/furuno.py | 82 +++++++------ xradar/io/backends/odim.py | 8 +- 10 files changed, 348 insertions(+), 47 deletions(-) create mode 100644 examples/notebooks/Furuno.ipynb diff --git a/docs/history.md b/docs/history.md index ffea69d..e79802f 100644 --- a/docs/history.md +++ b/docs/history.md @@ -2,6 +2,10 @@ ## 0.7.0 (2022-09-21) +* add Furuno SCN/SCNX importer ({pull}`30`) by [@kmuehlbauer](https://github.com/kmuehlbauer) + +## 0.7.0 (2022-09-21) + * Add zenodo badges to README.md ({pull}`22`) by [@mgrover1](https://github.com/mgrover1) * Fix version on RTD ({pull}`23`) by [@kmuehlbauer](https://github.com/kmuehlbauer) * Add minimal documentation for Datamodel ({pull}`24`) by [@kmuehlbauer](https://github.com/kmuehlbauer) diff --git a/docs/importers.md b/docs/importers.md index 0833312..18b16c7 100644 --- a/docs/importers.md +++ b/docs/importers.md @@ -41,3 +41,21 @@ more functions are applied on that {py:class}`xarray:xarray.Dataset`. With {class}`xradar.io.backends.odim.open_odim_datatree` all groups (eg. ``datasetN``) are extracted. From that the ``root`` group is processed. Everything is finally added as ParentNodes and ChildNodes to a {py:class}`datatree:datatree.Datatree`. + + +## Furuno SCN and SCNX + +### FurunoBackendEntrypoint + +The xarray backend {class}`xradar.io.backends.furuno.FurunoBackendEntrypoint` +opens the file with {class}`xradar.io.backends.furuno.FurunoStore`. +Furuno SCN and SCNX data files contain only one one sweep group, so the +group-keyword isn't used. Several private helper functions are used to +conveniently access data and metadata. Finally, the xarray machinery returns +a {py:class}`xarray:xarray.Dataset` with the sweep group. + +### open_furuno_datatree + +With {class}`xradar.io.backends.furuno.open_furuno_datatree` the single group +is 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 1928ded..9347978 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -27,4 +27,5 @@ notebooks/CfRadial1 notebooks/CfRadial1_full notebooks/ODIM_H5 notebooks/ODIM_H5_full +notebooks/Furuno ``` diff --git a/examples/notebooks/Furuno.ipynb b/examples/notebooks/Furuno.ipynb new file mode 100644 index 0000000..feb87a2 --- /dev/null +++ b/examples/notebooks/Furuno.ipynb @@ -0,0 +1,209 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "59447ad6-ac47-494e-b696-4335b36b205b", + "metadata": { + "tags": [] + }, + "source": [ + "# Furuno - Simple" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f96b5d8-2b96-4fd7-b8ba-166c34a8dcd2", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import xarray as xr\n", + "import xradar as xd\n", + "from urllib.parse import urljoin\n", + "from urllib.request import urlretrieve" + ] + }, + { + "cell_type": "markdown", + "id": "33d50be4-dfe5-4d99-a936-67a9a76bac94", + "metadata": {}, + "source": [ + "## Download\n", + "\n", + "Fetching Furuno radar data file from wradlib-data repository." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3c6d408-5ab2-43c3-afd1-b3a703ef3b24", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch_furuno_file():\n", + " fname = \"furuno_data.scn.gz\"\n", + " if not os.path.exists(fname):\n", + " base_url = \"https://raw.githubusercontent.com/wradlib/wradlib-data/main/furuno/\"\n", + " #filename = \"2006_20220324_000000_000.scnx.gz\"\n", + " filename = \"0080_20210730_160000_01_02.scn.gz\"\n", + " url = urljoin(base_url, filename)\n", + " urlretrieve(url, filename=fname)\n", + " return fname\n", + "\n", + "\n", + "filename = fetch_furuno_file()" + ] + }, + { + "cell_type": "markdown", + "id": "b987dcfd-5105-4483-932e-71b8002e5f09", + "metadata": {}, + "source": [ + "## xr.open_dataset\n", + "\n", + "Making use of the xarray `odim` backend. We also need to provide the group." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7675b518-18e4-4ea6-b101-f1bccf603902", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(\n", + " filename, engine=\"furuno\", backend_kwargs=dict(reindex_angle=False)\n", + ")\n", + "display(ds)" + ] + }, + { + "cell_type": "markdown", + "id": "01ec8c90-2da8-46ae-a0b5-0e1792a79bbe", + "metadata": {}, + "source": [ + "### Plot Time vs. Azimuth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24dfb5f8-de27-46f0-9246-722e5955a345", + "metadata": {}, + "outputs": [], + "source": [ + "ds.azimuth.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "222ca704-d6e1-49e6-84a0-de55df9fdf61", + "metadata": {}, + "source": [ + "### Plot Range vs. Time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34429720-e689-4786-99e3-5af9742d19ad", + "metadata": {}, + "outputs": [], + "source": [ + "ds.DBZH.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "d904cd09-8590-42e2-8dce-41d3949d313c", + "metadata": {}, + "source": [ + "### Plot Range vs. Azimuth\n", + "\n", + "We need to sort by azimuth and specify the y-coordinate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6479a374-25ab-42be-b53e-82849b6faffc", + "metadata": {}, + "outputs": [], + "source": [ + "ds.DBZH.sortby(\"azimuth\").plot(y=\"azimuth\")" + ] + }, + { + "cell_type": "markdown", + "id": "5deafdf1-9224-459c-8b10-7c501ae13234", + "metadata": {}, + "source": [ + "## open_odim_datatree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a65f0b48-7197-4f79-bc26-2045cfc59a4a", + "metadata": {}, + "outputs": [], + "source": [ + "dtree = xd.io.open_furuno_datatree(filename)\n", + "display(dtree)" + ] + }, + { + "cell_type": "markdown", + "id": "638bc6c0-3293-4661-a5d5-e9aaad14ffe9", + "metadata": {}, + "source": [ + "### Plot Sweep Range vs. Time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9db8500-7072-451b-84a1-f36767110e16", + "metadata": {}, + "outputs": [], + "source": [ + "dtree[\"sweep_0\"].ds.DBZH.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "1e44c0fd-02bc-4d2d-ac91-772bdcebe04b", + "metadata": {}, + "source": [ + "### Plot Sweep Range vs. Azimuth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7aab6b3-aeb0-4ed1-8397-8b505e63464a", + "metadata": {}, + "outputs": [], + "source": [ + "dtree[\"sweep_0\"].ds.DBZH.sortby(\"azimuth\").plot(y=\"azimuth\")" + ] + } + ], + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/conftest.py b/tests/conftest.py index 8930bf6..e93f393 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -28,3 +28,25 @@ def odim_file(tmp_path_factory): fname = os.path.join(fn, "odim_data.h5") urlretrieve(url, filename=fname) return fname + + +@pytest.fixture(scope="session") +def furuno_scn_file(tmp_path_factory): + base_url = "https://raw.githubusercontent.com/wradlib/wradlib-data/main/furuno/" + filename = "0080_20210730_160000_01_02.scn.gz" + url = urljoin(base_url, filename) + fn = tmp_path_factory.mktemp("furuno_data") + fname = os.path.join(fn, "furuno_data.scn.gz") + urlretrieve(url, filename=fname) + return fname + + +@pytest.fixture(scope="session") +def furuno_scnx_file(tmp_path_factory): + base_url = "https://raw.githubusercontent.com/wradlib/wradlib-data/main/furuno/" + filename = "2006_20220324_000000_000.scnx.gz" + url = urljoin(base_url, filename) + fn = tmp_path_factory.mktemp("furuno_data") + fname = os.path.join(fn, "furuno_data.scnx.gz") + urlretrieve(url, filename=fname) + return fname diff --git a/tests/test_io.py b/tests/test_io.py index ccd1745..aff0a4e 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -183,3 +183,37 @@ def test_open_odim_dataset(odim_file): backend_kwargs=dict(first_dim="auto"), ) assert dict(ds.dims) == {"azimuth": 360, "range": 280} + + +def test_open_furuno_scn_dataset(furuno_scn_file): + # open sweep group + ds = xr.open_dataset(furuno_scn_file, engine="furuno") + assert dict(ds.dims) == {"time": 1385, "range": 602} + assert set(ds.data_vars) & ( + sweep_dataset_vars | non_standard_sweep_dataset_vars + ) == {"KDP", "VRADH", "ZDR", "DBZH", "WRADH", "RHOHV", "PHIDP"} + + # open sweep group, auto + ds = xr.open_dataset( + furuno_scn_file, + engine="furuno", + backend_kwargs=dict(first_dim="auto"), + ) + assert dict(ds.dims) == {"azimuth": 1385, "range": 602} + + +def test_open_furuno_scnx_dataset(furuno_scnx_file): + # open sweep group + ds = xr.open_dataset(furuno_scnx_file, engine="furuno") + assert dict(ds.dims) == {"time": 720, "range": 936} + assert set(ds.data_vars) & ( + sweep_dataset_vars | non_standard_sweep_dataset_vars + ) == {"KDP", "VRADH", "ZDR", "DBZH", "WRADH", "RHOHV", "PHIDP"} + + # open sweep group, auto + ds = xr.open_dataset( + furuno_scnx_file, + engine="furuno", + backend_kwargs=dict(first_dim="auto"), + ) + assert dict(ds.dims) == {"azimuth": 720, "range": 936} diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 55c73b3..6953c5f 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -10,14 +10,14 @@ :maxdepth: 4 .. automodule:: xradar.io.backends.cfradial1 -.. automodule:: xradar.io.backends.furuno .. automodule:: xradar.io.backends.odim +.. automodule:: xradar.io.backends.furuno """ from .cfradial1 import * # noqa -from .odim import * # noqa from .furuno import * # noqa +from .odim 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 106ca4d..a3db177 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -12,13 +12,13 @@ """ +import struct +from collections import OrderedDict + import numpy as np import xarray as xr from datatree import DataTree -import struct -from collections import OrderedDict - def _maybe_decode(attr): try: @@ -178,7 +178,7 @@ def _assign_root(sweeps): # assign root attributes attrs = {} - attrs["Conventions"] = sweeps[0].attrs["Conventions"] + attrs["Conventions"] = sweeps[0].attrs.get("Conventions", "None") attrs.update( { "version": "None", @@ -197,7 +197,6 @@ def _assign_root(sweeps): return root - def _get_fmt_string(dictionary, retsub=False, byte_order="<"): """Get Format String from given dictionary. @@ -286,11 +285,9 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): # IRIS Data Types and corresponding python struct format characters # 4.2 Scalar Definitions, Page 23 # https://docs.python.org/3/library/struct.html#format-characters +# also used for Furuno data types SINT2 = {"fmt": "h", "dtype": "int16"} SINT4 = {"fmt": "i", "dtype": "int32"} UINT1 = {"fmt": "B", "dtype": "unit8"} UINT2 = {"fmt": "H", "dtype": "uint16"} UINT4 = {"fmt": "I", "dtype": "unint32"} - - - diff --git a/xradar/io/backends/furuno.py b/xradar/io/backends/furuno.py index 57d78ff..2ecda95 100644 --- a/xradar/io/backends/furuno.py +++ b/xradar/io/backends/furuno.py @@ -33,7 +33,7 @@ __all__ = [ "FurunoBackendEntrypoint", - #"open_furuno_datatree", + "open_furuno_datatree", ] __doc__ = __doc__.format("\n ".join(__all__)) @@ -45,26 +45,15 @@ from collections import OrderedDict import lat_lon_parser - 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.backends.common import ( - AbstractDataStore, - BackendArray, - BackendEntrypoint, -) +from xarray.core.utils import FrozenDict from xarray.core.variable import Variable -from xarray.core.utils import Frozen, FrozenDict, close_on_error, is_remote_uri -from xarray.backends.store import StoreBackendEntrypoint - -from .common import (SINT2, - SINT4, - UINT1, - UINT2, - UINT4, -_calculate_angle_res, _get_fmt_string, _unpack_dictionary, _attach_sweep_groups, _fix_angle, _maybe_decode, _reindex_angle) from ...model import ( get_altitude_attrs, @@ -77,6 +66,20 @@ moment_attrs, sweep_vars_mapping, ) +from .common import ( + SINT2, + SINT4, + UINT1, + UINT2, + UINT4, + _assign_root, + _attach_sweep_groups, + _calculate_angle_res, + _get_fmt_string, + _reindex_angle, + _unpack_dictionary, +) + def decode_time(data, time_struct=None): """Decode `YMDS_TIME` into datetime object.""" @@ -583,11 +586,7 @@ def open_store_variable(self, name, var): mapping = sweep_vars_mapping.get(name, {}) attrs = {key: mapping[key] for key in moment_attrs if key in mapping} if name in ["azimuth", "elevation"]: - attrs = ( - get_azimuth_attrs() - if name == "azimuth" - else get_elevation_attrs() - ) + attrs = get_azimuth_attrs() if name == "azimuth" else get_elevation_attrs() attrs["add_offset"] = add_offset attrs["scale_factor"] = scale_factor dims = (dim,) @@ -623,15 +622,12 @@ def open_store_coordinates(self): dtype="float32", ) range_attrs = get_range_attrs(rng) - #range_attrs["meters_to_center_of_first_gate"] = start_range + range_step / 2 - #range_attrs["meters_between_gates"] = range_step rng = Variable(("range",), rng, range_attrs) # making-up ray times time = self.ds.header["scan_start_time"] stop_time = self.ds.header.get("scan_stop_time", time) num_rays = self.ds.header["number_sweep_direction_data"] - total_seconds = (time - dt.datetime(1970, 1, 1)).total_seconds() # if no stop_time is available, get time from rotation speed if time == stop_time: @@ -650,30 +646,31 @@ def open_store_coordinates(self): diff = np.diff(raytimes) / 2.0 rtime = raytimes[:-1] + diff - # rtime_attrs = { - # "units": f"seconds since {time.isoformat()}Z", - # "standard_name": "time", - # } rtime_attrs = get_time_attrs(f"{time.isoformat()}Z") encoding = {} - #time_attrs = get_time_attrs("1970-01-01T00:00:00Z") rng = Variable(("range",), rng, range_attrs) time = Variable((dim,), rtime, rtime_attrs, encoding) - #time = Variable((), total_seconds, time_attrs, encoding) # get coordinates from Furuno File sweep_mode = "azimuth_surveillance" if dim == "azimuth" else "rhi" + sweep_number = 0 + prt_mode = "not_set" + follow_mode = "not_set" + lon, lat, alt = self.ds.site_coords coords = { "range": rng, "time": time, - #"rtime": rtime, + "sweep_mode": Variable((), sweep_mode), + "sweep_number": Variable((), sweep_number), + "prt_mode": Variable((), prt_mode), + "follow_mode": Variable((), follow_mode), + "fixed_angle": Variable((), self.ds.fixed_angle), "longitude": Variable((), lon, get_longitude_attrs()), "latitude": Variable((), lat, get_latitude_attrs()), "altitude": Variable((), alt, get_altitude_attrs()), - "sweep_mode": Variable((), sweep_mode), } return coords @@ -690,8 +687,8 @@ def get_variables(self): ) def get_attrs(self): - attributes = {"fixed_angle": float(self.ds.fixed_angle)} - return FrozenDict(attributes) + # attributes = {"fixed_angle": float(self.ds.fixed_angle)} + return FrozenDict() class FurunoBackendEntrypoint(BackendEntrypoint): @@ -751,6 +748,20 @@ def open_dataset( ds = ds.swap_dims({dim0: "time"}) ds = ds.sortby("time") + # 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}) + + # assign geo-coords + ds = ds.assign_coords( + { + "latitude": ds.latitude, + "longitude": ds.longitude, + "altitude": ds.altitude, + } + ) + return ds @@ -793,10 +804,9 @@ def open_furuno_datatree(filename_or_obj, **kwargs): ds = [xr.open_dataset(filename_or_obj, engine="furuno", **kwargs)] - ds.insert(0, xr.open_dataset(filename_or_obj, group="/")) + 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/io/backends/odim.py b/xradar/io/backends/odim.py index 08a61a6..bd01f64 100644 --- a/xradar/io/backends/odim.py +++ b/xradar/io/backends/odim.py @@ -65,7 +65,13 @@ sweep_vars_mapping, ) from ...util import has_import -from .common import _attach_sweep_groups, _fix_angle, _maybe_decode, _reindex_angle, _assign_root +from .common import ( + _assign_root, + _attach_sweep_groups, + _fix_angle, + _maybe_decode, + _reindex_angle, +) HDF5_LOCK = SerializableLock()