{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Precipitation downscaling with RainFARM\n\nThis example script shows how to use the stochastic downscaling method RainFARM\navailable in pysteps.\n\nRainFARM is a downscaling algorithm for rainfall fields developed by Rebora et\nal. (2006). The method can represent the realistic small-scale variability of the\ndownscaled precipitation field by means of Gaussian random fields.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport os\nfrom pprint import pprint\n\nfrom pysteps import io, rcparams\nfrom pysteps.utils import aggregate_fields_space, square_domain, to_rainrate\nfrom pysteps.downscaling import rainfarm\nfrom pysteps.visualization import plot_precip_field"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read the input data\n\nAs first step, we need to import the precipitation field that we are going\nto use in this example.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Import the example radar composite\nroot_path = rcparams.data_sources[\"mch\"][\"root_path\"]\nfilename = os.path.join(root_path, \"20160711\", \"AQC161932100V_00005.801.gif\")\nprecip, _, metadata = io.import_mch_gif(\n    filename, product=\"AQC\", unit=\"mm\", accutime=5.0\n)\n\n# Convert to mm/h\nprecip, metadata = to_rainrate(precip, metadata)\n\n# Reduce to a square domain\nprecip, metadata = square_domain(precip, metadata, \"crop\")\n\n# Nicely print the metadata\npprint(metadata)\n\n# Plot the original rainfall field\nplot_precip_field(precip, geodata=metadata)\nplt.show()\n\n# Assign the fill value to all the Nans\nprecip[~np.isfinite(precip)] = metadata[\"zerovalue\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Upscale the field\n\nTo test our downscaling method, we first need to upscale the original field to\na lower resolution. We are going to use an upscaling factor of 16 x.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "upscaling_factor = 16\nupscale_to = metadata[\"xpixelsize\"] * upscaling_factor  # upscaling factor : 16 x\nprecip_lr, metadata_lr = aggregate_fields_space(precip, metadata, upscale_to)\n\n# Plot the upscaled rainfall field\nplt.figure()\nplot_precip_field(precip_lr, geodata=metadata_lr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Downscale the field\n\nWe can now use RainFARM to generate stochastic realizations of the downscaled\nprecipitation field.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(5, 8))\n# Set the number of stochastic realizations\nnum_realizations = 5\n\n# Per realization, generate a stochastically downscaled precipitation field\n# and plot it.\n# The first time, the spectral slope alpha needs to be estimated. To illustrate\n# the sensitity of this parameter, we are going to plot some realizations with\n# half or double the estimated slope.\nalpha = None\nfor n in range(num_realizations):\n\n    # Spectral slope estimated from the upscaled field\n    precip_hr, alpha = rainfarm.downscale(\n        precip_lr, alpha=alpha, ds_factor=upscaling_factor, return_alpha=True\n    )\n    plt.subplot(num_realizations, 3, n * 3 + 2)\n    plot_precip_field(precip_hr, geodata=metadata, axis=\"off\", colorbar=False)\n    if n == 0:\n        plt.title(f\"alpha={alpha:.1f}\")\n\n    # Half the estimated slope\n    precip_hr = rainfarm.downscale(\n        precip_lr, alpha=alpha * 0.5, ds_factor=upscaling_factor\n    )\n    plt.subplot(num_realizations, 3, n * 3 + 1)\n    plot_precip_field(precip_hr, geodata=metadata, axis=\"off\", colorbar=False)\n    if n == 0:\n        plt.title(f\"alpha={alpha * 0.5:.1f}\")\n\n    # Double the estimated slope\n    precip_hr = rainfarm.downscale(\n        precip_lr, alpha=alpha * 2, ds_factor=upscaling_factor\n    )\n    plt.subplot(num_realizations, 3, n * 3 + 3)\n    plot_precip_field(precip_hr, geodata=metadata, axis=\"off\", colorbar=False)\n    if n == 0:\n        plt.title(f\"alpha={alpha * 2:.1f}\")\n\n    plt.subplots_adjust(wspace=0, hspace=0)\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Remarks\n\nCurrently, the pysteps implementation of RainFARM only covers spatial downscaling.\nThat is, it can improve the spatial resolution of a rainfall field. However, unlike\nthe original algorithm from Rebora et al. (2006), it cannot downscale the temporal\ndimension.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n\nRebora, N., L. Ferraris, J. von Hardenberg, and A. Provenzale, 2006: RainFARM:\nRainfall downscaling by a filtered autoregressive model. J. Hydrometeor., 7,\n724\u2013738.\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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}