|
4 | 4 | ===================== |
5 | 5 |
|
6 | 6 | This tutorial shows how to compute and plot an extrapolation nowcast using |
7 | | -Finnish radar data. |
| 7 | +MeteoSwiss radar data. |
8 | 8 |
|
9 | 9 | """ |
10 | 10 |
|
11 | | -from pylab import * |
12 | 11 | from datetime import datetime |
| 12 | +import matplotlib.pyplot as plt |
| 13 | +import numpy as np |
13 | 14 | from pprint import pprint |
14 | 15 | from pysteps import io, nowcasts, rcparams, verification |
15 | 16 | from pysteps.motion.lucaskanade import dense_lucaskanade |
16 | 17 | from pysteps.postprocessing import ensemblestats |
17 | 18 | from pysteps.utils import conversion, dimension, transformation |
18 | 19 | from pysteps.visualization import plot_precip_field |
19 | 20 |
|
20 | | -# Set nowcast parameters |
21 | | -n_ens_members = 20 |
22 | | -n_leadtimes = 6 |
23 | | -seed = 24 |
24 | 21 |
|
25 | 22 | ############################################################################### |
26 | 23 | # Read precipitation field |
27 | 24 | # ------------------------ |
28 | 25 | # |
29 | | -# First thing, the sequence of Swiss radar composites is imported, converted and |
30 | | -# transformed into units of dBR. |
| 26 | +# First, we will import the sequence of MeteoSwiss ("mch") radar composites. |
| 27 | +# You need the pysteps-data archive downloaded and the pystepsrc file |
| 28 | +# configured with the data_source paths pointing to data folders. |
31 | 29 |
|
| 30 | +# Selected case |
32 | 31 | date = datetime.strptime("201607112100", "%Y%m%d%H%M") |
33 | | -data_source = "mch" |
| 32 | +data_source = rcparams.data_sources["mch"] |
| 33 | +n_ens_members = 20 |
| 34 | +n_leadtimes = 6 |
| 35 | +seed = 24 |
| 36 | + |
| 37 | +############################################################################### |
| 38 | +# Load the data from the archive |
| 39 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
34 | 40 |
|
35 | | -# Load data source config |
36 | | -root_path = rcparams.data_sources[data_source]["root_path"] |
37 | | -path_fmt = rcparams.data_sources[data_source]["path_fmt"] |
38 | | -fn_pattern = rcparams.data_sources[data_source]["fn_pattern"] |
39 | | -fn_ext = rcparams.data_sources[data_source]["fn_ext"] |
40 | | -importer_name = rcparams.data_sources[data_source]["importer"] |
41 | | -importer_kwargs = rcparams.data_sources[data_source]["importer_kwargs"] |
42 | | -timestep = rcparams.data_sources[data_source]["timestep"] |
| 41 | +root_path = data_source["root_path"] |
| 42 | +path_fmt = data_source["path_fmt"] |
| 43 | +fn_pattern = data_source["fn_pattern"] |
| 44 | +fn_ext = data_source["fn_ext"] |
| 45 | +importer_name = data_source["importer"] |
| 46 | +importer_kwargs = data_source["importer_kwargs"] |
| 47 | +timestep = data_source["timestep"] |
43 | 48 |
|
44 | 49 | # Find the radar files in the archive |
45 | 50 | fns = io.find_by_date( |
|
58 | 63 |
|
59 | 64 | # Plot the rainfall field |
60 | 65 | plot_precip_field(R[-1, :, :], geodata=metadata) |
| 66 | +plt.show() |
61 | 67 |
|
62 | 68 | # Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h, |
63 | 69 | # set the fill value to -15 dBR |
|
101 | 107 | R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0] |
102 | 108 |
|
103 | 109 | # Plot some of the realizations |
104 | | -fig = figure() |
| 110 | +fig = plt.figure() |
105 | 111 | for i in range(4): |
106 | 112 | ax = fig.add_subplot(221 + i) |
107 | 113 | ax.set_title("Member %02d" % i) |
108 | 114 | plot_precip_field(R_f[i, -1, :, :], geodata=metadata, colorbar=False, axis="off") |
109 | | -tight_layout() |
| 115 | +plt.tight_layout() |
| 116 | +plt.show() |
110 | 117 |
|
111 | 118 | ############################################################################### |
112 | 119 | # Verification |
|
145 | 152 | # compute the exceedance probability of 0.1 mm/h from the ensemble |
146 | 153 | P_f = ensemblestats.excprob(R_f[:, -1, :, :], 0.1, ignore_nan=True) |
147 | 154 |
|
148 | | -# compute and plot the ROC curve |
| 155 | +############################################################################### |
| 156 | +# ROC curve |
| 157 | +# ~~~~~~~~~ |
| 158 | + |
149 | 159 | roc = verification.ROC_curve_init(0.1, n_prob_thrs=10) |
150 | 160 | verification.ROC_curve_accum(roc, P_f, R_o[-1, :, :]) |
151 | | -fig, ax = subplots() |
| 161 | +fig, ax = plt.subplots() |
152 | 162 | verification.plot_ROC(roc, ax, opt_prob_thr=True) |
153 | 163 | ax.set_title("ROC curve (+ %i min)" % (n_leadtimes * timestep)) |
| 164 | +plt.show() |
| 165 | + |
| 166 | +############################################################################### |
| 167 | +# Reliability diagram |
| 168 | +# ~~~~~~~~~~~~~~~~~~~ |
154 | 169 |
|
155 | | -# compute and plot the reliability diagram |
156 | 170 | reldiag = verification.reldiag_init(0.1) |
157 | 171 | verification.reldiag_accum(reldiag, P_f, R_o[-1, :, :]) |
158 | | -fig, ax = subplots() |
| 172 | +fig, ax = plt.subplots() |
159 | 173 | verification.plot_reldiag(reldiag, ax) |
160 | 174 | ax.set_title("Reliability diagram (+ %i min)" % (n_leadtimes * timestep)) |
| 175 | +plt.show() |
| 176 | + |
| 177 | +############################################################################### |
| 178 | +# Rank histogram |
| 179 | +# ~~~~~~~~~~~~~~ |
161 | 180 |
|
162 | | -# compute and plot the rank histogram |
163 | 181 | rankhist = verification.rankhist_init(R_f.shape[0], 0.1) |
164 | 182 | verification.rankhist_accum(rankhist, R_f[:, -1, :, :], R_o[-1, :, :]) |
165 | | -fig, ax = subplots() |
| 183 | +fig, ax = plt.subplots() |
166 | 184 | verification.plot_rankhist(rankhist, ax) |
167 | 185 | ax.set_title("Rank histogram (+ %i min)" % (n_leadtimes * timestep)) |
| 186 | +plt.show() |
168 | 187 |
|
169 | 188 | # sphinx_gallery_thumbnail_number = 5 |
0 commit comments