|
10 | 10 |
|
11 | 11 | from pylab import *
|
12 | 12 | from datetime import datetime
|
| 13 | +from pprint import pprint |
13 | 14 | from pysteps import io, nowcasts, rcparams
|
14 | 15 | from pysteps.motion.lucaskanade import dense_lucaskanade
|
15 | 16 | from pysteps.postprocessing.ensemblestats import excprob
|
|
29 | 30 | # transformed into units of dBR.
|
30 | 31 |
|
31 | 32 |
|
32 |
| - |
33 | 33 | date = datetime.strptime("201701311200", "%Y%m%d%H%M")
|
34 | 34 | data_source = "mch"
|
35 | 35 |
|
|
44 | 44 |
|
45 | 45 | # Find the radar files in the archive
|
46 | 46 | fns = io.find_by_date(
|
47 |
| - date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2, |
| 47 | + date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2 |
48 | 48 | )
|
49 | 49 |
|
50 | 50 | # Read the data from the archive
|
|
60 | 60 | # Plot the rainfall field
|
61 | 61 | plot_precip_field(R[-1, :, :], geodata=metadata)
|
62 | 62 |
|
63 |
| -# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h |
64 |
| -R = transformation.dB_transform(R, threshold=0.1, zerovalue=-15.0)[0] |
| 63 | +# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h, |
| 64 | +# set the fill value to -15 dBR |
| 65 | +R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0) |
65 | 66 |
|
66 | 67 | # Set missing values with the fill value
|
67 | 68 | R[~np.isfinite(R)] = -15.0
|
68 | 69 |
|
| 70 | +# Nicely print the metadata |
| 71 | +pprint(metadata) |
| 72 | + |
69 | 73 | ###############################################################################
|
70 | 74 | # Deterministic nowcast with S-PROG
|
71 | 75 | # ---------------------------------
|
|
96 | 100 | R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]
|
97 | 101 |
|
98 | 102 | # Plot the S-PROG forecast
|
99 |
| -plot_precip_field(R_f[-1, :, :], geodata=metadata, title="S-PROG") |
| 103 | +plot_precip_field( |
| 104 | + R_f[-1, :, :], |
| 105 | + geodata=metadata, |
| 106 | + title="S-PROG (+ %i min)" % (n_leadtimes * timestep), |
| 107 | +) |
100 | 108 |
|
101 | 109 | ###############################################################################
|
102 | 110 | # As we can see from the figure above, the forecast produced by S-PROG is a
|
103 | 111 | # smooth field. In other words, the forecast variance is lower than the
|
104 | 112 | # variance of the original observed field.
|
105 |
| -# However, given applications demand that the forecast retain the same |
| 113 | +# However, certain applications demand that the forecast retain the same |
106 | 114 | # statistical properties of the observations. In such cases, the S-PROG
|
107 |
| -# forecsats are of limited use and a stochatic approahc might be of more |
| 115 | +# forecasts are of limited use and a stochatic approach might be of more |
108 | 116 | # interest.
|
109 | 117 |
|
110 | 118 | ###############################################################################
|
111 | 119 | # Stochastic nowcast with STEPS
|
112 | 120 | # -----------------------------
|
113 | 121 | #
|
114 |
| -# |
115 | 122 | # The S-PROG approach is extended to include a stochastic term which represents
|
116 |
| -# the variance linked to the unpredictable development of precipitation. This |
| 123 | +# the variance associated to the unpredictable development of precipitation. This |
117 | 124 | # approach is known as STEPS (short-term ensemble prediction system).
|
118 | 125 |
|
119 | 126 | # The STEPES nowcast
|
|
141 | 148 |
|
142 | 149 | # Plot the ensemble mean
|
143 | 150 | R_f_mean = np.mean(R_f[:, -1, :, :], axis=0)
|
144 |
| -plot_precip_field(R_f_mean, geodata=metadata, title="Ensemble mean") |
| 151 | +plot_precip_field( |
| 152 | + R_f_mean, |
| 153 | + geodata=metadata, |
| 154 | + title="Ensemble mean (+ %i min)" % (n_leadtimes * timestep), |
| 155 | +) |
145 | 156 |
|
146 | 157 | ###############################################################################
|
147 | 158 | # The mean of the ensemble displays similar properties as the S-PROG
|
148 |
| -# forecast seen above, although the degree of smoothing strongly depends on |
| 159 | +# forecast seen above, although the degree of smoothing also depends on |
149 | 160 | # the ensemble size. In this sense, the S-PROG forecast can be seen as
|
150 | 161 | # the mean of an ensemble of infinite size.
|
151 | 162 |
|
152 | 163 | # Plot the first two realizations
|
153 | 164 | fig = figure()
|
154 |
| -for i in range(2): |
155 |
| - ax = fig.add_subplot(121 + i) |
| 165 | +for i in range(4): |
| 166 | + ax = fig.add_subplot(221 + i) |
156 | 167 | ax.set_title("Member %02d" % i)
|
157 | 168 | plot_precip_field(R_f[i, -1, :, :], geodata=metadata, colorbar=False, axis="off")
|
158 | 169 | tight_layout()
|
159 | 170 |
|
160 | 171 | ###############################################################################
|
161 | 172 | # As we can see from these two members of the ensemble, the stochastic forecast
|
162 | 173 | # mantains the same variance as in the observed rainfall field.
|
| 174 | +# STEPS also includes a stochatic perturbation of the motion field in order |
| 175 | +# to quantify the its uncertainty. |
163 | 176 |
|
164 | 177 | ###############################################################################
|
165 | 178 | # Finally, it is possible to derive probabilities from our ensemble forecast.
|
|
175 | 188 | type="prob",
|
176 | 189 | units="mm/h",
|
177 | 190 | probthr=0.5,
|
178 |
| - title="Exceedence probability", |
| 191 | + title="Exceedence probability (+ %i min)" % (n_leadtimes * timestep), |
179 | 192 | )
|
180 | 193 |
|
181 | 194 | # sphinx_gallery_thumbnail_number = 5
|
0 commit comments