notebooks/multisensor/wradlib_adjust_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden",
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Adjusting radar-base rainfall estimates by rain gauge observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Background\n",
    "\n",
    "There are various ways to correct specific errors and artifacts in radar-based quantitative precipitation estimates (*radar QPE*). Alternatively, you might want to correct your radar QPE regardless of the error source - by using ground truth, or, more specifically, rain gauge observations. Basically, you define the error of your radar QPE at a rain gauge location by the discrepancy between rain gauge observation (considered as \"the truth\") and radar QPE at that very location. Whether you consider this \"discrepancy\" as an *additive* or *multiplicative* error is somehow arbitrary - typically, it's a *mix* of both. If you quantify this error at various locations (i.e. rain gauges), you can go ahead and construct correction fields for your radar QPE. You might compute a single correction factor for your entire radar domain (which would e.g. make sense in case of hardware miscalibration), or you might want to compute a spatially variable correction field. This typically implies to interpolate the error in space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$\\omega radlib$ provides different error models and different spatial interpolation methods to address the adjustment problem. For details, please refer to $\\omega radlib's$ [library reference](http://wradlib.org/wradlib-docs/latest/adjust.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import wradlib.adjust as adjust\n",
    "import wradlib.verify as verify\n",
    "import wradlib.util as util\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example for the 1-dimensional case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Looking at the 1-D (instead of 2-D) case is more illustrative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Create synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First, we **create synthetic data**: \n",
    "- true rainfall, \n",
    "- point observations of the truth, \n",
    "- radar observations of the truth.\n",
    "\n",
    "The latter is disturbed by some kind of error, e.g. a combination between systemtic and random error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# gage and radar coordinates\n",
    "obs_coords = np.array([5, 10, 15, 20, 30, 45, 65, 70, 77, 90])\n",
    "radar_coords = np.arange(0, 101)\n",
    "\n",
    "# true rainfall\n",
    "truth = np.abs(1.5 + np.sin(0.075 * radar_coords)) + np.random.uniform(\n",
    "    -0.1, 0.1, len(radar_coords))\n",
    "\n",
    "# radar error\n",
    "erroradd = 0.7 * np.sin(0.2 * radar_coords + 10.)\n",
    "errormult = 0.75 + 0.015 * radar_coords\n",
    "noise = np.random.uniform(-0.05, 0.05, len(radar_coords))\n",
    "\n",
    "# radar observation\n",
    "radar = errormult * truth + erroradd + noise\n",
    "\n",
    "# gage observations are assumed to be perfect\n",
    "obs = truth[obs_coords]\n",
    "\n",
    "# add a missing value to observations (just for testing)\n",
    "obs[1] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Apply different adjustment methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- additive error, spatially variable (`AdjustAdd`)\n",
    "- multiplicative error, spatially variable (`AdjustMultiply`)\n",
    "- mixed error, spatially variable (`AdjustMixed`)\n",
    "- multiplicative error, spatially uniform (`AdjustMFB`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# number of neighbours to be used\n",
    "nnear_raws = 3\n",
    "\n",
    "# adjust the radar observation by additive model\n",
    "add_adjuster = adjust.AdjustAdd(obs_coords, radar_coords,\n",
    "                                nnear_raws=nnear_raws)\n",
    "add_adjusted = add_adjuster(obs, radar)\n",
    "\n",
    "# adjust the radar observation by multiplicative model\n",
    "mult_adjuster = adjust.AdjustMultiply(obs_coords, radar_coords,\n",
    "                                      nnear_raws=nnear_raws)\n",
    "mult_adjusted = mult_adjuster(obs, radar)\n",
    "\n",
    "# adjust the radar observation by AdjustMixed\n",
    "mixed_adjuster = adjust.AdjustMixed(obs_coords, radar_coords,\n",
    "                                    nnear_raws=nnear_raws)\n",
    "mixed_adjusted = mixed_adjuster(obs, radar)\n",
    "\n",
    "# adjust the radar observation by MFB\n",
    "mfb_adjuster = adjust.AdjustMFB(obs_coords, radar_coords,\n",
    "                                nnear_raws=nnear_raws, \n",
    "                                mfb_args = dict(method=\"median\"))\n",
    "mfb_adjusted = mfb_adjuster(obs, radar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plot adjustment results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Enlarge all label fonts\n",
    "font = {'size'   : 15}\n",
    "pl.rc('font', **font)\n",
    "\n",
    "pl.figure(figsize=(10,5))\n",
    "pl.plot(radar_coords, radar,         'k-', linewidth=2., linestyle=\"dashed\", label=\"Unadjusted radar\", )\n",
    "pl.plot(radar_coords, truth,         'k-', linewidth=2., label=\"True rainfall\", )\n",
    "pl.plot(obs_coords,   obs,           'o',  markersize=10.0, markerfacecolor=\"grey\", label=\"Gage observation\")\n",
    "pl.plot(radar_coords, add_adjusted,  '-',  color=\"red\",   label=\"Additive adjustment\")\n",
    "pl.plot(radar_coords, mult_adjusted, '-',  color=\"green\", label=\"Multiplicative adjustment\")\n",
    "pl.plot(radar_coords, mfb_adjusted,  '-',  color=\"orange\", label=\"Mean Field Bias adjustment\")\n",
    "pl.plot(radar_coords, mixed_adjusted,'-',  color=\"blue\", label=\"Mixed (mult./add.) adjustment\")\n",
    "pl.xlabel(\"Distance (km)\")\n",
    "pl.ylabel(\"Rainfall intensity (mm/h)\")\n",
    "leg = pl.legend(prop={'size': 10})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Verification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We use the `verify` module to compare the errors of different adjustment approaches.\n",
    "\n",
    "*Here, we compare the adjustment to the \"truth\". In practice, we would carry out a cross validation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Verification for this example\n",
    "rawerror = verify.ErrorMetrics(truth, radar)\n",
    "mfberror = verify.ErrorMetrics(truth, mfb_adjusted)\n",
    "adderror = verify.ErrorMetrics(truth, add_adjusted)\n",
    "multerror = verify.ErrorMetrics(truth, mult_adjusted)\n",
    "mixerror = verify.ErrorMetrics(truth, mixed_adjusted)\n",
    "\n",
    "# Verification reports\n",
    "maxval = 4.\n",
    "# Enlarge all label fonts\n",
    "font = {'size'   : 10}\n",
    "pl.rc('font', **font)\n",
    "fig = pl.figure(figsize=(14, 8))\n",
    "ax = fig.add_subplot(231, aspect=1.)\n",
    "rawerror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "ax.text(0.2, 0.9 * maxval, \"Unadjusted radar\")\n",
    "ax = fig.add_subplot(232, aspect=1.)\n",
    "adderror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "ax.text(0.2, 0.9 * maxval, \"Additive adjustment\")\n",
    "ax = fig.add_subplot(233, aspect=1.)\n",
    "multerror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "ax.text(0.2, 0.9 * maxval, \"Multiplicative adjustment\")\n",
    "ax = fig.add_subplot(234, aspect=1.)\n",
    "mixerror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "ax.text(0.2, 0.9 * maxval, \"Mixed (mult./add.) adjustment\")\n",
    "mixerror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "ax = fig.add_subplot(235, aspect=1.)\n",
    "mfberror.report(ax=ax, unit=\"mm\", maxval=maxval)\n",
    "txt = ax.text(0.2, 0.9 * maxval, \"Mean Field Bias adjustment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example for the 2-dimensional case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For the 2-D case, we follow the same approach as before:\n",
    "\n",
    "- create synthetic data: truth, rain gauge observations, radar-based rainfall estimates\n",
    "- apply adjustment methods\n",
    "- verification\n",
    "\n",
    "The way these synthetic data are created is totally arbitrary - it's just to show how the methods are applied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Create 2-D synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# grid axes\n",
    "xgrid = np.arange(0, 10)\n",
    "ygrid = np.arange(20, 30)\n",
    "\n",
    "# number of observations\n",
    "num_obs = 10\n",
    "\n",
    "# create grid\n",
    "gridshape = len(xgrid), len(ygrid)\n",
    "grid_coords = util.gridaspoints(ygrid, xgrid)\n",
    "\n",
    "# Synthetic true rainfall\n",
    "truth = np.abs(10. * np.sin(0.1 * grid_coords).sum(axis=1))\n",
    "\n",
    "# Creating radar data by perturbing truth with multiplicative and\n",
    "# additive error\n",
    "# YOU CAN EXPERIMENT WITH THE ERROR STRUCTURE\n",
    "radar = 0.6 * truth + 1. * np.random.uniform(low=-1., high=1,\n",
    "                                             size=len(truth))\n",
    "radar[radar < 0.] = 0.\n",
    "\n",
    "# indices for creating obs from raw (random placement of gauges)\n",
    "obs_ix = np.random.uniform(low=0, high=len(grid_coords),\n",
    "                           size=num_obs).astype('i4')\n",
    "\n",
    "# creating obs_coordinates\n",
    "obs_coords = grid_coords[obs_ix]\n",
    "\n",
    "# creating gauge observations from truth\n",
    "obs = truth[obs_ix]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Apply different adjustment methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Mean Field Bias Adjustment\n",
    "mfbadjuster = adjust.AdjustMFB(obs_coords, grid_coords)\n",
    "mfbadjusted = mfbadjuster(obs, radar)\n",
    "\n",
    "# Additive Error Model\n",
    "addadjuster = adjust.AdjustAdd(obs_coords, grid_coords)\n",
    "addadjusted = addadjuster(obs, radar)\n",
    "\n",
    "# Multiplicative Error Model\n",
    "multadjuster = adjust.AdjustMultiply(obs_coords, grid_coords)\n",
    "multadjusted = multadjuster(obs, radar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plot 2-D adjustment results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Two helper functions for repeated plotting tasks\n",
    "def scatterplot(x, y, title):\n",
    "    \"\"\"Quick and dirty helper function to produce scatter plots\n",
    "    \"\"\"\n",
    "    pl.scatter(x, y)\n",
    "    pl.plot([0, 1.2 * maxval], [0, 1.2 * maxval], '-', color='grey')\n",
    "    pl.xlabel(\"True rainfall (mm)\")\n",
    "    pl.ylabel(\"Estimated rainfall (mm)\")\n",
    "    pl.xlim(0, maxval + 0.1 * maxval)\n",
    "    pl.ylim(0, maxval + 0.1 * maxval)\n",
    "    pl.title(title)\n",
    "\n",
    "def gridplot(data, title):\n",
    "    \"\"\"Quick and dirty helper function to produce a grid plot\n",
    "    \"\"\"\n",
    "    xplot = np.append(xgrid, xgrid[-1] + 1.) - 0.5\n",
    "    yplot = np.append(ygrid, ygrid[-1] + 1.) - 0.5\n",
    "    grd = ax.pcolormesh(xplot, yplot, data.reshape(gridshape), vmin=0,\n",
    "                        vmax=maxval)\n",
    "    ax.scatter(obs_coords[:, 0], obs_coords[:, 1], c=obs.ravel(),\n",
    "               marker='s', s=50, vmin=0, vmax=maxval)\n",
    "    #pl.colorbar(grd, shrink=0.5)\n",
    "    pl.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Maximum value (used for normalisation of colorscales)\n",
    "maxval = np.max(np.concatenate((truth, radar, obs, addadjusted)).ravel())\n",
    "\n",
    "# open figure\n",
    "fig = pl.figure(figsize=(10, 6))\n",
    "\n",
    "# True rainfall\n",
    "ax = fig.add_subplot(231, aspect='equal')\n",
    "gridplot(truth, 'True rainfall')\n",
    "\n",
    "# Unadjusted radar rainfall\n",
    "ax = fig.add_subplot(232, aspect='equal')\n",
    "gridplot(radar, 'Radar rainfall')\n",
    "\n",
    "# Adjusted radar rainfall (MFB)\n",
    "ax = fig.add_subplot(234, aspect='equal')\n",
    "gridplot(mfbadjusted, 'Adjusted (MFB)')\n",
    "\n",
    "# Adjusted radar rainfall (additive)\n",
    "ax = fig.add_subplot(235, aspect='equal')\n",
    "gridplot(addadjusted, 'Adjusted (Add.)')\n",
    "\n",
    "# Adjusted radar rainfall (multiplicative)\n",
    "ax = fig.add_subplot(236, aspect='equal')\n",
    "gridplot(multadjusted, 'Adjusted (Mult.)')\n",
    "\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Open figure\n",
    "fig = pl.figure(figsize=(6, 6))\n",
    "\n",
    "# Scatter plot radar vs. observations\n",
    "ax = fig.add_subplot(221, aspect='equal')\n",
    "scatterplot(truth, radar, 'Radar vs. Truth (red: Gauges)')\n",
    "pl.plot(obs, radar[obs_ix], linestyle=\"None\", marker=\"o\", color=\"red\")\n",
    "\n",
    "# Adjusted (MFB) vs. radar (for control purposes)\n",
    "ax = fig.add_subplot(222, aspect='equal')\n",
    "scatterplot(truth, mfbadjusted, 'Adjusted (MFB) vs. Truth')\n",
    "\n",
    "# Adjusted (Add) vs. radar (for control purposes)\n",
    "ax = fig.add_subplot(223, aspect='equal')\n",
    "scatterplot(truth, addadjusted, 'Adjusted (Add.) vs. Truth')\n",
    "\n",
    "# Adjusted (Mult.) vs. radar (for control purposes)\n",
    "ax = fig.add_subplot(224, aspect='equal')\n",
    "scatterplot(truth, multadjusted, 'Adjusted (Mult.) vs. Truth')\n",
    "\n",
    "pl.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  },
  "livereveal": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}