From f46e5d16039cae51d3393d762869359631593a27 Mon Sep 17 00:00:00 2001 From: Hamid Ali Syed <35923822+syedhamidali@users.noreply.github.com> Date: Wed, 27 Sep 2023 07:39:40 -0400 Subject: [PATCH] Cf/Radial1 writer (#126) --- AUTHORS.md | 1 + CITATION.cff | 10 +- docs/exporters.md | 8 + docs/history.md | 3 +- examples/notebooks/CfRadial1.ipynb | 2 +- examples/notebooks/CfRadial1_Export.ipynb | 253 ++++++++++++++++++ tests/io/test_cfradial1.py | 28 ++ xradar/io/export/__init__.py | 1 + xradar/io/export/cfradial1.py | 311 ++++++++++++++++++++++ xradar/model.py | 2 +- 10 files changed, 613 insertions(+), 6 deletions(-) create mode 100644 examples/notebooks/CfRadial1_Export.ipynb create mode 100644 tests/io/test_cfradial1.py create mode 100644 xradar/io/export/cfradial1.py diff --git a/AUTHORS.md b/AUTHORS.md index c22a005..b0d3ad4 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -9,3 +9,4 @@ ## Contributors * Edouard Goudenhoofdt +* Hamid Ali Syed diff --git a/CITATION.cff b/CITATION.cff index b25b2ba..ce8a17a 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -3,7 +3,7 @@ cff-version: 1.0.3 message: If you use this software, please cite it using these metadata. # FIXME title as repository name might not be the best name, please make human readable -title: 'openradar/xradar: xradar v0.3.0' +title: 'openradar/xradar: xradar v0.4.0' doi: 10.5281/zenodo.7091737 # FIXME splitting of full names is error prone, please check if given/family name are correct authors: @@ -19,7 +19,11 @@ authors: family-names: Goudenhoofdt affiliation: Royal Meteorological Institute of Belgium orcid: https://orcid.org/0000-0002-6376-1515 -version: 0.3.0 -date-released: 2023-07-11 +- given-names: Hamid Ali + family-names: Syed + affiliation: Purdue University + orcid: https://orcid.org/0000-0002-7188-2544 +version: 0.4.0 +date-released: 2023-09-26 repository-code: https://github.com/openradar/xradar license: MIT diff --git a/docs/exporters.md b/docs/exporters.md index 07b66ff..fbba992 100644 --- a/docs/exporters.md +++ b/docs/exporters.md @@ -6,6 +6,7 @@ Currently xradar can export: - [](#odim_h5) - [](#cfradial2) +- [](#cfradial1) ## ODIM_H5 @@ -20,3 +21,10 @@ can be saved to an ODIM_H5 file (v2.2 at the moment). With {class}`~xradar.io.export.to_cfradial2` an xradar {py:class}`datatree:datatree.DataTree` can be saved to a CfRadial2-like file. + +## CfRadial1 + +### to_cfradial1 + +With {class}`~xradar.io.export.to_cfradial1` an xradar {py:class}`datatree:datatree.DataTree` +can be saved to a CfRadial1-like file. diff --git a/docs/history.md b/docs/history.md index 0561d91..05ad049 100644 --- a/docs/history.md +++ b/docs/history.md @@ -2,7 +2,8 @@ ## development version -* FIX: use datastore._group instead of variable["sweep_number"] ({issue}`121`) by [@aladinor](https://github.com/aladinor) , {pull}`123`) by [@kmuehlbauer](https://github.com/kmuehlbauer) +* ENH: Add cfradial1 exporter ({issue}`124`) by [@syedhamidali](https://github.com/syedhamidali), ({pull}`126`) by [@syedhamidali](https://github.com/syedhamidali), improved by [@kmuehlbauer](https://github.com/kmuehlbauer) +* FIX: use datastore._group instead of variable["sweep_number"] ({issue}`121`) by [@aladinor](https://github.com/aladinor) , ({pull}`123`) by [@kmuehlbauer](https://github.com/kmuehlbauer) * MIN: use "crs_wkt" instead of deprecated "spatial_ref" when adding CRS ({pull}`127`) by [@kmuehlbauer](https://github.com/kmuehlbauer) ## 0.3.0 (2023-07-11) diff --git a/examples/notebooks/CfRadial1.ipynb b/examples/notebooks/CfRadial1.ipynb index 07df193..5572950 100644 --- a/examples/notebooks/CfRadial1.ipynb +++ b/examples/notebooks/CfRadial1.ipynb @@ -268,7 +268,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.11.5" } }, "nbformat": 4, diff --git a/examples/notebooks/CfRadial1_Export.ipynb b/examples/notebooks/CfRadial1_Export.ipynb new file mode 100644 index 0000000..3bf7eea --- /dev/null +++ b/examples/notebooks/CfRadial1_Export.ipynb @@ -0,0 +1,253 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# CfRadial1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import tempfile\n", + "import cmweather\n", + "import numpy as np\n", + "import xarray as xr\n", + "import xradar as xd\n", + "import datatree as xt\n", + "from open_radar_data import DATASETS" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Download\n", + "\n", + "Fetching CfRadial1 radar data file from [open-radar-data](https://github.com/openradar/open-radar-data) repository." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filename = DATASETS.fetch(\"cfrad.20080604_002217_000_SPOL_v36_SUR.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radar = xd.io.open_cfradial1_datatree(filename, first_dim=\"auto\")\n", + "display(radar)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot Range vs. Time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radar.sweep_0.DBZ.plot(cmap=\"ChaseSpectral\", vmin=-10, vmax=70)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Georeference" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radar = radar.xradar.georeference()\n", + "display(radar)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot PPI" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radar[\"sweep_0\"][\"DBZ\"].plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=70)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Filter\n", + "\n", + "Apply basic reflectivity filter. This is just a demonstration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def ref_filter(dtree, sweep=\"sweep_0\", field=\"DBZ\"):\n", + " ds = dtree[sweep].where((dtree[sweep][field] >= -10) & (dtree[sweep][field] <= 70))\n", + " red_patch = ds.where(\n", + " (\n", + " (ds[field] >= ds[field].max().values - 0.5)\n", + " & (ds[field] <= ds[field].max().values + 0.5)\n", + " ),\n", + " drop=True,\n", + " )\n", + " rmin, rmax = int(red_patch.range.min().values - 150), int(\n", + " red_patch.range.max().values + 150\n", + " )\n", + " out_of_range_mask = (ds.range < rmin) | (ds.range > rmax)\n", + " ds[field] = ds[field].where(out_of_range_mask)\n", + " # Interpolate missing values using the slinear method along the 'range' dimension\n", + " ds[field] = ds[field].interpolate_na(dim=\"range\", method=\"slinear\")\n", + " dtree[sweep][f\"corr_{field}\"] = ds[field].copy()\n", + " return dtree[sweep]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "swp0 = ref_filter(radar, sweep=\"sweep_0\", field=\"DBZ\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "swp0.corr_DBZ.plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=70)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Filter full volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize an empty DataTree\n", + "result_tree = xt.DataTree()\n", + "\n", + "for sweep in radar.sweep_group_name.values:\n", + " corrected_data = ref_filter(radar, sweep, field=\"DBZ\")\n", + "\n", + " # Convert the xarray Dataset to a DataTree and add it to the result_tree\n", + " data_tree = xt.DataTree.from_dict(corrected_data.to_dict())\n", + "\n", + " # Copy the contents of data_tree into result_tree\n", + " for key, value in data_tree.items():\n", + " result_tree[key] = value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radar.sweep_6.corr_DBZ.plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=70)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Export\n", + "\n", + "Export to CfRadial1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "xd.io.to_cfradial1(dtree=radar, filename=\"cfradial1_qced.nc\", calibs=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "?xd.io.to_cfradial1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Note \n", + "\n", + "If `filename` is `None` in the `xd.io.to_cfradial1` function, it will automatically generate a
\n", + "filename using the instrument name and the first available timestamp from the data.\n" + ] + } + ], + "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": 2 +} diff --git a/tests/io/test_cfradial1.py b/tests/io/test_cfradial1.py new file mode 100644 index 0000000..2d3dfe6 --- /dev/null +++ b/tests/io/test_cfradial1.py @@ -0,0 +1,28 @@ +import tempfile + +import xarray as xr +from open_radar_data import DATASETS + +import xradar as xd + + +def test_compare_sweeps(): + # Fetch the radar data file + filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc") + + # Open the data tree + dtree = xd.io.open_cfradial1_datatree(filename) + + # Create a temporary file to store the modified data tree + with tempfile.NamedTemporaryFile(mode="w+b") as temp_file: + # Save the modified data tree to the temporary file + xd.io.to_cfradial1(dtree.copy(), temp_file.name, calibs=True) + + # Open the modified data tree + dtree1 = xd.io.open_cfradial1_datatree(temp_file.name) + + # Compare the values of the DataArrays for all sweeps + for sweep_num in range(9): # there are 9 sweeps in this file + xr.testing.assert_equal( + dtree[f"sweep_{sweep_num}"].ds, dtree1[f"sweep_{sweep_num}"].ds + ) diff --git a/xradar/io/export/__init__.py b/xradar/io/export/__init__.py index fb3f743..fc9c364 100644 --- a/xradar/io/export/__init__.py +++ b/xradar/io/export/__init__.py @@ -15,6 +15,7 @@ """ from .cfradial2 import * # noqa +from .cfradial1 import * # noqa from .odim import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/export/cfradial1.py b/xradar/io/export/cfradial1.py new file mode 100644 index 0000000..03415f5 --- /dev/null +++ b/xradar/io/export/cfradial1.py @@ -0,0 +1,311 @@ +#!/usr/bin/env python +# Copyright (c) 2023, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +""" + +CfRadial1 output +================ + +This sub-module contains the writer for export of CfRadial1-based radar +data. + +Example:: + + import xradar as xd + xd.io.to_cfradial1(dtree, filename, calibs=True) + +.. autosummary:: + :nosignatures: + :toctree: generated/ + + {} + +""" + +__all__ = [ + "to_cfradial1", +] + +import numpy as np +import xarray as xr + + +def _calib_mapper(calib_params): + """ + Map calibration parameters to a new dataset format. + + Parameters + ---------- + calib_params: xarray.Dataset + Calibration parameters dataset. + + Returns + ------- + xarray.Dataset + New dataset with mapped calibration parameters. + """ + new_data_vars = {} + for var in calib_params.data_vars: + data_array = calib_params[var] + new_data_vars["r_calib_" + var] = xr.DataArray( + data=data_array.data[np.newaxis, ...], + dims=["r_calib"] + list(data_array.dims), + coords={"r_calib": [0]}, + attrs=data_array.attrs, + ) + radar_calib_renamed = xr.Dataset(new_data_vars) + dummy_ds = radar_calib_renamed.rename_vars({"r_calib": "fake_coord"}) + del dummy_ds["fake_coord"] + return dummy_ds + + +def _main_info_mapper(dtree): + """ + Map main radar information from a radar datatree dataset. + + Parameters + ---------- + dtree: datatree.DataTree + Radar datatree. + + Returns + ------- + xarray.Dataset + Dataset containing the mapped radar information. + """ + dataset = ( + dtree.root.to_dataset() + .drop_vars("sweep_group_name", errors="ignore") + .rename({"sweep_fixed_angle": "fixed_angle"}) + ) + return dataset + + +def _variable_mapper(dtree, dim0=None): + """ + Map radar variables for different sweep groups. + + Parameters + ---------- + dtree: datatree.DataTree + Radar datatree. + dim0: str + Either `azimuth` or `elevation` + + Returns + ------- + xarray.Dataset + Dataset containing mapped radar variables. + """ + + sweep_info = _sweep_info_mapper(dtree) + vol_info = _main_info_mapper(dtree).drop("fixed_angle") + sweep_datasets = [] + for grp in dtree.groups: + if "sweep" in grp: + data = dtree[grp] + + # handling first dimension + if dim0 is None: + dim0 = ( + "elevation" + if data.sweep_mode.load().astype(str) == "rhi" + else "azimuth" + ) + if dim0 not in data.dims: + dim0 = "time" + assert dim0 in data.dims + + # swap dims, if needed + if dim0 != "time" and dim0 in data.dims: + data = data.swap_dims({dim0: "time"}) + + # sort in any case + data = data.sortby("time") + + data = data.drop_vars(["x", "y", "z"], errors="ignore") + + # Convert to a dataset and append to the list + sweep_datasets.append(data.to_dataset()) + + result_dataset = xr.concat( + sweep_datasets, + dim="time", + compat="no_conflicts", + join="right", + combine_attrs="drop_conflicts", + ) + + # Check if specific variables exist before dropping them + drop_variables = [ + "sweep_fixed_angle", + "sweep_number", + "sweep_mode", + "prt_mode", + "follow_mode", + ] + result_dataset = result_dataset.drop_vars(drop_variables, errors="ignore") + + drop_coords = ["latitude", "longitude", "altitude", "spatial_ref", "crs_wkt"] + result_dataset = result_dataset.drop_vars(drop_coords, errors="ignore") + + result_dataset = result_dataset.update(sweep_info) + sweep_indices = calculate_sweep_indices(dtree, result_dataset) + result_dataset = result_dataset.update(sweep_indices) + result_dataset = result_dataset.reset_coords(["elevation", "azimuth"]) + result_dataset = result_dataset.update(vol_info) + return result_dataset + + +def _sweep_info_mapper(dtree): + """ + Extract specified sweep information variables from a radar datatree + + Parameters + ---------- + dtree: datatree.DataTree + Radar datatree. + + Returns + ------- + xarray.Dataset + Dataset containing the specified sweep information variables. + """ + dataset = xr.Dataset() + + sweep_vars = [ + "sweep_number", + "sweep_mode", + "polarization_mode", + "prt_mode", + "follow_mode", + "sweep_fixed_angle", + "sweep_start_ray_index", + "sweep_end_ray_index", + ] + + for var_name in sweep_vars: + var_data_list = [ + np.unique(dtree[s][var_name].values[np.newaxis, ...]) + if var_name in dtree[s] + else np.array([np.nan]) + for s in dtree.groups + if "sweep" in s + ] + + var_attrs_list = [ + dtree[s][var_name].attrs if var_name in dtree[s] else {} + for s in dtree.groups + if "sweep" in s + ] + + if not var_data_list: + var_data = np.array([np.nan]) + else: + var_data = np.concatenate(var_data_list) + + var_attrs = {} + for attrs in var_attrs_list: + var_attrs.update(attrs) + + var_data_array = xr.DataArray(var_data, dims=("sweep",), attrs=var_attrs) + dataset[var_name] = var_data_array + + dataset = dataset.rename({"sweep_fixed_angle": "fixed_angle"}) + + return dataset + + +def calculate_sweep_indices(dtree, dataset=None): + """ + Calculate sweep start and end ray indices for elevation + values in a radar dataset. + + Parameters + ---------- + dtree: datatree.DataTree + Radar datatree containing elevation values for different sweep groups. + dataset: xarray.Dataset, optional + An optional dataset to which the calculated indices will be added. + If None, a new dataset will be created. + + Returns: + xarray.Dataset + Dataset with sweep start and end ray indices. + """ + if dataset is None: + dataset = xr.Dataset() + + sweep_start_ray_index = [] + sweep_end_ray_index = [] + + cumulative_size = 0 + + try: + for group_name in dtree.groups: + if "sweep" in group_name: + ele_size = dtree[group_name].elevation.size + sweep_start_ray_index.append(cumulative_size) + sweep_end_ray_index.append(cumulative_size + ele_size - 1) + cumulative_size += ele_size + + except KeyError as e: + print( + f"Error: The sweep group '{e.args[0]}' was not found in radar datatree. Skipping..." + ) + + dataset["sweep_start_ray_index"] = xr.DataArray( + sweep_start_ray_index, + dims=("sweep",), + attrs={"standard_name": "index_of_first_ray_in_sweep"}, + ) + + dataset["sweep_end_ray_index"] = xr.DataArray( + sweep_end_ray_index, + dims=("sweep",), + attrs={"standard_name": "index_of_last_ray_in_sweep"}, + ) + + return dataset + + +def to_cfradial1(dtree=None, filename=None, calibs=True): + """ + Convert a radar dtreeume dataset to the CFRadial1 format + and save it to a file. + + Parameters + ---------- + dtree: datatree.DataTree + Radar datatree object. + filename: str, optional + The name of the output netCDF file. + calibs: Bool, optional + calibration parameters + """ + dataset = _variable_mapper(dtree) + + # Check if radar_parameters, radar_calibration, and + # georeferencing_correction exist in dtree + if calibs: + if "radar_calibration" in dtree: + calib_params = dtree["radar_calibration"].to_dataset() + calibs = _calib_mapper(calib_params) + dataset.update(calibs) + + if "radar_parameters" in dtree: + radar_params = dtree["radar_parameters"].to_dataset() + dataset.update(radar_params) + + if "georeferencing_correction" in dtree: + radar_georef = dtree["georeferencing_correction"].to_dataset() + dataset.update(radar_georef) + + dataset.attrs = dtree.attrs + + if filename is None: + time = str(dataset.time[0].dt.strftime("%Y%m%d_%H%M%S").values) + filename = f"cfrad1_{dataset.instrument_name}_{time}.nc" + + dataset.to_netcdf(filename, format="netcdf4") diff --git a/xradar/model.py b/xradar/model.py index 86bb0b0..6f0ab3a 100644 --- a/xradar/model.py +++ b/xradar/model.py @@ -995,7 +995,7 @@ def determine_cfradial2_sweep_variables(obj, optional, dim0): def conform_cfradial2_sweep_group(obj, optional, dim0=None): if dim0 is None: # handling first dimension - dim0 = "elevation" if obj.sweep_mode.load() == "rhi" else "azimuth" + dim0 = "elevation" if obj.sweep_mode.load().astype(str) == "rhi" else "azimuth" if dim0 not in obj.dims: dim0 = "time" assert dim0 in obj.dims