{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Generation of stochastic noise\n\nThis example script shows how to run the stochastic noise field generators\nincluded in pysteps.\n\nThese noise fields are used as perturbation terms during an extrapolation\nnowcast in order to represent the uncertainty in the evolution of the rainfall\nfield.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib import cm, pyplot as plt\nimport numpy as np\nimport os\nfrom pprint import pprint\nfrom pysteps import io, rcparams\nfrom pysteps.noise.fftgenerators import initialize_param_2d_fft_filter\nfrom pysteps.noise.fftgenerators import initialize_nonparam_2d_fft_filter\nfrom pysteps.noise.fftgenerators import generate_noise_2d_fft_filter\nfrom pysteps.utils import conversion, rapsd, transformation\nfrom pysteps.visualization import plot_precip_field, plot_spectrum1d"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read precipitation field\n\nFirst thing,  the radar composite is imported and transformed in units\nof dB.\nThis image will be used to train the Fourier filters that are necessary to\nproduce the fields of spatially correlated noise.\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\")\nR, _, metadata = io.import_mch_gif(filename, product=\"AQC\", unit=\"mm\", accutime=5.0)\n\n# Convert to mm/h\nR, metadata = conversion.to_rainrate(R, metadata)\n\n# Nicely print the metadata\npprint(metadata)\n\n# Plot the rainfall field\nplot_precip_field(R, geodata=metadata)\nplt.show()\n\n# Log-transform the data\nR, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)\n\n# Assign the fill value to all the Nans\nR[~np.isfinite(R)] = metadata[\"zerovalue\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Parametric filter\n\nIn the parametric approach, a power-law model is used to approximate the power\nspectral density (PSD) of a given rainfall field.\n\nThe parametric model uses  a  piece-wise  linear  function  with  two  spectral\nslopes (beta1 and beta2) and one breaking point\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Fit the parametric PSD to the observation\nFp = initialize_param_2d_fft_filter(R)\n\n# Compute the observed and fitted 1D PSD\nL = np.max(Fp[\"input_shape\"])\nif L % 2 == 1:\n    wn = np.arange(0, int(L / 2) + 1)\nelse:\n    wn = np.arange(0, int(L / 2))\nR_, freq = rapsd(R, fft_method=np.fft, return_freq=True)\nf = np.exp(Fp[\"model\"](np.log(wn), *Fp[\"pars\"]))\n\n# Extract the scaling break in km, beta1 and beta2\nw0 = L / np.exp(Fp[\"pars\"][0])\nb1 = Fp[\"pars\"][2]\nb2 = Fp[\"pars\"][3]\n\n# Plot the observed power spectrum and the model\nfig, ax = plt.subplots()\nplot_scales = [512, 256, 128, 64, 32, 16, 8, 4]\nplot_spectrum1d(\n    freq,\n    R_,\n    x_units=\"km\",\n    y_units=\"dBR\",\n    color=\"k\",\n    ax=ax,\n    label=\"Observed\",\n    wavelength_ticks=plot_scales,\n)\nplot_spectrum1d(\n    freq,\n    f,\n    x_units=\"km\",\n    y_units=\"dBR\",\n    color=\"r\",\n    ax=ax,\n    label=\"Fit\",\n    wavelength_ticks=plot_scales,\n)\nplt.legend()\nax.set_title(\n    \"Radially averaged log-power spectrum of R\\n\"\n    r\"$\\omega_0=%.0f km, \\beta_1=%.1f, \\beta_2=%.1f$\" % (w0, b1, b2)\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Nonparametric filter\n\nIn the nonparametric approach,  the Fourier filter is obtained directly\nfrom the power spectrum of the observed precipitation field R.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Fnp = initialize_nonparam_2d_fft_filter(R)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Noise generator\n\nThe parametric and nonparametric filters obtained above can now be used\nto produce N realizations of random fields of prescribed power spectrum,\nhence with the same correlation structure as the initial rainfall field.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "seed = 42\nnum_realizations = 3\n\n# Generate noise\nNp = []\nNnp = []\nfor k in range(num_realizations):\n    Np.append(generate_noise_2d_fft_filter(Fp, seed=seed + k))\n    Nnp.append(generate_noise_2d_fft_filter(Fnp, seed=seed + k))\n\n# Plot the generated noise fields\n\nfig, ax = plt.subplots(nrows=2, ncols=3)\n\n# parametric noise\nax[0, 0].imshow(Np[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[0, 1].imshow(Np[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[0, 2].imshow(Np[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)\n\n# nonparametric noise\nax[1, 0].imshow(Nnp[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 1].imshow(Nnp[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 2].imshow(Nnp[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)\n\nfor i in range(2):\n    for j in range(3):\n        ax[i, j].set_xticks([])\n        ax[i, j].set_yticks([])\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above figure highlights the main limitation of the parametric approach\n(top row), that is, the assumption of an isotropic power law scaling\nrelationship, meaning that anisotropic structures such as rainfall bands\ncannot be represented.\n\nInstead, the nonparametric approach (bottom row) allows generating\nperturbation fields with anisotropic  structures, but it also requires a\nlarger sample size and is sensitive to the quality of the input data, e.g.\nthe presence of residual clutter in the radar image.\n\nIn addition, both techniques assume spatial stationarity of the covariance\nstructure of the field.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# 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
}