{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Advection correction\n\nThis tutorial shows how to use the optical flow routines of pysteps to implement\nthe advection correction procedure described in Anagnostou and Krajewski (1999).\n\nAdvection correction is a temporal interpolation procedure that is often used\nwhen estimating rainfall accumulations to correct for the shift of rainfall patterns\nbetween consecutive radar rainfall maps. This shift becomes particularly \nsignificant for long radar scanning cycles and in presence of fast moving\nprecipitation features.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The code for the advection correction using pysteps was originally\n          written by `Daniel Wolfensberger <https://github.com/wolfidan>`_.</p></div>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom pysteps import io, motion, rcparams\nfrom pysteps.utils import conversion, dimension\nfrom pysteps.visualization import plot_precip_field\nfrom scipy.ndimage import map_coordinates"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read the radar input images\n\nFirst, we import a sequence of 36 images of 5-minute radar composites\nthat we will use to produce a 3-hour rainfall accumulation map.\nWe will keep only one frame every 10 minutes, to simulate a longer scanning\ncycle and thus better highlight the need for advection correction.\n\nYou need the pysteps-data archive downloaded and the pystepsrc file\nconfigured with the data_source paths pointing to data folders.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Selected case\ndate = datetime.strptime(\"201607112100\", \"%Y%m%d%H%M\")\ndata_source = rcparams.data_sources[\"mch\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Load the data from the archive\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "root_path = data_source[\"root_path\"]\npath_fmt = data_source[\"path_fmt\"]\nfn_pattern = data_source[\"fn_pattern\"]\nfn_ext = data_source[\"fn_ext\"]\nimporter_name = data_source[\"importer\"]\nimporter_kwargs = data_source[\"importer_kwargs\"]\ntimestep = data_source[\"timestep\"]\n\n# Find the input files from the archive\nfns = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=35\n)\n\n# Read the radar composites\nimporter = io.get_method(importer_name, \"importer\")\nR, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Convert to mm/h\nR, metadata = conversion.to_rainrate(R, metadata)\n\n# Upscale to 2 km (simply to reduce the memory demand)\nR, metadata = dimension.aggregate_fields_space(R, metadata, 2000)\n\n# Keep only one frame every 10 minutes (i.e., every 2 timesteps)\n# (to highlight the need for advection correction)\nR = R[::2]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Advection correction\n\nNow we need to implement the advection correction for a pair of successive\nradar images. The procedure is based on the algorithm described in Anagnostou\nand Krajewski (Appendix A, 1999).\n\nTo evaluate the advection occurred between two successive radar images, we are\ngoing to use the Lucas-Kanade optical flow routine available in pysteps.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def advection_correction(R, T=5, t=1):\n    \"\"\"\n    R = np.array([qpe_previous, qpe_current])\n    T = time between two observations (5 min)\n    t = interpolation timestep (1 min)\n    \"\"\"\n\n    # Evaluate advection\n    oflow_method = motion.get_method(\"LK\")\n    fd_kwargs = {\"buffer_mask\": 10}  # avoid edge effects\n    V = oflow_method(np.log(R), fd_kwargs=fd_kwargs)\n\n    # Perform temporal interpolation\n    Rd = np.zeros((R[0].shape))\n    x, y = np.meshgrid(\n        np.arange(R[0].shape[1], dtype=float), np.arange(R[0].shape[0], dtype=float)\n    )\n    for i in range(t, T + t, t):\n\n        pos1 = (y - i / T * V[1], x - i / T * V[0])\n        R1 = map_coordinates(R[0], pos1, order=1)\n\n        pos2 = (y + (T - i) / T * V[1], x + (T - i) / T * V[0])\n        R2 = map_coordinates(R[1], pos2, order=1)\n\n        Rd += (T - i) * R1 + i * R2\n\n    return t / T ** 2 * Rd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we apply the advection correction to the whole sequence of radar\nimages and produce the rainfall accumulation map.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "R_ac = R[0].copy()\nfor i in range(R.shape[0] - 1):\n    R_ac += advection_correction(R[i : (i + 2)], T=10, t=1)\nR_ac /= R.shape[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Results\n\nWe compare the two accumulation maps. The first map on the left is\ncomputed without advection correction and we can therefore see that the shift\nbetween successive images 10 minutes apart produces irregular accumulations.\nConversely, the rainfall accumulation of the right is produced using advection\ncorrection to account for this spatial shift. The final result is a smoother\nrainfall accumulation map.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(9, 4))\nplt.subplot(121)\nplot_precip_field(R.mean(axis=0), title=\"3-h rainfall accumulation\")\nplt.subplot(122)\nplot_precip_field(R_ac, title=\"Same with advection correction\")\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Reference\n\nAnagnostou, E. N., and W. F. Krajewski. 1999. \"Real-Time Radar Rainfall\nEstimation. Part I: Algorithm Formulation.\" Journal of Atmospheric and\nOceanic Technology 16: 189\u201397.\nhttps://doi.org/10.1175/1520-0426(1999)016<0189:RTRREP>2.0.CO;2\n\n"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}