{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b887302e-3f3e-4d95-acaf-4c1be3d945ac",
   "metadata": {},
   "source": [
    "# Angle Reindexing\n",
    "\n",
    "As a legacy from wradlib we have complex code for angle reindexing in xradar's codebase.\n",
    "\n",
    "## High precision angle coordinates\n",
    "\n",
    "As radar angle coordinates (`azimuth` or `elevation`) are measured constantly by different techniques of detection of the antenna pointing direction the values in the data files are mostly floating point numbers. In many cases these floating point numbers are not rounded to a certain decimal but keep the full possible range of the used dtype.\n",
    "\n",
    "Problems of that:\n",
    "\n",
    "- 1D angle coordinate arrays yield no equidistant vector.\n",
    "- 1D angle coordinate arrays are not equivalent for different timesteps but same scan setup\n",
    "\n",
    "## Missing rays, duplicate or additional rays\n",
    "\n",
    "Sometimes rays (even sectors) are missing from the dataset, sometimes there are duplicate rays. Another problem with radar data are additional rays, which I call \"antenna hickup\" (two rays measured with within one resolution interval).\n",
    "\n",
    "## What is angle reindexing?\n",
    "\n",
    "Angle reindexing takes care of these problems by trying to determine the wanted layout from the radar metadata and the angle coordinates. With that newly created angle coordinate xarray machinery is used to reindex the radar moment data to that by nearest neighbor lookup (up to a tolerance). Missing rays will be filled with NaN.\n",
    "\n",
    "## Why should it be used?\n",
    "\n",
    "For most operations this is not a real problem. It will turn into a problem, if you want to stack your xarray.Dataset radar data on a third dimension (eg. `time`, by using `open_mfdataset`). Then all coordinates need alignment to keep things simple and manageable (eg. `azimuth=[0.5, 1.5, 2.5,..., 359.5]`)\n",
    "\n",
    "## How should we treat it?\n",
    "\n",
    "Currently the reindexing code relies on some internals which make things a bit hard to maintain. My suggestion would be to disentangle the reindexing code from the internals but feed the needed values as parameters. Then every reader can call this per activated `reindex_angle` kwarg."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f525ad9-98a2-4fe4-8a24-83ab67541918",
   "metadata": {},
   "source": [
    "## Angle Reindexing Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffbca852-7358-44c9-b382-99362cbcdf58",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import xradar as xd\n",
    "from open_radar_data import DATASETS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bb792db-463f-49a4-818a-e24f614c4c4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = DATASETS.fetch(\"DWD-Vol-2_99999_20180601054047_00.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dddf7f00-1080-4081-b059-66fb06fa01df",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fix_angle(ds):\n",
    "    angle_dict = xd.util.extract_angle_parameters(ds)\n",
    "    display(angle_dict)\n",
    "    start_ang = angle_dict[\"start_angle\"]\n",
    "    stop_ang = angle_dict[\"stop_angle\"]\n",
    "    angle_res = angle_dict[\"angle_res\"]\n",
    "    direction = angle_dict[\"direction\"]\n",
    "\n",
    "    # first find exact duplicates and remove\n",
    "    ds = xd.util.remove_duplicate_rays(ds)\n",
    "\n",
    "    # second reindex according to retrieved parameters\n",
    "    ds = xd.util.reindex_angle(\n",
    "        ds, start_ang, stop_ang, angle_res, direction, method=\"nearest\"\n",
    "    )\n",
    "\n",
    "    return ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6662525-ecc8-44c0-96d3-51ff43d07d3d",
   "metadata": {},
   "source": [
    "### Read example data with one additional ray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9670d6d0-227a-4dbc-a91a-81206ffdffb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds0 = xr.open_dataset(filename, group=\"sweep_7\", engine=\"gamic\", first_dim=\"auto\")\n",
    "display(ds0.load())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a8eda80-da3e-4158-8e5b-f4b6741f9a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds0.DBTH.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "688cea77-55db-4b90-905a-d826abf0d401",
   "metadata": {},
   "source": [
    "### Prepare sweep with several sections removed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce60e202-da51-41f2-9950-8e849f8db03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_in = xr.concat(\n",
    "    [\n",
    "        ds0.isel(azimuth=slice(0, 100)),\n",
    "        ds0.isel(azimuth=slice(150, 200)),\n",
    "        ds0.isel(azimuth=slice(243, 300)),\n",
    "        ds0.isel(azimuth=slice(330, 361)),\n",
    "    ],\n",
    "    \"azimuth\",\n",
    "    data_vars=\"minimal\",\n",
    ")\n",
    "display(ds_in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6814e7d-2896-4b08-8ade-4db94f571f24",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_in.DBTH.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bcea677-db47-43bd-aaa8-8ba4199faa0a",
   "metadata": {},
   "source": [
    "### Reindex angle\n",
    "\n",
    "First output is the extracted angle/time dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c46977b6-864d-4d75-8dba-c64f9599f2a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out = fix_angle(ds_in)\n",
    "display(ds_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a28522ee-b978-4a7c-a9b4-9696d2e6e276",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out.time.plot(marker=\".\")\n",
    "plt.gca().grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9679aa64-3cc9-411e-bce9-8a3e8c2e9144",
   "metadata": {},
   "source": [
    "We can observe that the dataset is aligned to it's expected number of rays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d161523b-71f1-4ee3-a894-3b2342dd12a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out.DBTH.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4f56ce8-8234-4525-b2f1-a1281dae5951",
   "metadata": {},
   "source": [
    "### Fix timestamps\n",
    "\n",
    "As reindexing instantiates the variables/coordinates added rays with `NaN`/`NaT` we need to take care of the coordinates.\n",
    "The second angle (`elevation` in this case is already treated while reindexing by inserting it's median value, the time coordinate needs special handling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b15dd0ed-b4bf-4655-947d-24282c31a763",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out2 = ds_out.copy(deep=True)\n",
    "ds_out2 = ds_out2.pipe(xd.util.ipol_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dfc6e11-2d7c-4c49-9980-01a3eb34c6f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out2.time.plot(marker=\".\")\n",
    "plt.gca().grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7c83de6-c012-42a6-b860-9c1def5b7f61",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_out2.DBTH.sortby(\"time\").plot(y=\"time\")"
   ]
  }
 ],
 "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
}