{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# STEPS nowcast\n\nThis tutorial shows how to compute and plot an ensemble nowcast using Swiss\nradar data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom datetime import datetime\nfrom pprint import pprint\nfrom pysteps import io, nowcasts, rcparams\nfrom pysteps.motion.lucaskanade import dense_lucaskanade\nfrom pysteps.postprocessing.ensemblestats import excprob\nfrom pysteps.utils import conversion, dimension, transformation\nfrom pysteps.visualization import plot_precip_field\n\n# Set nowcast parameters\nn_ens_members = 20\nn_leadtimes = 6\nseed = 24" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read precipitation field\n\nFirst thing, the sequence of Swiss radar composites is imported, converted and\ntransformed into units of dBR.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "date = datetime.strptime(\"201701311200\", \"%Y%m%d%H%M\")\ndata_source = \"mch\"\n\n# Load data source config\nroot_path = rcparams.data_sources[data_source][\"root_path\"]\npath_fmt = rcparams.data_sources[data_source][\"path_fmt\"]\nfn_pattern = rcparams.data_sources[data_source][\"fn_pattern\"]\nfn_ext = rcparams.data_sources[data_source][\"fn_ext\"]\nimporter_name = rcparams.data_sources[data_source][\"importer\"]\nimporter_kwargs = rcparams.data_sources[data_source][\"importer_kwargs\"]\ntimestep = rcparams.data_sources[data_source][\"timestep\"]\n\n# Find the radar files in the archive\nfns = io.find_by_date(\n date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2\n)\n\n# Read the data from the archive\nimporter = io.get_method(importer_name, \"importer\")\nR, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Convert to rain rate\nR, metadata = conversion.to_rainrate(R, metadata)\n\n# Upscale data to 2 km to limit memory usage\nR, metadata = dimension.aggregate_fields_space(R, metadata, 2000)\n\n# Plot the rainfall field\nplot_precip_field(R[-1, :, :], geodata=metadata)\nplt.show()\n\n# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,\n# set the fill value to -15 dBR\nR, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)\n\n# Set missing values with the fill value\nR[~np.isfinite(R)] = -15.0\n\n# Nicely print the metadata\npprint(metadata)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Deterministic nowcast with S-PROG\n\nFirst, the motiong field is estimated using a local tracking approach based\non the Lucas-Kanade optical flow.\nThe motion field can then be used to generate a deterministic nowcast with\nthe S-PROG model, which implements a scale filtering appraoch in order to\nprogressively remove the unpredictable spatial scales during the forecast.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Estimate the motion field\nV = dense_lucaskanade(R)\n\n# The S-PROG nowcast\nnowcast_method = nowcasts.get_method(\"sprog\")\nR_f = nowcast_method(\n R[-3:, :, :],\n V,\n n_leadtimes,\n n_cascade_levels=6,\n R_thr=-10.0,\n)\n\n# Back-transform to rain rate\nR_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]\n\n# Plot the S-PROG forecast\nplot_precip_field(\n R_f[-1, :, :],\n geodata=metadata,\n title=\"S-PROG (+ %i min)\" % (n_leadtimes * timestep),\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from the figure above, the forecast produced by S-PROG is a\nsmooth field. In other words, the forecast variance is lower than the\nvariance of the original observed field.\nHowever, certain applications demand that the forecast retain the same\nstatistical properties of the observations. In such cases, the S-PROG\nforecasts are of limited use and a stochatic approach might be of more\ninterest.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stochastic nowcast with STEPS\n\nThe S-PROG approach is extended to include a stochastic term which represents\nthe variance associated to the unpredictable development of precipitation. This\napproach is known as STEPS (short-term ensemble prediction system).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# The STEPS nowcast\nnowcast_method = nowcasts.get_method(\"steps\")\nR_f = nowcast_method(\n R[-3:, :, :],\n V,\n n_leadtimes,\n n_ens_members,\n n_cascade_levels=6,\n R_thr=-10.0,\n kmperpixel=2.0,\n timestep=timestep,\n noise_method=\"nonparametric\",\n vel_pert_method=\"bps\",\n mask_method=\"incremental\",\n seed=seed,\n)\n\n# Back-transform to rain rates\nR_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]\n\n\n# Plot the ensemble mean\nR_f_mean = np.mean(R_f[:, -1, :, :], axis=0)\nplot_precip_field(\n R_f_mean,\n geodata=metadata,\n title=\"Ensemble mean (+ %i min)\" % (n_leadtimes * timestep),\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mean of the ensemble displays similar properties as the S-PROG\nforecast seen above, although the degree of smoothing also depends on\nthe ensemble size. In this sense, the S-PROG forecast can be seen as\nthe mean of an ensemble of infinite size.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plot some of the realizations\nfig = plt.figure()\nfor i in range(4):\n ax = fig.add_subplot(221 + i)\n ax = plot_precip_field(\n R_f[i, -1, :, :], geodata=metadata, colorbar=False, axis=\"off\"\n )\n ax.set_title(\"Member %02d\" % i)\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from these two members of the ensemble, the stochastic forecast\nmantains the same variance as in the observed rainfall field.\nSTEPS also includes a stochatic perturbation of the motion field in order\nto quantify the its uncertainty.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, it is possible to derive probabilities from our ensemble forecast.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Compute exceedence probabilities for a 0.5 mm/h threshold\nP = excprob(R_f[:, -1, :, :], 0.5)\n\n# Plot the field of probabilities\nplot_precip_field(\n P,\n geodata=metadata,\n ptype=\"prob\",\n units=\"mm/h\",\n probthr=0.5,\n title=\"Exceedence probability (+ %i min)\" % (n_leadtimes * timestep),\n)\nplt.show()\n\n# sphinx_gallery_thumbnail_number = 5" ] } ], "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 }