{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# LINDA nowcasts\n\nThis example shows how to compute and plot a deterministic and ensemble LINDA\nnowcasts using Swiss radar data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\nimport warnings\n\nwarnings.simplefilter(\"ignore\")\n\nimport matplotlib.pyplot as plt\n\nfrom pysteps import io, rcparams\nfrom pysteps.motion.lucaskanade import dense_lucaskanade\nfrom pysteps.nowcasts import linda, sprog, steps\nfrom pysteps.utils import conversion, dimension, transformation\nfrom pysteps.visualization import plot_precip_field"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read the input rain rate fields\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "date = datetime.strptime(\"201701311200\", \"%Y%m%d%H%M\")\ndata_source = \"mch\"\n\n# Read the data source information from rcparams\ndatasource_params = rcparams.data_sources[data_source]\n\n# Find the radar files in the archive\nfns = io.find_by_date(\n    date,\n    datasource_params[\"root_path\"],\n    datasource_params[\"path_fmt\"],\n    datasource_params[\"fn_pattern\"],\n    datasource_params[\"fn_ext\"],\n    datasource_params[\"timestep\"],\n    num_prev_files=2,\n)\n\n# Read the data from the archive\nimporter = io.get_method(datasource_params[\"importer\"], \"importer\")\nreflectivity, _, metadata = io.read_timeseries(\n    fns, importer, **datasource_params[\"importer_kwargs\"]\n)\n\n# Convert reflectivity to rain rate\nrainrate, metadata = conversion.to_rainrate(reflectivity, metadata)\n\n# Upscale data to 2 km to reduce computation time\nrainrate, metadata = dimension.aggregate_fields_space(rainrate, metadata, 2000)\n\n# Plot the most recent rain rate field\nplt.figure()\nplot_precip_field(rainrate[-1, :, :])\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Estimate the advection field\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# The advection field is estimated using the Lucas-Kanade optical flow\nadvection = dense_lucaskanade(rainrate, verbose=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Deterministic nowcast\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute 30-minute LINDA nowcast with 8 parallel workers\n# Restrict the number of features to 15 to reduce computation time\nnowcast_linda = linda.forecast(\n    rainrate,\n    advection,\n    6,\n    max_num_features=15,\n    add_perturbations=False,\n    num_workers=8,\n    measure_time=True,\n)[0]\n\n# Compute S-PROG nowcast for comparison\nrainrate_db, _ = transformation.dB_transform(\n    rainrate, metadata, threshold=0.1, zerovalue=-15.0\n)\nnowcast_sprog = sprog.forecast(\n    rainrate_db[-3:, :, :],\n    advection,\n    6,\n    n_cascade_levels=6,\n    R_thr=-10.0,\n)\n\n# Convert reflectivity nowcast to rain rate\nnowcast_sprog = transformation.dB_transform(\n    nowcast_sprog, threshold=-10.0, inverse=True\n)[0]\n\n# Plot the nowcasts\nfig = plt.figure(figsize=(9, 4))\nax = fig.add_subplot(1, 2, 1)\nplot_precip_field(\n    nowcast_linda[-1, :, :],\n    title=\"LINDA (+ 30 min)\",\n)\n\nax = fig.add_subplot(1, 2, 2)\nplot_precip_field(\n    nowcast_sprog[-1, :, :],\n    title=\"S-PROG (+ 30 min)\",\n)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above figure shows that the filtering scheme implemented in LINDA preserves\nsmall-scale and band-shaped features better than S-PROG. This is because the\nformer uses a localized elliptical convolution kernel instead of the\ncascade-based autoregressive process, where the parameters are estimated over\nthe whole domain.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Probabilistic nowcast\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute 30-minute LINDA nowcast ensemble with 40 members and 8 parallel workers\nnowcast_linda = linda.forecast(\n    rainrate,\n    advection,\n    6,\n    max_num_features=15,\n    add_perturbations=True,\n    num_ens_members=40,\n    num_workers=8,\n    measure_time=True,\n)[0]\n\n# Compute 40-member STEPS nowcast for comparison\nnowcast_steps = steps.forecast(\n    rainrate_db[-3:, :, :],\n    advection,\n    6,\n    40,\n    n_cascade_levels=6,\n    R_thr=-10.0,\n    mask_method=\"incremental\",\n    kmperpixel=2.0,\n    timestep=datasource_params[\"timestep\"],\n    vel_pert_method=None,\n)\n\n# Convert reflectivity nowcast to rain rate\nnowcast_steps = transformation.dB_transform(\n    nowcast_steps, threshold=-10.0, inverse=True\n)[0]\n\n# Plot two ensemble members of both nowcasts\nfig = plt.figure()\nfor i in range(2):\n    ax = fig.add_subplot(2, 2, i + 1)\n    ax = plot_precip_field(\n        nowcast_linda[i, -1, :, :], geodata=metadata, colorbar=False, axis=\"off\"\n    )\n    ax.set_title(f\"LINDA Member {i+1}\")\n\nfor i in range(2):\n    ax = fig.add_subplot(2, 2, 3 + i)\n    ax = plot_precip_field(\n        nowcast_steps[i, -1, :, :], geodata=metadata, colorbar=False, axis=\"off\"\n    )\n    ax.set_title(f\"STEPS Member {i+1}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above figure shows the main difference between LINDA and STEPS. In\naddition to the convolution kernel, another improvement in LINDA is a\nlocalized perturbation generator using the short-space Fourier transform\n(SSFT) and a spatially variable marginal distribution. As a result, the\nLINDA ensemble members preserve the anisotropic and small-scale structures\nconsiderably better than STEPS.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.tight_layout()\nplt.show()"
      ]
    }
  ],
  "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
}