notebooks/attenuation/wradlib_attenuation.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": [
    "# Attenuation correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In this example, we will compare different approaches to **constrain the gate-by-gate retrieval** of path-integrated attenuation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Rainfall-induced attenuation is a major source of underestimation for radar-based precipitation estimation at C-band and X-band. Unconstrained forward gate-by-gate correction is known to be inherently unstable and thus not suited for unsupervised quality control procedures. Ideally, reference measurements (e.g. from microwave links) should be used to constrain gate-by-gate procedures. However, such attenuation references are usually not available. $\\omega radlib$ provides a pragmatic approach to constrain gate-by-gate correction procedures, inspired by the work of [Kraemer et al., 2008](http://wradlib.org/wradlib-docs/latest/zreferences.html#kraemer2008). It turned out that these procedures can effectively reduce the error introduced by attenuation, and, at the same time, minimize instability issues [(Jacobi et al., 2016)](http://wradlib.org/wradlib-docs/latest/zreferences.html#jacobi2015)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "nbsphinx": "hidden",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import wradlib as wrl\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()\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Example Event: June 2nd, 2008, SW-Germany"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's have a look at the situation in South-West Germany on June 2nd, 2008, at 16:55 UTC, as observed by the DWD C-band radar on mount Feldberg.\n",
    "\n",
    "The data can be read by the following lines and then visualized by [wradlib.vis.plot_ppi()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.vis.plot_ppi.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fpath = 'dx/raa00-dx_10908-0806021655-fbg---bin.gz'\n",
    "f = wrl.util.get_wradlib_data_file(fpath)\n",
    "data, attrs = wrl.io.readDX(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(10,8))\n",
    "ax, cf = wrl.vis.plot_ppi(data, cmap=\"spectral\")\n",
    "pl.xlabel(\"Easting from radar (km)\")\n",
    "pl.ylabel(\"Northing from radar (km)\")\n",
    "pl.title(\"Radar Feldberg, 2008-06-02 16:55 UTC\")\n",
    "cb = pl.colorbar(cf, shrink=0.8)\n",
    "cb.set_label(\"dBZ\")\n",
    "pl.plot([0,105.6],[0,73.4],\"-\", color=\"white\", lw=2)\n",
    "pl.xlim(-128,128)\n",
    "pl.ylim(-128,128)\n",
    "pl.grid(color=\"grey\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We see a set of convective cells with high rainfall intensity in the NE-sector of the Feldberg radar. Let us examine the reflectivity profile along **three beams which at azimuths 53-55 degree** (as marked by the white line in the PPI above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# just a little helper function\n",
    "def plot_beams(data, mybeams, sub=111):\n",
    "    ax = fig.add_subplot(sub)\n",
    "    labelsize=13\n",
    "    for beam in range(mybeams.start, mybeams.stop):\n",
    "        pl.plot(data[beam], label=\"{0} deg\".format(beam))\n",
    "    pl.grid()\n",
    "    pl.text(0.99, 0.88, \"Reflectivity along beams\", \n",
    "            horizontalalignment='right', \n",
    "            transform = ax.transAxes, fontsize=\"large\")\n",
    "    pl.xlabel(\"range (km)\", fontsize=\"large\")\n",
    "    pl.ylabel(\"Reflectivity (dBZ)\", fontsize=\"large\")\n",
    "    pl.legend(loc=\"upper left\")\n",
    "    ax.tick_params(axis='x', labelsize=labelsize)\n",
    "    ax.tick_params(axis='y', labelsize=labelsize)\n",
    "    pl.xlim(0,128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mybeams = slice(53,56)\n",
    "fig = pl.figure(figsize=(10,3))\n",
    "plot_beams(data, mybeams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examining different attenuation correction procedures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hitschfeld and Bordan (1954): Unconstrained gate-by-gate retrieval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "First, we examine the behaviour of the \"classical\" unconstrained forward correction which is typically referred to [Hitschfeld et al., 1954)](http://wradlib.org/wradlib-docs/latest/zreferences.html#hitschfeld1954), although Hitschfeld and Bordan themselves rejected this approach. The Path Integrated Attenuation (PIA) according to this approach can be obtained as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "pia_hibo = wrl.atten.correctAttenuationHB(\n",
    "    data, \n",
    "    coefficients = dict(a=8.e-5, b=0.731, gate_length=1.0),\n",
    "    mode=\"warn\", \n",
    "    thrs=59.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the coefficients dictionary, we can pass the power law parameters of the A(Z) relation as well as the gate length (in km). If we pass \"warn\" as the mode argument, we will obtain a warning log in case the corrected reflectivity exceeds the value of argument ``thrs`` (dBZ).\n",
    "\n",
    "Plotting the result below the reflectivity profile, we obtain the following figure.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# another little helper function\n",
    "def plot_pia(pia, sub=111, title=None):\n",
    "    ax = fig.add_subplot(sub)\n",
    "    labelsize=13\n",
    "    pl.plot(pia.T)\n",
    "    pl.grid()\n",
    "    pl.ylim(0,30)\n",
    "    pl.ylabel(\"PIA (dB)\", fontsize=\"large\")\n",
    "    pl.text(0.01, 0.88, title, \n",
    "            transform = ax.transAxes, fontsize=\"large\")\n",
    "    ax.tick_params(axis='x', labelsize=labelsize)\n",
    "    ax.tick_params(axis='y', labelsize=labelsize)\n",
    "    pl.xlim(0,128)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,6))\n",
    "plot_beams(data, mybeams, 211)\n",
    "plot_pia(pia_hibo[mybeams], 212, \n",
    "         \"PIA according to Hitchfeld and Bordan\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Apparently, slight differences in the reflectivity profile can cause a dramatic change in the behaviour. While at 54 and 55 degrees, the retrieval of PIA appears to be fairly stable, the profile of PIA for 53 degree demonstrates a case of instability.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Harrison et al. (2000): Cap PIA at a factor of two"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "[Harrison et al., 2000](http://wradlib.org/wradlib-docs/latest/zreferences.html#harrison2000) suggested to simply cap PIA in case it would cause a correction of rainfall intensity by more than a factor of two. Depending on the parameters of the Z(R) relationship, that would correpond to PIA values between 4 and 5 dB (4.8 dB if we assume exponent b=1.6).\n",
    "\n",
    "One way to implement this approach would be the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "pia_harrison = wrl.atten.correctAttenuationHB(\n",
    "    data, \n",
    "    coefficients = dict(a=4.57e-5, b=0.731, gate_length=1.0),\n",
    "    mode=\"warn\", \n",
    "    thrs=59.)\n",
    "pia_harrison[pia_harrison > 4.8] = 4.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "And the results would look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,6))\n",
    "\n",
    "mybeams = slice(53,56)\n",
    "labelsize=13\n",
    "\n",
    "plot_beams(data, mybeams, 211)\n",
    "plot_pia(pia_harrison[mybeams], 212, \n",
    "         \"PIA according to Harrison\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Kraemer et al. (2008): Iterative estimation of A(Z) relationsip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "[Kraemer et al., 2008](http://wradlib.org/wradlib-docs/latest/zreferences.html#kraemer2008) suggested to iteratively determine the power law parameters of the A(Z). In particular, the power law coefficient is interatively decreased until the attenuation correction does not lead to reflectivity values above a given threshold (Kraemer suggested 59 dBZ). Using $\\omega radlib$, this would be called by using the function [wradlib.atten.correctAttenuationConstrained2()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.atten.correctAttenuationConstrained2.html) with a specific ``constraints`` argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "pia_kraemer = wrl.atten.correctAttenuationConstrained2(\n",
    "    data,\n",
    "    a_max=1.67e-4, \n",
    "    a_min=2.33e-5,\n",
    "    n_a=100,\n",
    "    b_max=0.7,\n",
    "    b_min=0.65,\n",
    "    n_b=6,\n",
    "    gate_length=1.,\n",
    "    constraints=[wrl.atten.constraint_dBZ],\n",
    "    constraint_args=[[59.0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In brief, this call specifies ranges of the power parameters a and b of the A(Z) relation. Beginning from the maximum values (``a_max`` and ``b_max``), the function searches for values of ``a`` and ``b`` so that the corrected reflectivity will not exceed the dBZ constraint of 59 dBZ. Compared to the previous results, the corresponding profiles of PIA look like this:   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,6))\n",
    "plot_beams(data, mybeams, 211)\n",
    "plot_pia(pia_kraemer[mybeams], 212, \n",
    "         \"PIA according to Kraemer\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Generalised Kraemer procedure: adding additional constraints "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The function [wradlib.atten.correctAttenuationConstrained2()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.atten.correctAttenuationConstrained2.html) allows us to pass any kind of constraint function or lists of constraint functions via the argument ``constraints``. The arguments of these functions are passed via a nested list as argument ``constraint_args``. For example, [Jacobi et al., 2016](http://wradlib.org/wradlib-docs/latest/zreferences.html#jacobi2015) suggested to constrain *both* the corrected reflectivity (by a maximum of 59 dBZ) *and* the resulting path-intgrated attenuation PIA (by a maximum of 20 dB):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "pia_mkraemer = wrl.atten.correctAttenuationConstrained2(\n",
    "    data,\n",
    "    a_max=1.67e-4,\n",
    "    a_min=2.33e-5,\n",
    "    n_a=100,\n",
    "    b_max=0.7,\n",
    "    b_min=0.65,\n",
    "    n_b=6,\n",
    "    gate_length=1.,\n",
    "    constraints=[wrl.atten.constraint_dBZ, \n",
    "                 wrl.atten.constraint_pia], \n",
    "    constraint_args=[[59.0],[20.0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,6))\n",
    "plot_beams(data, mybeams, 211)\n",
    "plot_pia(pia_mkraemer[mybeams], 212, \n",
    "         \"PIA according to modified Kraemer\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparison of all Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Plotting all of the above methods ([Hitschfeld and Bordan](#Hitschfeld-and-Bordan), [Harrison](#Harrison), [Kraemer](#Kraemer), [Modified Kraemer](#Modified-Kraemer) allows for a better comparison of their behaviour. Please refer to [Jacobi et al., 2016](http://wradlib.org/wradlib-docs/latest/zreferences.html#jacobi2015) for an in-depth discussion of this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,12))\n",
    "\n",
    "plot_beams(data, mybeams, 511)\n",
    "plot_pia(pia_hibo[mybeams], 512, \n",
    "         \"PIA according to Hitschfeld and Bordan\" )\n",
    "plot_pia(pia_harrison[mybeams], 513, \n",
    "         \"PIA according to Harrison\")\n",
    "plot_pia(pia_kraemer[mybeams], 514, \n",
    "         \"PIA according to Kraemer\")\n",
    "plot_pia(pia_mkraemer[mybeams], 515, \n",
    "         \"PIA according to modified Kraemer\")\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
}