examples/my_first_nowcast.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L_dntwSQBnbK"
   },
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/pySTEPS/pysteps/blob/master/examples/my_first_nowcast.ipynb)\n",
    "\n",
    "# My first precipitation nowcast\n",
    "\n",
    "In this example, we will use pysteps to compute and plot an extrapolation nowcast using the NSSL's Multi-Radar/Multi-Sensor System\n",
    "([MRMS](https://www.nssl.noaa.gov/projects/mrms/)) rain rate product.\n",
    "\n",
    "The MRMS precipitation product is available every 2 minutes, over the contiguous US. \n",
    "Each precipitation composite has 3500 x 7000 grid points, separated 1 km from each other.\n",
    "\n",
    "## Set-up Colab environment\n",
    "\n",
    "**Important**: In colab, execute this section one cell at a time. Trying to excecute all the cells at once may results in cells being skipped and some dependencies not being installed.\n",
    "\n",
    "First, let's set up our working environment. Note that these steps are only needed to work with google colab. \n",
    "\n",
    "To install pysteps locally, you can follow [these instructions](https://pysteps.readthedocs.io/en/latest/user_guide/install_pysteps.html).\n",
    "\n",
    "First, let's install the latest Pysteps version from the Python Package Index (PyPI) using pip. This will also install the minimal dependencies needed to run pysteps. \n",
    "\n",
    "#### Install optional dependencies\n",
    "\n",
    "Now, let's install the optional dependendies that will allow us to plot and read the example data.\n",
    "- pygrib: to read the MRMS data grib format\n",
    "- pyproj: needed by pygrib\n",
    "\n",
    "**NOTE:** Do not import pysteps in this notebook until the following optional dependencies are loaded. Otherwise, pysteps will assume that they are not installed and some of its functionalities won't work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mFx4hq_DBtp-"
   },
   "outputs": [],
   "source": [
    "# These libraries are needed for the pygrib library in Colab.\n",
    "# Note that is needed if you install pygrib using pip.\n",
    "# If you use conda, the libraries will be installed automatically.\n",
    "! apt-get install libeccodes-dev libproj-dev\n",
    "\n",
    "# Install the python packages\n",
    "! pip install pyproj\n",
    "! pip install pygrib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6BF2paxnTuGB"
   },
   "outputs": [],
   "source": [
    "# Uninstall existing shapely\n",
    "# We will re-install shapely in the next step by ignoring the binary\n",
    "# wheels to make it compatible with other modules that depend on\n",
    "# GEOS, such as Cartopy (used here).\n",
    "!pip uninstall --yes shapely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "7x8Hx_4hE_BU"
   },
   "outputs": [],
   "source": [
    "# To install cartopy in Colab using pip, we need to install the library\n",
    "# dependencies first.\n",
    "\n",
    "!apt-get install -qq libgdal-dev libgeos-dev\n",
    "!pip install shapely --no-binary shapely\n",
    "!pip install cartopy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ybD55ZJhmdYa"
   },
   "source": [
    "#### Install pysteps\n",
    "\n",
    "Now that all dependencies are installed, we can install pysteps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "VA7zp3nRmhfF"
   },
   "outputs": [],
   "source": [
    "# ! pip install git+https://github.com/pySTEPS/pysteps\n",
    "! pip install pysteps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-AkfR6JSBujn"
   },
   "source": [
    "## Getting the example data\n",
    "\n",
    "Now that we have the environment ready, let's install the example data and configure the pysteps's default parameters by following [this tutorial](https://pysteps.readthedocs.io/en/latest/user_guide/example_data.html).\n",
    "\n",
    "First, we will use the [pysteps.datasets.download_pysteps_data()](https://pysteps.readthedocs.io/en/latest/generated/pysteps.datasets.download_pysteps_data.html) function to download the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "vri-R_ZVGihj"
   },
   "outputs": [],
   "source": [
    "# Import the helper functions\n",
    "from pysteps.datasets import download_pysteps_data, create_default_pystepsrc\n",
    "\n",
    "# Download the pysteps data in the \"pysteps_data\"\n",
    "download_pysteps_data(\"pysteps_data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wdKfjliCKXhx"
   },
   "source": [
    "Next, we need to create a default configuration file that points to the downloaded data. \n",
    "By default, pysteps will place the configuration file in `$HOME/.pysteps` (unix and Mac OS X) or `$USERPROFILE/pysteps` (windows).\n",
    "To quickly create a configuration file, we will use the [pysteps.datasets.create_default_pystepsrc()](https://pysteps.readthedocs.io/en/latest/generated/pysteps.datasets.create_default_pystepsrc.html#pysteps.datasets.create_default_pystepsrc) helper function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "pGdKHa36H5JX"
   },
   "outputs": [],
   "source": [
    "# If the configuration file is placed in one of the default locations\n",
    "# (https://pysteps.readthedocs.io/en/latest/user_guide/set_pystepsrc.html#configuration-file-lookup)\n",
    "# it will be loaded automatically when pysteps is imported.\n",
    "config_file_path = create_default_pystepsrc(\"pysteps_data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DAFUJgR5K1CS"
   },
   "source": [
    "Since pysteps was already initialized in this notebook, we need to load the new configuration file and update the default configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "tMIbQLPAK42h"
   },
   "outputs": [],
   "source": [
    "# Import pysteps and load the new configuration file\n",
    "import pysteps\n",
    "\n",
    "_ = pysteps.load_config_file(config_file_path, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SzSqp1DFJ0M9"
   },
   "source": [
    "Let's see what the default parameters look like (these are stored in the\n",
    "[pystepsrc file](https://pysteps.readthedocs.io/en/latest/user_guide/set_pystepsrc.html)). We will be using them to load the MRMS data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6Gr65nH4BnbP"
   },
   "outputs": [],
   "source": [
    "# The default parameters are stored in pysteps.rcparams.\n",
    "from pprint import pprint\n",
    "\n",
    "pprint(pysteps.rcparams.data_sources[\"mrms\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9M_buv7WBnbf"
   },
   "source": [
    "This should have printed the following lines:\n",
    "\n",
    "- `fn_ext`: 'grib2' -- The file extension\n",
    "- `fn_pattern`: 'PrecipRate_00.00_%Y%m%d-%H%M%S' -- The file naming convention of the MRMS data.\n",
    "- `importer`: 'mrms_grib' -- The name of the importer for the MRMS data.\n",
    "- `importer_kwargs`: {} -- Extra options provided to the importer. None in this example.\n",
    "- `path_fmt`: '%Y/%m/%d' -- The folder structure in which the files are stored. Here, year/month/day/filename.\n",
    "- `root_path`: '/content/pysteps_data/mrms' -- The root path of the MRMS-data.\n",
    "- `timestep`: 2 -- The temporal interval of the (radar) rainfall data\n",
    "\n",
    "Note that the default `timestep` parameter is 2 minutes, which corresponds to the time interval at which the MRMS product is available.\n",
    "\n",
    "## Load the MRMS example data\n",
    "\n",
    "Now that we have installed the example data, let's import the example MRMS dataset using the [load_dataset()](https://pysteps.readthedocs.io/en/latest/generated/pysteps.datasets.load_dataset.html) helper function from the `pysteps.datasets` module.\n",
    "\n",
    "We import 1 hour and 10 minutes of data, which corresponds to a sequence of 35 frames of 2-D precipitation composites.\n",
    "Note that importing the data takes approximately 30 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-8Q4e58VBnbl"
   },
   "outputs": [],
   "source": [
    "from pysteps.datasets import load_dataset\n",
    "\n",
    "# We'll import the time module to measure the time the importer needed\n",
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "# Import the data\n",
    "precipitation, metadata, timestep = load_dataset(\n",
    "    \"mrms\", frames=35\n",
    ")  # precipitation in mm/h\n",
    "\n",
    "end_time = time.time()\n",
    "\n",
    "print(\"Precipitation data imported\")\n",
    "print(\"Importing the data took \", (end_time - start_time), \" seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "btiTxYYMBnby"
   },
   "source": [
    "Let's have a look at the values returned by the `load_dataset()` function. \n",
    "\n",
    "- `precipitation`: A numpy array with (time, latitude, longitude) dimensions.\n",
    "- `metadata`: A dictionary with additional information (pixel sizes, map projections, etc.).\n",
    "- `timestep`: Time separation between each sample (in minutes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WqUHbJ_qBnb3"
   },
   "outputs": [],
   "source": [
    "# Let's inspect the shape of the imported data array\n",
    "precipitation.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xa8woT0ABncD"
   },
   "source": [
    "Note that the shape of the precipitation is 4 times smaller than the raw MRMS data (3500 x 7000).\n",
    "The `load_dataset()` function uses the default parameters from `importers` to read the data. By default, the MRMS importer upscales the data 4x. That is, from ~1km resolution to ~4km. It also uses single precision to reduce the memory requirements.\n",
    "Thanks to the upscaling, the memory footprint of this example dataset is ~200Mb instead of the 3.1Gb of the raw (3500 x 7000) data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "22O2YXrfBncG"
   },
   "outputs": [],
   "source": [
    "timestep  # In minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "J8_4hwcXBncT"
   },
   "outputs": [],
   "source": [
    "pprint(metadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uQREORtJBnch"
   },
   "source": [
    "# Time to make a nowcast\n",
    "\n",
    "So far, we have 1 hour and 10 minutes of precipitation images, separated 2 minutes from each other.\n",
    "But, how do we use that data to run a precipitation forecast? \n",
    "\n",
    "A simple way is by extrapolating the precipitation field, assuming it will continue to move as observed in the recent past, and without changes in intensity. This is commonly known as *Lagrangian persistence*.\n",
    "\n",
    "The first step to run our nowcast based on Lagrangian persistence, is the estimation of the motion field from a sequence of past precipitation observations.\n",
    "We use the Lucas-Kanade (LK) optical flow method implemented in pysteps.\n",
    "This method follows a local tracking approach that relies on the OpenCV package.\n",
    "Local features are tracked in a sequence of two or more radar images.\n",
    "The scheme includes a final interpolation step to produce a smooth field of motion vectors.\n",
    "Other optical flow methods are also available in pysteps. \n",
    "Check the full list [here](https://pysteps.readthedocs.io/en/latest/pysteps_reference/motion.html).\n",
    "\n",
    "Now let's use the first 5 precipitation images (10 min) to estimate the motion field of the radar pattern and the remaining 30 images (1h) to evaluate the quality of our forecast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "jcb2Sf6xBnck"
   },
   "outputs": [],
   "source": [
    "# precipitation[0:5] -> Used to find motion (past data). Let's call it training precip.\n",
    "train_precip = precipitation[0:5]\n",
    "\n",
    "# precipitation[5:] -> Used to evaluate forecasts (future data, not available in \"real\" forecast situation)\n",
    "# Let's call it observed precipitation because we will use it to compare our forecast with the actual observations.\n",
    "observed_precip = precipitation[3:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xt1TbB0RBncu"
   },
   "source": [
    "Let's see what this 'training' precipitation event looks like using the [pysteps.visualization.plot_precip_field](https://pysteps.readthedocs.io/en/latest/generated/pysteps.visualization.precipfields.plot_precip_field.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bmNYLo1jBncw"
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from pysteps.visualization import plot_precip_field\n",
    "\n",
    "# Set a figure size that looks nice ;)\n",
    "plt.figure(figsize=(9, 5), dpi=100)\n",
    "\n",
    "# Plot the last rainfall field in the \"training\" data.\n",
    "# train_precip[-1] -> Last available composite for nowcasting.\n",
    "plot_precip_field(train_precip[-1], geodata=metadata, axis=\"off\")\n",
    "plt.show()  # (This line is actually not needed if you are using jupyter notebooks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NVRfJm11Bnc7"
   },
   "source": [
    "Did you note the **shaded grey** regions? Those are the regions were no valid observations where available to estimate the precipitation (e.g., due to ground clutter, no radar coverage, or radar beam blockage).\n",
    "Those regions need to be handled with care when we run our nowcast.\n",
    "\n",
    "### Data exploration\n",
    "\n",
    "Before we produce a forecast, let's explore the precipitation data. In particular, let's see how the distribution of the rain rate values looks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WER6RttPBnc9"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Let's define some plotting default parameters for the next plots\n",
    "# Note: This is not strictly needed.\n",
    "plt.rc(\"figure\", figsize=(4, 4))\n",
    "plt.rc(\"figure\", dpi=100)\n",
    "plt.rc(\"font\", size=14)  # controls default text sizes\n",
    "plt.rc(\"axes\", titlesize=14)  # fontsize of the axes title\n",
    "plt.rc(\"axes\", labelsize=14)  # fontsize of the x and y labels\n",
    "plt.rc(\"xtick\", labelsize=14)  # fontsize of the tick labels\n",
    "plt.rc(\"ytick\", labelsize=14)  # fontsize of the tick labels\n",
    "\n",
    "# Let's use the last available composite for nowcasting from the \"training\" data (train_precip[-1])\n",
    "# Also, we will discard any invalid value.\n",
    "valid_precip_values = train_precip[-1][~np.isnan(train_precip[-1])]\n",
    "\n",
    "# Plot the histogram\n",
    "bins = np.concatenate(([-0.01, 0.01], np.linspace(1, 40, 39)))\n",
    "plt.hist(valid_precip_values, bins=bins, log=True, edgecolor=\"black\")\n",
    "plt.autoscale(tight=True, axis=\"x\")\n",
    "plt.xlabel(\"Rainfall intensity [mm/h]\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.title(\"Precipitation rain rate histogram in mm/h units\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "O6TvIXS3BndH"
   },
   "source": [
    "The histogram shows that rain rate values have a non-Gaussian and asymmetric distribution that is bounded at zero. Also, the probability of occurrence decays extremely fast with increasing rain rate values (note the logarithmic y-axis).\n",
    "\n",
    "\n",
    "For better performance of the motion estimation algorithms, we can convert the rain rate values (in mm/h) to a more log-normal distribution  of rain rates by applying the following logarithmic transformation:\n",
    "\n",
    "\\begin{equation}\n",
    "R\\rightarrow\n",
    "\\begin{cases}\n",
    "    10\\log_{10}R, & \\text{if } R\\geq 0.1\\text{mm h$^{-1}$} \\\\\n",
    "    -15,          & \\text{otherwise}\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "The transformed precipitation corresponds to logarithmic rain rates in units of dBR. The value of −15 dBR is equivalent to assigning a rain rate of approximately 0.03 mm h$^{−1}$ to the zeros. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hgA4PeapBndK"
   },
   "outputs": [],
   "source": [
    "from pysteps.utils import transformation\n",
    "\n",
    "# Log-transform the data to dBR.\n",
    "# The threshold of 0.1 mm/h sets the fill value to -15 dBR.\n",
    "train_precip_dbr, metadata_dbr = transformation.dB_transform(\n",
    "    train_precip, metadata, threshold=0.1, zerovalue=-15.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Nx3VESBlBndU"
   },
   "source": [
    "Let's inspect the resulting **transformed precipitation** distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "rYS5hBIGBndX"
   },
   "outputs": [],
   "source": [
    "# Only use the valid data!\n",
    "valid_precip_dbr = train_precip_dbr[-1][~np.isnan(train_precip_dbr[-1])]\n",
    "\n",
    "plt.figure(figsize=(4, 4), dpi=100)\n",
    "\n",
    "# Plot the histogram\n",
    "counts, bins, _ = plt.hist(valid_precip_dbr, bins=40, log=True, edgecolor=\"black\")\n",
    "plt.autoscale(tight=True, axis=\"x\")\n",
    "plt.xlabel(\"Rainfall intensity [dB]\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.title(\"Precipitation rain rate histogram in dB units\")\n",
    "\n",
    "# Let's add a lognormal distribution that fits that data to the plot.\n",
    "import scipy\n",
    "\n",
    "bin_center = (bins[1:] + bins[:-1]) * 0.5\n",
    "bin_width = np.diff(bins)\n",
    "\n",
    "# We will only use one composite to fit the function to speed up things.\n",
    "# First, remove the no precip areas.\"\n",
    "precip_to_fit = valid_precip_dbr[valid_precip_dbr > -15]\n",
    "\n",
    "fit_params = scipy.stats.lognorm.fit(precip_to_fit)\n",
    "\n",
    "fitted_pdf = scipy.stats.lognorm.pdf(bin_center, *fit_params)\n",
    "\n",
    "# Multiply pdf by the bin width and the total number of grid points: pdf -> total counts per bin.\n",
    "fitted_pdf = fitted_pdf * bin_width * precip_to_fit.size\n",
    "\n",
    "# Plot the log-normal fit\n",
    "plt.plot(bin_center, fitted_pdf, label=\"Fitted log-normal\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZocO5zqUBndg"
   },
   "source": [
    "That looks more like a log-normal distribution. Note the large peak at -15dB. That peak corresponds to \"zero\" (below threshold) precipitation. The jump with no data in between -15 and -10 dB is caused by the precision of the data, which we had set to 1 decimal. Hence, the lowest precipitation intensities (above zero) are 0.1 mm/h (= -10 dB).\n",
    "\n",
    "## Compute the nowcast\n",
    "\n",
    "These are the minimal steps to compute a short-term forecast using Lagrangian extrapolation of the precipitation patterns:\n",
    " \n",
    " 1. Estimate the precipitation motion field.\n",
    " 1. Use the motion field to advect the most recent radar rainfall field and produce an extrapolation forecast.\n",
    "\n",
    "### Estimate the motion field\n",
    "\n",
    "Now we can estimate the motion field. Here we use a local feature-tracking approach (Lucas-Kanade).\n",
    "However, check the other methods available in the [pysteps.motion](https://pysteps.readthedocs.io/en/latest/pysteps_reference/motion.html) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mnACmX_0Bndi"
   },
   "outputs": [],
   "source": [
    "# Estimate the motion field with Lucas-Kanade\n",
    "from pysteps import motion\n",
    "from pysteps.visualization import plot_precip_field, quiver\n",
    "\n",
    "# Import the Lucas-Kanade optical flow algorithm\n",
    "oflow_method = motion.get_method(\"LK\")\n",
    "\n",
    "# Estimate the motion field from the training data (in dBR)\n",
    "motion_field = oflow_method(train_precip_dbr)\n",
    "\n",
    "## Plot the motion field.\n",
    "# Use a figure size that looks nice ;)\n",
    "plt.figure(figsize=(9, 5), dpi=100)\n",
    "plt.title(\"Estimated motion field with the Lukas-Kanade algorithm\")\n",
    "\n",
    "# Plot the last rainfall field in the \"training\" data.\n",
    "# Remember to use the mm/h precipitation data since plot_precip_field assumes\n",
    "# mm/h by default. You can change this behavior using the \"units\" keyword.\n",
    "plot_precip_field(train_precip[-1], geodata=metadata, axis=\"off\")\n",
    "\n",
    "# Plot the motion field vectors\n",
    "quiver(motion_field, geodata=metadata, step=40)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YObddRFCBnd1"
   },
   "source": [
    "### Extrapolate the observations\n",
    "\n",
    "We have all ingredients to make an extrapolation nowcast now. \n",
    "The final step is to advect the most recent radar rainfall field along the estimated motion field, producing an extrapolation forecast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "erSLAzvNBnd3"
   },
   "outputs": [],
   "source": [
    "from pysteps import nowcasts\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "# Extrapolate the last radar observation\n",
    "extrapolate = nowcasts.get_method(\"extrapolation\")\n",
    "\n",
    "# You can use the precipitation observations directly in mm/h for this step.\n",
    "last_observation = train_precip[-1]\n",
    "\n",
    "last_observation[~np.isfinite(last_observation)] = metadata[\"zerovalue\"]\n",
    "\n",
    "# We set the number of leadtimes (the length of the forecast horizon) to the\n",
    "# length of the observed/verification preipitation data. In this way, we'll get\n",
    "# a forecast that covers these time intervals.\n",
    "n_leadtimes = observed_precip.shape[0]\n",
    "\n",
    "# Advect the most recent radar rainfall field and make the nowcast.\n",
    "precip_forecast = extrapolate(train_precip[-1], motion_field, n_leadtimes)\n",
    "\n",
    "# This shows the shape of the resulting array with [time intervals, rows, cols]\n",
    "print(\"The shape of the resulting array is: \", precip_forecast.shape)\n",
    "\n",
    "end = time.time()\n",
    "print(\"Advecting the radar rainfall fields took \", (end - start), \" seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "csy5s-yRBneB"
   },
   "source": [
    "Let's inspect the last forecast time (hence this is the forecast rainfall an hour ahead)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "MUiS5-HPBneD"
   },
   "outputs": [],
   "source": [
    "# Plot precipitation at the end of the forecast period.\n",
    "plt.figure(figsize=(9, 5), dpi=100)\n",
    "plot_precip_field(precip_forecast[-1], geodata=metadata, axis=\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mQEseXvhBneI"
   },
   "source": [
    "## Evaluate the forecast quality\n",
    "\n",
    "Many verification methods are already present in pysteps (see a complete list [here](https://pysteps.readthedocs.io/en/latest/pysteps_reference/verification.html)). We just have to import them. \n",
    "\n",
    "Here, we will evaluate our forecast using the Fractions Skill Score (FSS). \n",
    "This metric provides an intuitive assessment of the dependency of forecast skill on spatial scale and intensity. This makes the FSS an ideal skill score for high-resolution precipitation forecasts.\n",
    "\n",
    "More precisely, the FSS is a neighborhood spatial verification method that directly compares the fractional coverage of events in windows surrounding the observations and forecasts.\n",
    "The FSS varies from 0 (total mismatch) to 1 (perfect forecast).\n",
    "For most situations, an FSS value of > 0.5 serves as a good indicator of a useful forecast ([Roberts and Lean, 2008](https://journals.ametsoc.org/doi/full/10.1175/2007MWR2123.1) and [Skok and Roberts, 2016](https://rmets.onlinelibrary.wiley.com/doi/full/10.1002/qj.2849)). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "No3qBjqSBneK"
   },
   "outputs": [],
   "source": [
    "from pysteps import verification\n",
    "\n",
    "fss = verification.get_method(\"FSS\")\n",
    "\n",
    "# Compute fractions skill score (FSS) for all lead times for different scales using a 1 mm/h detection threshold.\n",
    "scales = [\n",
    "    2,\n",
    "    4,\n",
    "    8,\n",
    "    16,\n",
    "    32,\n",
    "    64,\n",
    "]  # In grid points.\n",
    "\n",
    "scales_in_km = np.array(scales) * 4\n",
    "\n",
    "# Set the threshold\n",
    "thr = 1.0  # in mm/h\n",
    "\n",
    "score = []\n",
    "\n",
    "# Calculate the FSS for every lead time and all predefined scales.\n",
    "for i in range(n_leadtimes):\n",
    "    score_ = []\n",
    "    for scale in scales:\n",
    "        score_.append(\n",
    "            fss(precip_forecast[i, :, :], observed_precip[i, :, :], thr, scale)\n",
    "        )\n",
    "    score.append(score_)\n",
    "\n",
    "# Now plot it\n",
    "plt.figure()\n",
    "x = np.arange(1, n_leadtimes + 1) * timestep\n",
    "plt.plot(x, score, lw=2.0)\n",
    "plt.xlabel(\"Lead time [min]\")\n",
    "plt.ylabel(\"FSS ( > 1.0 mm/h ) \")\n",
    "plt.title(\"Fractions Skill Score\")\n",
    "plt.legend(\n",
    "    scales_in_km,\n",
    "    title=\"Scale [km]\",\n",
    "    loc=\"center left\",\n",
    "    bbox_to_anchor=(1.01, 0.5),\n",
    "    bbox_transform=plt.gca().transAxes,\n",
    ")\n",
    "plt.autoscale(axis=\"x\", tight=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the FSS decreases with increasing lead time.\n",
    "This is expected, as the forecasting quality slowly decreases when we forecast further ahead.\n",
    "Upscaling the forecast, however, clearly leads to higher skill (up to longer ahead) compared to the forecast on the highest resolutions.\n",
    "\n",
    "## Concluding remarks\n",
    "Congratulations, you have successfully made your first nowcast using the pysteps library!\n",
    "This was a simple extrapolation-based nowcast and a lot more advanced options are possible too, see [the pysteps examples gallery](https://pysteps.readthedocs.io/en/latest/auto_examples/index.html) for some nice examples."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "my_first_nowcast.ipynb",
   "provenance": []
  },
  "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.2"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}