Skip to content

Commit 11dbae5

Browse files
pulkkinspulkkins
authored
Anvil (#150)
1 parent 4382e96 commit 11dbae5

File tree

8 files changed

+696
-8
lines changed

8 files changed

+696
-8
lines changed

‎doc/source/pysteps_reference/nowcasts.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ Implementation of deterministic and ensemble nowcasting methods.
66

77

88
.. automodule:: pysteps.nowcasts.interface
9+
.. automodule:: pysteps.nowcasts.anvil
910
.. automodule:: pysteps.nowcasts.extrapolation
1011
.. automodule:: pysteps.nowcasts.sprog
1112
.. automodule:: pysteps.nowcasts.sseps

‎doc/source/references.bib‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,15 @@ @ARTICLE{PCH2019b
116116
YEAR = 2019
117117
}
118118

119+
@ARTICLE{PCLH2020,
120+
AUTHOR = "S. Pulkkinen and V. Chandrasekar and A. von Lerber and A.-M. Harri",
121+
TITLE = "Nowcasting of Convective Rainfall Using Volumetric Radar Observations",
122+
JOURNAL = "IEEE Transactions on Geoscience and Remote Sensing",
123+
DOI = "10.1109/TGRS.2020.2984594",
124+
PAGES = "1--15",
125+
YEAR = 2020
126+
}
127+
119128
@INCOLLECTION{PGPO1994,
120129
AUTHOR = "M. Proesmans and L. van Gool and E. Pauwels and A. Oosterlinck",
121130
TITLE = "Determination of optical flow and its discontinuities using non-linear diffusion",

‎examples/anvil_nowcast.py‎

Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
1+
# coding: utf-8
2+
3+
"""
4+
ANVIL nowcast
5+
=============
6+
7+
This example demonstrates how to use ANVIL and the advantages compared to
8+
extrapolation nowcast and S-PROG.
9+
10+
Load the libraries.
11+
"""
12+
from datetime import datetime, timedelta
13+
import warnings
14+
warnings.simplefilter("ignore")
15+
import matplotlib.pyplot as plt
16+
import numpy as np
17+
from pysteps import motion, io, rcparams, utils
18+
from pysteps.nowcasts import anvil, extrapolation, sprog
19+
from pysteps.visualization import plot_precip_field
20+
21+
###############################################################################
22+
# Read the input data
23+
# -------------------
24+
#
25+
# ANVIL was originally developed to use vertically integrated liquid (VIL) as
26+
# the input data, but the model allows using any two-dimensional input fields.
27+
# Here we use a composite of rain rates.
28+
29+
date = datetime.strptime("201505151620", "%Y%m%d%H%M")
30+
31+
# Read the data source information from rcparams
32+
data_source = rcparams.data_sources["mch"]
33+
34+
root_path = data_source["root_path"]
35+
path_fmt = data_source["path_fmt"]
36+
fn_pattern = data_source["fn_pattern"]
37+
fn_ext = data_source["fn_ext"]
38+
importer_name = data_source["importer"]
39+
importer_kwargs = data_source["importer_kwargs"]
40+
41+
# Find the input files in the archive. Use history length of 5 timesteps
42+
filenames = io.archive.find_by_date(date, root_path, path_fmt, fn_pattern,
43+
fn_ext, timestep=5, num_prev_files=5)
44+
45+
# Read the input time series
46+
importer = io.get_method(importer_name, "importer")
47+
rainrate_field, quality, metadata = io.read_timeseries(filenames, importer,
48+
**importer_kwargs)
49+
50+
# Convert to rain rate (mm/h)
51+
rainrate_field, metadata = utils.to_rainrate(rainrate_field, metadata)
52+
53+
################################################################################
54+
# Compute the advection field
55+
# ---------------------------
56+
#
57+
# Apply the Lucas-Kanade method with the parameters given in Pulkkinen et al.
58+
# (2020) to compute the advection field.
59+
60+
fd_kwargs = {}
61+
fd_kwargs["max_corners"] = 1000
62+
fd_kwargs["quality_level"] = 0.01
63+
fd_kwargs["min_distance"] = 2
64+
fd_kwargs["block_size"] = 8
65+
66+
lk_kwargs = {}
67+
lk_kwargs["winsize"] = (15, 15)
68+
69+
oflow_kwargs = {}
70+
oflow_kwargs["fd_kwargs"] = fd_kwargs
71+
oflow_kwargs["lk_kwargs"] = lk_kwargs
72+
oflow_kwargs["decl_scale"] = 10
73+
74+
oflow = motion.get_method("lucaskanade")
75+
76+
# transform the input data to logarithmic scale
77+
rainrate_field_log, _ = utils.transformation.dB_transform(rainrate_field,
78+
metadata=metadata)
79+
velocity = oflow(rainrate_field_log, **oflow_kwargs)
80+
81+
###############################################################################
82+
# Compute the nowcasts and threshold rain rates below 0.5 mm/h
83+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
84+
forecast_extrap = extrapolation.forecast(rainrate_field[-1], velocity, 3,
85+
extrap_kwargs={"allow_nonfinite_values": True})
86+
forecast_extrap[forecast_extrap < 0.5] = 0.0
87+
88+
rainrate_field_finite = rainrate_field.copy()
89+
rainrate_field_finite[~np.isfinite(rainrate_field_finite)] = 0.0
90+
forecast_sprog = sprog.forecast(rainrate_field_finite[-3:], velocity, 3,
91+
n_cascade_levels=8, R_thr=0.5)
92+
forecast_sprog[~np.isfinite(forecast_extrap)] = np.nan
93+
forecast_sprog[forecast_sprog < 0.5] = 0.0
94+
95+
forecast_anvil = anvil.forecast(rainrate_field[-4:], None, velocity, 3,
96+
ar_window_radius=25, ar_order=2)
97+
forecast_anvil[forecast_anvil < 0.5] = 0.0
98+
99+
###############################################################################
100+
# Read the reference observation field and threshold rain rates below 0.5 mm/h
101+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
102+
filenames = io.archive.find_by_date(date, root_path, path_fmt, fn_pattern,
103+
fn_ext, timestep=5, num_next_files=3)
104+
105+
refobs_field, quality, metadata = io.read_timeseries(filenames, importer,
106+
**importer_kwargs)
107+
108+
refobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata)
109+
refobs_field[refobs_field < 0.5] = 0.0
110+
111+
###############################################################################
112+
# Plot the extrapolation, S-PROG and ANVIL nowcasts.
113+
# --------------------------------------------------
114+
#
115+
# For comparison, the observed rain rate fields are also plotted. Growth and
116+
# decay areas are marked with red and blue circles, respectively.
117+
def plot_growth_decay_circles(ax):
118+
circle = plt.Circle((360, 300), 25, color="b", clip_on=False, fill=False,
119+
zorder=1e9)
120+
ax.add_artist(circle)
121+
circle = plt.Circle((420, 350), 30, color="b", clip_on=False, fill=False,
122+
zorder=1e9)
123+
ax.add_artist(circle)
124+
circle = plt.Circle((405, 380), 30, color="b", clip_on=False, fill=False,
125+
zorder=1e9)
126+
ax.add_artist(circle)
127+
circle = plt.Circle((420, 500), 25, color="b", clip_on=False, fill=False,
128+
zorder=1e9)
129+
ax.add_artist(circle)
130+
circle = plt.Circle((480, 535), 30, color="b", clip_on=False, fill=False,
131+
zorder=1e9)
132+
ax.add_artist(circle)
133+
circle = plt.Circle((330, 470), 35, color="b", clip_on=False, fill=False,
134+
zorder=1e9)
135+
ax.add_artist(circle)
136+
circle = plt.Circle((505, 205), 30, color="b", clip_on=False, fill=False,
137+
zorder=1e9)
138+
ax.add_artist(circle)
139+
circle = plt.Circle((440, 180), 30, color="r", clip_on=False, fill=False,
140+
zorder=1e9)
141+
ax.add_artist(circle)
142+
circle = plt.Circle((590, 240), 30, color="r", clip_on=False, fill=False,
143+
zorder=1e9)
144+
ax.add_artist(circle)
145+
circle = plt.Circle((585, 160), 15, color="r", clip_on=False, fill=False,
146+
zorder=1e9)
147+
ax.add_artist(circle)
148+
149+
fig = plt.figure(figsize=(10, 13))
150+
151+
ax = fig.add_subplot(321)
152+
rainrate_field[-1][rainrate_field[-1] < 0.5] = 0.0
153+
plot_precip_field(rainrate_field[-1])
154+
plot_growth_decay_circles(ax)
155+
ax.set_title("Obs. %s" % str(date))
156+
157+
ax = fig.add_subplot(322)
158+
plot_precip_field(refobs_field)
159+
plot_growth_decay_circles(ax)
160+
ax.set_title("Obs. %s" % str(date + timedelta(minutes=15)))
161+
162+
ax = fig.add_subplot(323)
163+
plot_precip_field(forecast_extrap[-1])
164+
plot_growth_decay_circles(ax)
165+
ax.set_title("Extrapolation +15 minutes")
166+
167+
ax = fig.add_subplot(324)
168+
plot_precip_field(forecast_sprog[-1])
169+
plot_growth_decay_circles(ax)
170+
ax.set_title("S-PROG (with post-processing)\n +15 minutes")
171+
172+
ax = fig.add_subplot(325)
173+
plot_precip_field(forecast_anvil[-1])
174+
plot_growth_decay_circles(ax)
175+
ax.set_title("ANVIL +15 minutes")
176+
177+
plt.show()
178+
179+
###############################################################################
180+
# Remarks
181+
# -------
182+
#
183+
# The extrapolation nowcast is static, i.e. it does not predict any growth or
184+
# decay. While S-PROG is to some extent able to predict growth and decay, this
185+
# this comes with loss of small-scale features. In addition, statistical
186+
# post-processing needs to be applied to correct the bias and incorrect wet-area
187+
# ratio introduced by the autoregressive process. ANVIL is able to do both:
188+
# predict growth and decay and preserve the small-scale structure in a way that
189+
# post-processing is not necessary.

‎pysteps/cascade/decomposition.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -259,7 +259,8 @@ def recompose_fft(decomp, **kwargs):
259259
mu = decomp["means"]
260260
sigma = decomp["stds"]
261261

262-
if not decomp["normalized"] and not decomp["compact_output"]:
262+
if not decomp["normalized"] and not \
263+
(decomp["domain"] == "spectral" and decomp["compact_output"]):
263264
return np.sum(levels, axis=0)
264265
else:
265266
if decomp["compact_output"]:

0 commit comments

Comments
 (0)