{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# ANVIL nowcast\n\nThis example demonstrates how to use ANVIL and the advantages compared to\nextrapolation nowcast and S-PROG.\n\nLoad the libraries.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime, timedelta\nimport warnings\n\nwarnings.simplefilter(\"ignore\")\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pysteps import motion, io, rcparams, utils\nfrom pysteps.nowcasts import anvil, extrapolation, sprog\nfrom pysteps.utils import transformation\nfrom pysteps.visualization import plot_precip_field"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read the input data\n\nANVIL was originally developed to use vertically integrated liquid (VIL) as\nthe input data, but the model allows using any two-dimensional input fields.\nHere we use a composite of rain rates.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "date = datetime.strptime(\"201505151620\", \"%Y%m%d%H%M\")\n\n# Read the data source information from rcparams\ndata_source = rcparams.data_sources[\"mch\"]\n\nroot_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\"]\n\n# Find the input files in the archive. Use history length of 5 timesteps\nfilenames = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_prev_files=5\n)\n\n# Read the input time series\nimporter = io.get_method(importer_name, \"importer\")\nrainrate_field, quality, metadata = io.read_timeseries(\n    filenames, importer, **importer_kwargs\n)\n\n# Convert to rain rate (mm/h)\nrainrate_field, metadata = utils.to_rainrate(rainrate_field, metadata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute the advection field\n\nApply the Lucas-Kanade method with the parameters given in Pulkkinen et al.\n(2020) to compute the advection field.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fd_kwargs = {}\nfd_kwargs[\"max_corners\"] = 1000\nfd_kwargs[\"quality_level\"] = 0.01\nfd_kwargs[\"min_distance\"] = 2\nfd_kwargs[\"block_size\"] = 8\n\nlk_kwargs = {}\nlk_kwargs[\"winsize\"] = (15, 15)\n\noflow_kwargs = {}\noflow_kwargs[\"fd_kwargs\"] = fd_kwargs\noflow_kwargs[\"lk_kwargs\"] = lk_kwargs\noflow_kwargs[\"decl_scale\"] = 10\n\noflow = motion.get_method(\"lucaskanade\")\n\n# transform the input data to logarithmic scale\nrainrate_field_log, _ = utils.transformation.dB_transform(\n    rainrate_field, metadata=metadata\n)\nvelocity = oflow(rainrate_field_log, **oflow_kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Compute the nowcasts and threshold rain rates below 0.5 mm/h\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "forecast_extrap = extrapolation.forecast(\n    rainrate_field[-1], velocity, 3, extrap_kwargs={\"allow_nonfinite_values\": True}\n)\nforecast_extrap[forecast_extrap < 0.5] = 0.0\n\n# log-transform the data and the threshold value to dBR units for S-PROG\nrainrate_field_db, _ = transformation.dB_transform(\n    rainrate_field, metadata, threshold=0.1, zerovalue=-15.0\n)\nrainrate_thr, _ = transformation.dB_transform(\n    np.array([0.5]), metadata, threshold=0.1, zerovalue=-15.0\n)\nforecast_sprog = sprog.forecast(\n    rainrate_field_db[-3:], velocity, 3, n_cascade_levels=8, R_thr=rainrate_thr[0]\n)\nforecast_sprog, _ = transformation.dB_transform(\n    forecast_sprog, threshold=-10.0, inverse=True\n)\nforecast_sprog[forecast_sprog < 0.5] = 0.0\n\nforecast_anvil = anvil.forecast(\n    rainrate_field[-4:], velocity, 3, ar_window_radius=25, ar_order=2\n)\nforecast_anvil[forecast_anvil < 0.5] = 0.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Read the reference observation field and threshold rain rates below 0.5 mm/h\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "filenames = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=3\n)\n\nrefobs_field, quality, metadata = io.read_timeseries(\n    filenames, importer, **importer_kwargs\n)\n\nrefobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata)\nrefobs_field[refobs_field < 0.5] = 0.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the extrapolation, S-PROG and ANVIL nowcasts.\n\nFor comparison, the observed rain rate fields are also plotted. Growth and\ndecay areas are marked with red and blue circles, respectively.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_growth_decay_circles(ax):\n    circle = plt.Circle(\n        (360, 300), 25, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (420, 350), 30, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (405, 380), 30, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (420, 500), 25, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (480, 535), 30, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (330, 470), 35, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (505, 205), 30, color=\"b\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (440, 180), 30, color=\"r\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (590, 240), 30, color=\"r\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n    circle = plt.Circle(\n        (585, 160), 15, color=\"r\", clip_on=False, fill=False, zorder=1e9\n    )\n    ax.add_artist(circle)\n\n\nfig = plt.figure(figsize=(10, 13))\n\nax = fig.add_subplot(321)\nrainrate_field[-1][rainrate_field[-1] < 0.5] = 0.0\nplot_precip_field(rainrate_field[-1])\nplot_growth_decay_circles(ax)\nax.set_title(\"Obs. %s\" % str(date))\n\nax = fig.add_subplot(322)\nplot_precip_field(refobs_field)\nplot_growth_decay_circles(ax)\nax.set_title(\"Obs. %s\" % str(date + timedelta(minutes=15)))\n\nax = fig.add_subplot(323)\nplot_precip_field(forecast_extrap[-1])\nplot_growth_decay_circles(ax)\nax.set_title(\"Extrapolation +15 minutes\")\n\nax = fig.add_subplot(324)\nplot_precip_field(forecast_sprog[-1])\nplot_growth_decay_circles(ax)\nax.set_title(\"S-PROG (with post-processing)\\n +15 minutes\")\n\nax = fig.add_subplot(325)\nplot_precip_field(forecast_anvil[-1])\nplot_growth_decay_circles(ax)\nax.set_title(\"ANVIL +15 minutes\")\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Remarks\n\nThe extrapolation nowcast is static, i.e. it does not predict any growth or\ndecay. While S-PROG is to some extent able to predict growth and decay, this\nthis comes with loss of small-scale features. In addition, statistical\npost-processing needs to be applied to correct the bias and incorrect wet-area\nratio introduced by the autoregressive process. ANVIL is able to do both:\npredict growth and decay and preserve the small-scale structure in a way that\npost-processing is not necessary.\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
}