Note
Go to the end to download the full example code.
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).
import matplotlib.pyplot as plt
import numpy as np
from pysteps.nowcasts.lagrangian_probability import forecast
from pysteps.visualization import plot_precip_field
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.
# 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()
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.
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()
Let’s plot one single leadtime in more detail using the pysteps visualization functionality.
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()
/home/docs/checkouts/readthedocs.org/user_builds/pysteps/envs/latest/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(
Verification#
# 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()
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.
# sphinx_gallery_thumbnail_number = 3
Total running time of the script: (0 minutes 2.988 seconds)