.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/probability_forecast.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_probability_forecast.py: Probability forecasts ===================== This example script shows how to forecast the probability of exceeding an intensity threshold. The method is based on the local Lagrangian approach described in Germann and Zawadzki (2004). .. GENERATED FROM PYTHON SOURCE LINES 12-19 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from pysteps.nowcasts.lagrangian_probability import forecast from pysteps.visualization import plot_precip_field .. GENERATED FROM PYTHON SOURCE LINES 20-28 Numerical example ----------------- First, we use some dummy data to show the basic principle of this approach. The probability forecast is produced by sampling a spatial neighborhood that is increased as a function of lead time. As a result, the edges of the yellow square becomes more and more smooth as t increases. This represents the strong loss of predictability with lead time of any extrapolation nowcast. .. GENERATED FROM PYTHON SOURCE LINES 28-48 .. code-block:: Python # parameters precip = np.zeros((100, 100)) precip[10:50, 10:50] = 1 velocity = np.ones((2, *precip.shape)) timesteps = [0, 2, 6, 12] thr = 0.5 slope = 1 # pixels / timestep # compute probability forecast out = forecast(precip, velocity, timesteps, thr, slope=slope) # plot for n, frame in enumerate(out): plt.subplot(2, 2, n + 1) plt.imshow(frame, interpolation="nearest", vmin=0, vmax=1) plt.title(f"t={timesteps[n]}") plt.xticks([]) plt.yticks([]) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_probability_forecast_001.png :alt: t=0, t=2, t=6, t=12 :srcset: /auto_examples/images/sphx_glr_probability_forecast_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 49-56 Real-data example ----------------- We now apply the same method to real data. We use a slope of 1 km / minute as suggested by Germann and Zawadzki (2004), meaning that after 30 minutes, the probabilities are computed by using all pixels within a neighborhood of 30 kilometers. .. GENERATED FROM PYTHON SOURCE LINES 56-103 .. code-block:: Python from datetime import datetime from pysteps import io, rcparams, utils from pysteps.motion.lucaskanade import dense_lucaskanade from pysteps.verification import reldiag_init, reldiag_accum, plot_reldiag # data source source = rcparams.data_sources["mch"] root = rcparams.data_sources["mch"]["root_path"] fmt = rcparams.data_sources["mch"]["path_fmt"] pattern = rcparams.data_sources["mch"]["fn_pattern"] ext = rcparams.data_sources["mch"]["fn_ext"] timestep = rcparams.data_sources["mch"]["timestep"] importer_name = rcparams.data_sources["mch"]["importer"] importer_kwargs = rcparams.data_sources["mch"]["importer_kwargs"] # read precip field date = datetime.strptime("201607112100", "%Y%m%d%H%M") fns = io.find_by_date(date, root, fmt, pattern, ext, timestep, num_prev_files=2) importer = io.get_method(importer_name, "importer") precip, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) precip, metadata = utils.to_rainrate(precip, metadata) # precip[np.isnan(precip)] = 0 # motion motion = dense_lucaskanade(precip) # parameters nleadtimes = 6 thr = 1 # mm / h slope = 1 * timestep # km / min # compute probability forecast extrap_kwargs = dict(allow_nonfinite_values=True) fct = forecast( precip[-1], motion, nleadtimes, thr, slope=slope, extrap_kwargs=extrap_kwargs ) # plot for n, frame in enumerate(fct): plt.subplot(2, 3, n + 1) plt.imshow(frame, interpolation="nearest", vmin=0, vmax=1) plt.xticks([]) plt.yticks([]) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_probability_forecast_002.png :alt: probability forecast :srcset: /auto_examples/images/sphx_glr_probability_forecast_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 104-106 Let's plot one single leadtime in more detail using the pysteps visualization functionality. .. GENERATED FROM PYTHON SOURCE LINES 106-118 .. code-block:: Python plt.close() # Plot the field of probabilities plot_precip_field( fct[2], geodata=metadata, ptype="prob", probthr=thr, title="Exceedence probability (+ %i min)" % (nleadtimes * timestep), ) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_probability_forecast_003.png :alt: Exceedence probability (+ 30 min) :srcset: /auto_examples/images/sphx_glr_probability_forecast_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/pysteps/envs/v1.17.0/lib/python3.10/site-packages/pysteps/visualization/utils.py:439: UserWarning: cartopy package is required for the get_geogrid function but it is not installed. Ignoring geographical information. warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 119-121 Verification ------------ .. GENERATED FROM PYTHON SOURCE LINES 121-140 .. code-block:: Python # verifying observations importer = io.get_method(importer_name, "importer") fns = io.find_by_date( date, root, fmt, pattern, ext, timestep, num_next_files=nleadtimes ) obs, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) obs, metadata = utils.to_rainrate(obs, metadata) obs[np.isnan(obs)] = 0 # reliability diagram reldiag = reldiag_init(thr) reldiag_accum(reldiag, fct, obs[1:]) fig, ax = plt.subplots() plot_reldiag(reldiag, ax) ax.set_title("Reliability diagram") plt.show() .. image-sg:: /auto_examples/images/sphx_glr_probability_forecast_004.png :alt: Reliability diagram :srcset: /auto_examples/images/sphx_glr_probability_forecast_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-147 References ---------- Germann, U. and I. Zawadzki, 2004: Scale Dependence of the Predictability of Precipitation from Continental Radar Images. Part II: Probability Forecasts. Journal of Applied Meteorology, 43(1), 74-89. .. GENERATED FROM PYTHON SOURCE LINES 147-149 .. code-block:: Python # sphinx_gallery_thumbnail_number = 3 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.920 seconds) .. _sphx_glr_download_auto_examples_probability_forecast.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: probability_forecast.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: probability_forecast.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: probability_forecast.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_