notebooks/verification/wradlib_verify_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of the $\\omega radlib$ documentation: http://wradlib.org/wradlib-docs.\n",
    "\n",
    "Copyright (c) 2018, $\\omega radlib$ developers.\n",
    "Distributed under the MIT License. See LICENSE.txt for more info."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flake8: noqa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Routine verification measures for radar-based precipitation estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import wradlib\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract bin values from a polar radar data set at rain gage locations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read polar data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = wradlib.util.get_wradlib_data_file('misc/polar_R_tur.gz')\n",
    "data = np.loadtxt(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define site coordinates (lon/lat) and polar coordinate system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.arange(1,129)\n",
    "az = np.linspace(0,360,361)[0:-1]\n",
    "sitecoords = (9.7839, 48.5861)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make up two rain gauge locations (say we want to work in Gaus Krueger zone 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the projection via epsg-code\n",
    "proj = wradlib.georef.epsg_to_osr(31467)\n",
    "# Coordinates of the rain gages in Gauss-Krueger 3 coordinates\n",
    "x, y = np.array([3557880, 3557890]), np.array([5383379, 5383375])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now extract the radar values at those bins that are closest to our rain gauges\n",
    "\n",
    "For this purppose, we use the PolarNeighbours class from wraldib's verify module. Here, we extract the 9 nearest bins..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polarneighbs = wradlib.verify.PolarNeighbours(r, az, sitecoords, proj, x, y, nnear=9)\n",
    "radar_at_gages = polarneighbs.extract(data)\n",
    "print( \"Radar values at rain gauge #1: %r\" % radar_at_gages[0].tolist())\n",
    "print( \"Radar values at rain gauge #2: %r\" % radar_at_gages[1].tolist())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Retrieve the bin coordinates (all of them or those at the rain gauges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binx, biny = polarneighbs.get_bincoords()\n",
    "binx_nn, biny_nn = polarneighbs.get_bincoords_at_points()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the entire radar domain and zoom into the surrounding of the rain gauge locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,12))\n",
    "ax = fig.add_subplot(121)\n",
    "ax.plot(binx, biny, 'r+')\n",
    "ax.plot(binx_nn, biny_nn, 'b+', markersize=10)\n",
    "ax.plot(x, y, 'bo')\n",
    "ax.axis('tight')\n",
    "ax.set_aspect(\"equal\")\n",
    "pl.title(\"Full view\")\n",
    "ax = fig.add_subplot(122)\n",
    "ax.plot(binx, biny, 'r+')\n",
    "ax.plot(binx_nn, biny_nn, 'b+', markersize=10)\n",
    "ax.plot(x, y, 'bo')\n",
    "pl.xlim(binx_nn.min()-5, binx_nn.max()+5)\n",
    "pl.ylim(biny_nn.min()-7, biny_nn.max()+8)\n",
    "ax.set_aspect(\"equal\")\n",
    "txt = pl.title(\"Zoom into rain gauge locations\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a verification report\n",
    "\n",
    "In this example, we make up a true Kdp profile and verify our reconstructed Kdp."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create synthetic data and reconstruct KDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Synthetic truth\n",
    "dr = 0.5\n",
    "r = np.arange(0, 100, dr)\n",
    "kdp_true = np.sin(0.3*r)\n",
    "kdp_true[kdp_true<0] = 0.\n",
    "phidp_true = np.cumsum(kdp_true)*2*dr\n",
    "# Synthetic observation of PhiDP with a random noise and gaps\n",
    "phidp_raw = phidp_true + np.random.uniform(-2, 2, len(phidp_true))\n",
    "gaps = np.random.uniform(0, len(r), 20).astype(\"int\")\n",
    "phidp_raw[gaps] = np.nan\n",
    "\n",
    "# Reconstruct PhiDP and KDP\n",
    "phidp_re, kdp_re = wradlib.dp.process_raw_phidp_vulpiani(phidp_raw, dr=dr)\n",
    "\n",
    "# Plot results\n",
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(211)\n",
    "pl.plot(kdp_true, \"g-\", label=\"True KDP\")\n",
    "pl.plot(kdp_re, \"r-\", label=\"Reconstructed KDP\")\n",
    "pl.grid()\n",
    "lg = pl.legend()\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "pl.plot(r, phidp_true, \"b--\", label=\"True PhiDP\")\n",
    "pl.plot(r, np.ma.masked_invalid(phidp_raw), \"b-\", label=\"Raw PhiDP\")\n",
    "pl.plot(r, phidp_re, \"g-\", label=\"Reconstructed PhiDP\")\n",
    "pl.grid()\n",
    "lg = pl.legend(loc=\"lower right\")\n",
    "txt = pl.xlabel(\"Range (km)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create verification report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "metrics = wradlib.verify.ErrorMetrics(kdp_true, kdp_re)\n",
    "metrics.pprint() \n",
    "ax = metrics.plot()\n",
    "xlabel = ax.set_xlabel(\"True KDP (deg/km)\")\n",
    "ylabel = ax.set_ylabel(\"Reconstructed KDP (deg/km)\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}