{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Probability forecasts\n\nThis example script shows how to forecast the probability of exceeding an\nintensity threshold.\n\nThe method is based on the local Lagrangian approach described in Germann and\nZawadzki (2004).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom pysteps.nowcasts.lagrangian_probability import forecast\nfrom pysteps.visualization import plot_precip_field" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical example\n\nFirst, we use some dummy data to show the basic principle of this approach.\nThe probability forecast is produced by sampling a spatial neighborhood that is\nincreased as a function of lead time. As a result, the edges of\nthe yellow square becomes more and more smooth as t increases. This represents\nthe strong loss of predictability with lead time of any extrapolation nowcast.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# parameters\nprecip = np.zeros((100, 100))\nprecip[10:50, 10:50] = 1\nvelocity = np.ones((2, *precip.shape))\ntimesteps = [0, 2, 6, 12]\nthr = 0.5\nslope = 1 # pixels / timestep\n\n# compute probability forecast\nout = forecast(precip, velocity, timesteps, thr, slope=slope)\n# plot\nfor n, frame in enumerate(out):\n plt.subplot(2, 2, n + 1)\n plt.imshow(frame, interpolation=\"nearest\", vmin=0, vmax=1)\n plt.title(f\"t={timesteps[n]}\")\n plt.xticks([])\n plt.yticks([])\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Real-data example\n\nWe now apply the same method to real data. We use a slope of 1 km / minute\nas suggested by Germann and Zawadzki (2004), meaning that after 30 minutes,\nthe probabilities are computed by using all pixels within a neighborhood of 30\nkilometers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from datetime import datetime\n\nfrom pysteps import io, rcparams, utils\nfrom pysteps.motion.lucaskanade import dense_lucaskanade\nfrom pysteps.verification import reldiag_init, reldiag_accum, plot_reldiag\n\n# data source\nsource = rcparams.data_sources[\"mch\"]\nroot = rcparams.data_sources[\"mch\"][\"root_path\"]\nfmt = rcparams.data_sources[\"mch\"][\"path_fmt\"]\npattern = rcparams.data_sources[\"mch\"][\"fn_pattern\"]\next = rcparams.data_sources[\"mch\"][\"fn_ext\"]\ntimestep = rcparams.data_sources[\"mch\"][\"timestep\"]\nimporter_name = rcparams.data_sources[\"mch\"][\"importer\"]\nimporter_kwargs = rcparams.data_sources[\"mch\"][\"importer_kwargs\"]\n\n# read precip field\ndate = datetime.strptime(\"201607112100\", \"%Y%m%d%H%M\")\nfns = io.find_by_date(date, root, fmt, pattern, ext, timestep, num_prev_files=2)\nimporter = io.get_method(importer_name, \"importer\")\nprecip, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\nprecip, metadata = utils.to_rainrate(precip, metadata)\n# precip[np.isnan(precip)] = 0\n\n# motion\nmotion = dense_lucaskanade(precip)\n\n# parameters\nnleadtimes = 6\nthr = 1 # mm / h\nslope = 1 * timestep # km / min\n\n# compute probability forecast\nextrap_kwargs = dict(allow_nonfinite_values=True)\nfct = forecast(\n precip[-1], motion, nleadtimes, thr, slope=slope, extrap_kwargs=extrap_kwargs\n)\n\n# plot\nfor n, frame in enumerate(fct):\n plt.subplot(2, 3, n + 1)\n plt.imshow(frame, interpolation=\"nearest\", vmin=0, vmax=1)\n plt.xticks([])\n plt.yticks([])\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot one single leadtime in more detail using the pysteps visualization\nfunctionality.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.close()\n# Plot the field of probabilities\nplot_precip_field(\n fct[2],\n geodata=metadata,\n ptype=\"prob\",\n probthr=thr,\n title=\"Exceedence probability (+ %i min)\" % (nleadtimes * timestep),\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verification\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# verifying observations\nimporter = io.get_method(importer_name, \"importer\")\nfns = io.find_by_date(\n date, root, fmt, pattern, ext, timestep, num_next_files=nleadtimes\n)\nobs, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\nobs, metadata = utils.to_rainrate(obs, metadata)\nobs[np.isnan(obs)] = 0\n\n# reliability diagram\nreldiag = reldiag_init(thr)\nreldiag_accum(reldiag, fct, obs[1:])\nfig, ax = plt.subplots()\nplot_reldiag(reldiag, ax)\nax.set_title(\"Reliability diagram\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\nGermann, U. and I. Zawadzki, 2004:\nScale Dependence of the Predictability of Precipitation from Continental\nRadar Images. Part II: Probability Forecasts.\nJournal of Applied Meteorology, 43(1), 74-89.\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 }