{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Thunderstorm Detection and Tracking - DATing\n\nThis example shows how to use the thunderstorm DATing module. The example is based on\nMeteoSwiss radar data and uses the Cartesian composite of maximum reflectivity on a\n1 km grid. All default values are tuned to this grid, but can be modified.\nThe first section demonstrates thunderstorm cell detection and how to plot contours.\nThe second section demonstrates detection and tracking in combination,\nas well as how to plot the resulting tracks.\nThis module was implemented following the procedures used in the TRT Thunderstorms\nRadar Tracking algorithm (:cite:`TRT2004`) used operationally at MeteoSwiss.\nModifications include advecting the identified thunderstorms with the optical flow\nobtained from pysteps, as well as additional options in the thresholding.\n\n## References\n:cite:`TRT2004`\n\n@author: mfeldman\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Import all required functions\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\nfrom pprint import pprint\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom pysteps import io, rcparams\nfrom pysteps.feature import tstorm as tstorm_detect\nfrom pysteps.tracking import tdating as tstorm_dating\nfrom pysteps.utils import to_reflectivity\nfrom pysteps.visualization import plot_precip_field, plot_track, plot_cart_contour"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Read the radar input images\n\nA series of 20 files containing Swiss Cartesian gridded rain rates are imported. Since\nthe algorithm is tuned to Swiss max-reflectivity data, the rain rates are transformed\nto reflectivity fields using the 'to_reflectivity' utility in pysteps.utils.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Select the input data\ndate = datetime.strptime(\"201607112100\", \"%Y%m%d%H%M\")\ndata_source = rcparams.data_sources[\"mch\"]\n\n# Extract corresponding settings\nroot_path = data_source[\"root_path\"]\npath_fmt = data_source[\"path_fmt\"]\nfn_pattern = data_source[\"fn_pattern\"]\nfn_ext = data_source[\"fn_ext\"]\nimporter_name = data_source[\"importer\"]\nimporter_kwargs = data_source[\"importer_kwargs\"]\ntimestep = data_source[\"timestep\"]\n\n# Load the data from the archive\nfns = io.archive.find_by_date(\n    date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_next_files=20\n)\nimporter = io.get_method(importer_name, \"importer\")\nR, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Convert to reflectivity (it is possible to give the a- and b- parameters of the\n# Marshall-Palmer relationship here: zr_a = and zr_b =).\nZ, metadata = to_reflectivity(R, metadata)\n\n# Extract the list of timestamps\ntimelist = metadata[\"timestamps\"]\n\npprint(metadata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Example of thunderstorm identification in a single timestep\nThe function tstorm_detect.detection requires a 2-D input image, all further inputs are\noptional.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "input_image = Z[2, :, :].copy()\ntime = timelist[2]\ncells_id, labels = tstorm_detect.detection(input_image, time=time)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Properties of one of the identified cells:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(cells_id.iloc[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Example of thunderstorm tracking over a timeseries\nThe tstorm-dating function requires the entire pre-loaded time series.\nThe first two timesteps are required to initialize the\nflow prediction and are not used to compute tracks.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "track_list, cell_list, label_list = tstorm_dating.dating(\n    input_video=Z, timelist=timelist\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Plotting the results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Plot precipitation field\nplot_precip_field(Z[2, :, :], geodata=metadata, units=metadata[\"unit\"])\nplt.xlabel(\"Swiss easting [m]\")\nplt.ylabel(\"Swiss northing [m]\")\n\n# Add the identified cells\nplot_cart_contour(cells_id.cont, geodata=metadata)\n\n# Filter the tracks to only contain cells existing in this timestep\nIDs = cells_id.ID.values\ntrack_filt = []\nfor track in track_list:\n    if np.unique(track.ID) in IDs:\n        track_filt.append(track)\n\n# Add their tracks\nplot_track(track_filt, geodata=metadata)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Evaluating temporal behaviour of cell\nMaximum reflectivity of cells in time\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Make an empty list\ntlen = []\n# Get a list of colors that we will use for the plot\ncolor = iter(plt.cm.ocean(np.linspace(0, 0.8, len(track_filt))))\n# Now, loop through all the tracks and plot the maximum reflectivity of the cell\n# in time.\nfor track in track_filt:\n    plt.plot(np.arange(len(track)), track.max_ref, c=next(color))\n    tlen.append(len(track))\nplt.xticks(np.arange(max(tlen) + 1), labels=np.arange(max(tlen) + 1) * 5)\nplt.ylabel(\"Maximum reflectivity (dBZ)\")\nplt.xlabel(\"Time since cell detection (min)\")\nplt.legend(IDs, loc=\"lower right\", ncol=3, title=\"Track number\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The size of the thunderstorm cells in time\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Make an empty list\ntlen = []\n# Get a list of colors that we will use for the plot\ncolor = iter(plt.cm.ocean(np.linspace(0, 0.8, len(track_filt))))\n# Now, loop through all the tracks and plot the cell size of the thunderstorms\n# in time.\nfor track in track_filt:\n    size = []\n    for ID, t in track.iterrows():\n        size.append(len(t.x))\n    plt.plot(np.arange(len(track)), size, c=next(color))\n    tlen.append(len(track))\nplt.xticks(np.arange(max(tlen) + 1), labels=np.arange(max(tlen) + 1) * 5)\nplt.ylabel(\"Thunderstorm cell size (pixels)\")\nplt.xlabel(\"Time since cell detection (min)\")\nplt.legend(IDs, loc=\"upper left\", ncol=3, title=\"Track number\")\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
}