{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nHandling of no-data in Lucas-Kanade\n===================================\n\nAreas of missing data in radar images are typically caused by visibility limits\nsuch as beam blockage and the radar coverage itself. These artifacts can mislead \nthe echo tracking algorithms. For instance, precipitation leaving the domain\nmight be erroneously detected as having nearly stationary velocity.\n\nThis example shows how the Lucas-Kanade algorithm can be tuned to avoid the \nerroneous interpretation of velocities near the maximum range of the radars by \nbuffering the no-data mask in the radar image in order to exclude all vectors \ndetected nearby no-data areas.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\nfrom matplotlib import cm, colors\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom pysteps import io, motion, nowcasts, rcparams, verification\nfrom pysteps.utils import conversion, transformation\nfrom pysteps.visualization import plot_precip_field, quiver"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read the radar input images\n---------------------------\n\nFirst, we will import the sequence of radar composites.\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\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 two input files from the archive\nfns = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_prev_files=1\n)\n\n# Read the radar composites\nimporter = io.get_method(importer_name, \"importer\")\nR, quality, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\ndel quality  # Not used"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Preprocess the data\n~~~~~~~~~~~~~~~~~~~\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Convert to mm/h\nR, metadata = conversion.to_rainrate(R, metadata)\n\n# Keep the reference frame in mm/h and its mask (for plotting purposes)\nref_mm = R[0, :, :].copy()\nmask = np.ones(ref_mm.shape)\nmask[~np.isnan(ref_mm)] = np.nan\n\n# Log-transform the data [dBR]\nR, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)\n\n# Keep the reference frame in dBR (for plotting purposes)\nref_dbr = R[0].copy()\nref_dbr[ref_dbr < -10] = np.nan\n\n# Plot the reference field\nplot_precip_field(ref_mm, title=\"Reference field\")\ncircle = plt.Circle((620, 400), 100, color=\"b\", clip_on=False, fill=False)\nplt.gca().add_artist(circle)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Notice the \"half-in, half-out\" precipitation area within the blue circle.\nAs we are going to show next, the tracking algorithm can erroneously interpret\nprecipitation leaving the domain as stationary motion.\n\nAlso note that the radar image includes NaNs in areas of missing data.\nThese are used by the optical flow algorithm to define the radar mask.\n\nSparse Lucas-Kanade\n-------------------\n\nBy setting the optional argument 'dense=False' in 'x,y,u,v = LK_optflow(.....)',\nthe LK algorithm returns the motion vectors detected by the Lucas-Kanade scheme\nwithout interpolating them on the grid.\nThis allows us to better identify the presence of wrongly detected\nstationary motion in areas where precipitation is leaving the domain (look\nfor the red dots within the blue circle in the figure below).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get Lucas-Kanade optical flow method\nLK_optflow = motion.get_method(\"LK\")\n\n# Mask invalid values\nR = np.ma.masked_invalid(R)\nR.data[R.mask] = np.nan\n\n# Use default settings (i.e., no buffering of the radar mask)\nx, y, u, v = LK_optflow(R, dense=False, buffer_mask=0, quality_level_ST=0.1)\nplt.imshow(ref_dbr, cmap=plt.get_cmap(\"Greys\"))\nplt.imshow(mask, cmap=colors.ListedColormap([\"black\"]), alpha=0.5)\nplt.quiver(x, y, u, v, color=\"red\", angles=\"xy\", scale_units=\"xy\", scale=0.2)\ncircle = plt.Circle((620, 245), 100, color=\"b\", clip_on=False, fill=False)\nplt.gca().add_artist(circle)\nplt.title(\"buffer_mask = 0 (default)\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By default, the LK algorithm considers missing values as no precipitation, i.e.,\nno-data are the same as no-echoes. As a result, the fixed boundaries produced\nby precipitation in contact with no-data areas are interpreted as stationary motion.\nOne way to mitigate this effect of the boundaries is to introduce a slight buffer\nof the no-data mask so that the algorithm will ignore all the portions of the\nradar domain that are nearby no-data areas.\nThis is achieved by setting the keyword argument 'buffer_mask = 20' in\n'x,y,u,v = LK_optflow(.....)'.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# with buffer\nx, y, u, v = LK_optflow(R, dense=False, buffer_mask=20, quality_level_ST=0.2)\nplt.imshow(ref_dbr, cmap=plt.get_cmap(\"Greys\"))\nplt.imshow(mask, cmap=colors.ListedColormap([\"black\"]), alpha=0.5)\nplt.quiver(x, y, u, v, color=\"red\", angles=\"xy\", scale_units=\"xy\", scale=0.2)\ncircle = plt.Circle((620, 245), 100, color=\"b\", clip_on=False, fill=False)\nplt.gca().add_artist(circle)\nplt.title(\"buffer_mask = 20\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Dense Lucas-Kanade\n------------------\n\nThe above displacement vectors produced by the Lucas-Kanade method are now\ninterpolated to produce a full field of motion (i.e., 'dense=True').\nBy comparing the velocity of the motion fields, we can easily notice\nthe negative bias that is introduced by the the erroneous interpretation of\nvelocities near the maximum range of the radars.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "UV1 = LK_optflow(R, dense=True, buffer_mask=0, quality_level_ST=0.1)\nUV2 = LK_optflow(R, dense=True, buffer_mask=20, quality_level_ST=0.2)\n\nV1 = np.sqrt(UV1[0] ** 2 + UV1[1] ** 2)\nV2 = np.sqrt(UV2[0] ** 2 + UV2[1] ** 2)\n\nplt.imshow((V1 - V2) / V2, cmap=cm.RdBu_r, vmin=-0.1, vmax=0.1)\nplt.colorbar(fraction=0.04, pad=0.04)\nplt.title(\"Relative difference in motion speed\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Notice that the default motion field can be significantly slower (more than 10%\nslower) because of the presence of wrong stationary motion vectors at the edges.\n\nForecast skill\n--------------\n\nWe are now going to evaluate the benefit of buffering the radar mask by computing\nthe forecast skill in terms of the Spearman correlation coefficient.\nThe extrapolation forecasts are computed using the dense UV motion fields\nestimted above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Get the advection routine and extrapolate the last radar frame by 12 time steps\n# (i.e., 1 hour lead time)\nextrapolate = nowcasts.get_method(\"extrapolation\")\nR[~np.isfinite(R)] = metadata[\"zerovalue\"]\nR_f1 = extrapolate(R[-1], UV1, 12)\nR_f2 = extrapolate(R[-1], UV2, 12)\n\n# Back-transform to rain rate\nR_f1 = transformation.dB_transform(R_f1, threshold=-10.0, inverse=True)[0]\nR_f2 = transformation.dB_transform(R_f2, threshold=-10.0, inverse=True)[0]\n\n# Find the veriyfing observations in the archive\nfns = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=12\n)\n\n# Read and convert the radar composites\nR_o, _, metadata_o = io.read_timeseries(fns, importer, **importer_kwargs)\nR_o, metadata_o = conversion.to_rainrate(R_o, metadata_o)\n\n# Compute Spearman correlation\nskill = verification.get_method(\"corr_s\")\nscore_1 = []\nscore_2 = []\nfor i in range(12):\n    score_1.append(skill(R_f1[i, :, :], R_o[i + 1, :, :])[\"corr_s\"])\n    score_2.append(skill(R_f2[i, :, :], R_o[i + 1, :, :])[\"corr_s\"])\n\nx = (np.arange(12) + 1) * 5  # [min]\nplt.plot(x, score_1, label=\"no mask buffer\")\nplt.plot(x, score_2, label=\"with mask buffer\")\nplt.legend()\nplt.xlabel(\"Lead time [min]\")\nplt.ylabel(\"Corr. coeff. []\")\nplt.title(\"Spearman correlation\")\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As expected, the corrected motion field produces better forecast skill already\nwithin the first hour into the nowcast.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2"
      ]
    }
  ],
  "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.7.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}