{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Ensemble verification\n\nIn this tutorial we perform a verification of a probabilistic extrapolation nowcast \nusing MeteoSwiss radar data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from datetime import datetime\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pprint import pprint\nfrom pysteps import io, nowcasts, rcparams, verification\nfrom pysteps.motion.lucaskanade import dense_lucaskanade\nfrom pysteps.postprocessing import ensemblestats\nfrom pysteps.utils import conversion, dimension, transformation\nfrom pysteps.visualization import plot_precip_field" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read precipitation field\n\nFirst, we will import the sequence of MeteoSwiss (\"mch\") radar composites.\nYou need the pysteps-data archive downloaded and the pystepsrc file\nconfigured with the data_source paths pointing to data folders.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Selected case\ndate = datetime.strptime(\"201607112100\", \"%Y%m%d%H%M\")\ndata_source = rcparams.data_sources[\"mch\"]\nn_ens_members = 20\nn_leadtimes = 6\nseed = 24" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the data from the archive\n\nThe data are upscaled to 2 km resolution to limit the memory usage and thus\nbe able to afford a larger number of ensemble members.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "root_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# 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\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": [ "## Forecast\n\nWe use the STEPS approach to produce a ensemble nowcast of precipitation fields.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Estimate the motion field\nV = dense_lucaskanade(R)\n\n# Perform the ensemble nowcast with STEPS\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 decomp_method=\"fft\",\n bandpass_filter_method=\"gaussian\",\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# Plot some of the realizations\nfig = plt.figure()\nfor i in range(4):\n ax = fig.add_subplot(221 + i)\n ax.set_title(\"Member %02d\" % i)\n plot_precip_field(R_f[i, -1, :, :], geodata=metadata, colorbar=False, axis=\"off\")\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verification\n\nPysteps includes a number of verification metrics to help users to analyze\nthe general characteristics of the nowcasts in terms of consistency and\nquality (or goodness).\nHere, we will verify our probabilistic forecasts using the ROC curve,\nreliability diagrams, and rank histograms, as implemented in the verification\nmodule of pysteps.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Find the files containing the verifying observations\nfns = io.archive.find_by_date(\n date,\n root_path,\n path_fmt,\n fn_pattern,\n fn_ext,\n timestep,\n 0,\n num_next_files=n_leadtimes,\n)\n\n# Read the observations\nR_o, _, metadata_o = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Convert to mm/h\nR_o, metadata_o = conversion.to_rainrate(R_o, metadata_o)\n\n# Upscale data to 2 km\nR_o, metadata_o = dimension.aggregate_fields_space(R_o, metadata_o, 2000)\n\n# Compute the verification for the last lead time\n\n# compute the exceedance probability of 0.1 mm/h from the ensemble\nP_f = ensemblestats.excprob(R_f[:, -1, :, :], 0.1, ignore_nan=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ROC curve\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "roc = verification.ROC_curve_init(0.1, n_prob_thrs=10)\nverification.ROC_curve_accum(roc, P_f, R_o[-1, :, :])\nfig, ax = plt.subplots()\nverification.plot_ROC(roc, ax, opt_prob_thr=True)\nax.set_title(\"ROC curve (+%i min)\" % (n_leadtimes * timestep))\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reliability diagram\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reldiag = verification.reldiag_init(0.1)\nverification.reldiag_accum(reldiag, P_f, R_o[-1, :, :])\nfig, ax = plt.subplots()\nverification.plot_reldiag(reldiag, ax)\nax.set_title(\"Reliability diagram (+%i min)\" % (n_leadtimes * timestep))\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rank histogram\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rankhist = verification.rankhist_init(R_f.shape[0], 0.1)\nverification.rankhist_accum(rankhist, R_f[:, -1, :, :], R_o[-1, :, :])\nfig, ax = plt.subplots()\nverification.plot_rankhist(rankhist, ax)\nax.set_title(\"Rank histogram (+%i min)\" % (n_leadtimes * timestep))\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 }