Cf/Radial1 writer (#126)

AUTHORS.md

@@ -9,3 +9,4 @@
 ## Contributors
 
 * Edouard Goudenhoofdt <edouard.goudenhoofdt@meteo.be>
+* Hamid Ali Syed <hamidsyed37@gmail.com>

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

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.

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)

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,

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<br>\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
+}

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
+            )

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("_")]