{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Extrapolation nowcast\n\nThis tutorial shows how to compute and plot an extrapolation nowcast using \nFinnish radar data.\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\nfrom pprint import pprint\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\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(\"201609281600\", \"%Y%m%d%H%M\")\ndata_source = rcparams.data_sources[\"fmi\"]\nn_leadtimes = 12"
      ]
    },
    {
      "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, num_prev_files=2\n)\n\n# Read the radar composites\nimporter = io.get_method(importer_name, \"importer\")\nZ, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Convert to rain rate\nR, metadata = conversion.to_rainrate(Z, metadata)\n\n# Plot the rainfall field\nplot_precip_field(R[-1, :, :], geodata=metadata)\nplt.show()\n\n# Store the last frame for plotting it later later\nR_ = R[-1, :, :].copy()\n\n# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,\n# set the fill value to -15 dBR\nR, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)\n\n# Nicely print the metadata\npprint(metadata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute the nowcast\n\nThe extrapolation nowcast is based on the estimation of the motion field,\nwhich is here performed using a local tracking approach (Lucas-Kanade).\nThe most recent radar rainfall field is then simply advected along this motion\nfield in oder to produce an extrapolation forecast.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Estimate the motion field with Lucas-Kanade\noflow_method = motion.get_method(\"LK\")\nV = oflow_method(R[-3:, :, :])\n\n# Extrapolate the last radar observation\nextrapolate = nowcasts.get_method(\"extrapolation\")\nR[~np.isfinite(R)] = metadata[\"zerovalue\"]\nR_f = extrapolate(R[-1, :, :], V, n_leadtimes)\n\n# Back-transform to rain rate\nR_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]\n\n# Plot the motion field\nplot_precip_field(R_, geodata=metadata)\nquiver(V, geodata=metadata, step=50)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify with FSS\n\nThe fractions skill score (FSS) provides an intuitive assessment of the\ndependency of skill on spatial scale and intensity, which makes it an ideal\nskill score for high-resolution precipitation forecasts.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Find observations in the data archive\nfns = io.archive.find_by_date(\n    date,\n    root_path,\n    path_fmt,\n    fn_pattern,\n    fn_ext,\n    timestep,\n    num_prev_files=0,\n    num_next_files=n_leadtimes,\n)\n# Read 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, 223.0, 1.53)\n\n# Compute fractions skill score (FSS) for all lead times, a set of scales and 1 mm/h\nfss = verification.get_method(\"FSS\")\nscales = [2, 4, 8, 16, 32, 64, 128, 256, 512]\nthr = 1.0\nscore = []\nfor i in range(n_leadtimes):\n    score_ = []\n    for scale in scales:\n        score_.append(fss(R_f[i, :, :], R_o[i + 1, :, :], thr, scale))\n    score.append(score_)\n\nplt.figure()\nx = np.arange(1, n_leadtimes + 1) * timestep\nplt.plot(x, score)\nplt.legend(scales, title=\"Scale [km]\")\nplt.xlabel(\"Lead time [min]\")\nplt.ylabel(\"FSS ( > 1.0 mm/h ) \")\nplt.title(\"Fractions skill score\")\nplt.show()\n\n# sphinx_gallery_thumbnail_number = 3"
      ]
    }
  ],
  "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
}