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.
0 commit comments