{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Cascade decomposition\n\nThis example script shows how to compute and plot the cascade decompositon of \na single radar precipitation field in pysteps.\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.cascade.bandpass_filters import filter_gaussian\nfrom pysteps import io, rcparams\nfrom pysteps.cascade.decomposition import decomposition_fft\nfrom pysteps.utils import conversion, transformation\nfrom pysteps.visualization import plot_precip_field"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read precipitation field\n\nFirst thing,  the radar composite is imported and transformed in units\nof dB.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Import the example radar composite\nroot_path = rcparams.data_sources[\"fmi\"][\"root_path\"]\nfilename = os.path.join(\n    root_path, \"20160928\", \"201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz\"\n)\nR, _, metadata = io.import_fmi_pgm(filename, gzipped=True)\n\n# Convert to rain rate\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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2D Fourier spectrum\n\nCompute and plot the 2D Fourier power spectrum of the precipitaton field.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set Nans as the fill value\nR[~np.isfinite(R)] = metadata[\"zerovalue\"]\n\n# Compute the Fourier transform of the input field\nF = abs(np.fft.fftshift(np.fft.fft2(R)))\n\n# Plot the power spectrum\nM, N = F.shape\nfig, ax = plt.subplots()\nim = ax.imshow(\n    np.log(F ** 2), vmin=4, vmax=24, cmap=cm.jet, extent=(-N / 2, N / 2, -M / 2, M / 2)\n)\ncb = fig.colorbar(im)\nax.set_xlabel(\"Wavenumber $k_x$\")\nax.set_ylabel(\"Wavenumber $k_y$\")\nax.set_title(\"Log-power spectrum of R\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Cascade decomposition\n\nFirst, construct a set of Gaussian bandpass filters and plot the corresponding\n1D filters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "num_cascade_levels = 7\n\n# Construct the Gaussian bandpass filters\nfilter = filter_gaussian(R.shape, num_cascade_levels)\n\n# Plot the bandpass filter weights\nL = max(N, M)\nfig, ax = plt.subplots()\nfor k in range(num_cascade_levels):\n    ax.semilogx(\n        np.linspace(0, L / 2, len(filter[\"weights_1d\"][k, :])),\n        filter[\"weights_1d\"][k, :],\n        \"k-\",\n        basex=pow(0.5 * L / 3, 1.0 / (num_cascade_levels - 2)),\n    )\nax.set_xlim(1, L / 2)\nax.set_ylim(0, 1)\nxt = np.hstack([[1.0], filter[\"central_wavenumbers\"][1:]])\nax.set_xticks(xt)\nax.set_xticklabels([\"%.2f\" % cf for cf in filter[\"central_wavenumbers\"]])\nax.set_xlabel(\"Radial wavenumber $|\\mathbf{k}|$\")\nax.set_ylabel(\"Normalized weight\")\nax.set_title(\"Bandpass filter weights\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, apply the 2D Gaussian filters to decompose the radar rainfall field\ninto a set of cascade levels of decreasing spatial scale and plot them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "decomp = decomposition_fft(R, filter, compute_stats=True)\n\n# Plot the normalized cascade levels\nfor i in range(num_cascade_levels):\n    mu = decomp[\"means\"][i]\n    sigma = decomp[\"stds\"][i]\n    decomp[\"cascade_levels\"][i] = (decomp[\"cascade_levels\"][i] - mu) / sigma\n\nfig, ax = plt.subplots(nrows=2, ncols=4)\n\nax[0, 0].imshow(R, cmap=cm.RdBu_r, vmin=-5, vmax=5)\nax[0, 1].imshow(decomp[\"cascade_levels\"][0], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[0, 2].imshow(decomp[\"cascade_levels\"][1], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[0, 3].imshow(decomp[\"cascade_levels\"][2], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 0].imshow(decomp[\"cascade_levels\"][3], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 1].imshow(decomp[\"cascade_levels\"][4], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 2].imshow(decomp[\"cascade_levels\"][5], cmap=cm.RdBu_r, vmin=-3, vmax=3)\nax[1, 3].imshow(decomp[\"cascade_levels\"][6], cmap=cm.RdBu_r, vmin=-3, vmax=3)\n\nax[0, 0].set_title(\"Observed\")\nax[0, 1].set_title(\"Level 1\")\nax[0, 2].set_title(\"Level 2\")\nax[0, 3].set_title(\"Level 3\")\nax[1, 0].set_title(\"Level 4\")\nax[1, 1].set_title(\"Level 5\")\nax[1, 2].set_title(\"Level 6\")\nax[1, 3].set_title(\"Level 7\")\n\nfor i in range(2):\n    for j in range(4):\n        ax[i, j].set_xticks([])\n        ax[i, j].set_yticks([])\nplt.tight_layout()\nplt.show()\n\n# sphinx_gallery_thumbnail_number = 4"
      ]
    }
  ],
  "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
}