|
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