|
11 | 11 |
|
12 | 12 | # global imports |
13 | 13 | import logging |
14 | | -logging.basicConfig(level=logging.INFO) |
15 | 14 | import numpy as np |
16 | 15 | import copy |
17 | 16 | import datetime |
18 | 17 | import matplotlib.pyplot as plt |
19 | 18 | import pandas as pd |
20 | 19 |
|
21 | 20 | # Local imports |
22 | | -from ..common.utils import read_df, get_qpe_files, get_qpe_files_multiple_dirs, perfscores |
| 21 | +from ..common.utils import ( |
| 22 | + read_df, |
| 23 | + get_qpe_files, |
| 24 | + get_qpe_files_multiple_dirs, |
| 25 | + perfscores, |
| 26 | +) |
23 | 27 | from ..common.utils import timestamp_from_datetime, nearest_time |
24 | 28 | from ..common.lookup import get_lookup |
25 | 29 | from ..common.io_data import read_cart |
26 | 30 | from ..common.graphics import score_plot, qpe_scatterplot |
27 | 31 | from ..common.retrieve_data import retrieve_CPCCV |
28 | 32 |
|
29 | | -def evaluation(qpefolder, gaugepattern, list_models = None, |
30 | | - outputfolder = './', t0 = None, t1 = None, |
31 | | - bounds10 = [0,2,10,100], bounds60 = [0,1,10,100]): |
32 | | - |
| 33 | +logging.basicConfig(level=logging.INFO) |
| 34 | + |
| 35 | + |
| 36 | +def evaluation( |
| 37 | + qpefolder, |
| 38 | + gaugepattern, |
| 39 | + list_models=None, |
| 40 | + outputfolder="./", |
| 41 | + t0=None, |
| 42 | + t1=None, |
| 43 | + bounds10=[0, 2, 10, 100], |
| 44 | + bounds60=[0, 1, 10, 100], |
| 45 | +): |
33 | 46 | """ |
34 | | - PErforms an evaluation of QPE products with reference gauge data |
35 | | - |
36 | | - Parameters |
37 | | - ---------- |
38 | | - qpefolder : str |
39 | | - Main directory where the QPE data is stored, each model corresponding |
40 | | - to a subfolder, as given by qpe_compute.py |
41 | | - gaugepattern : str |
42 | | - The pattern of gauge files that contain the gauge data. |
43 | | - on CSCS: '/store/msrad/radar/radar_database/gauge/*.csv.gz' |
44 | | - t0: datetime.datetime instance |
45 | | - starting time of the time range, default is first timestep available |
46 | | - t1 : datetime.datetime instance |
47 | | - end time of the time range, default is last timestep available |
48 | | - bounds10 : list of float |
49 | | - list of precipitation bounds for which to compute scores separately |
50 | | - at 10 min time resolution |
51 | | - [0,2,10,100] will give scores in range [0-2], [2-10] and [10-100] |
52 | | - bounds60 : list of float |
53 | | - list of precipitation bounds for which to compute scores separately |
54 | | - at hourly time resolution |
55 | | - [0,1,10,100] will give scores in range [0-1], [1-10] and [10-100] |
56 | | - list_models : list of str |
57 | | - list of models to use in the evaluation, default is to use all |
58 | | - subfolders (models) available in qpefolder |
59 | | - |
60 | | - |
| 47 | + PErforms an evaluation of QPE products with reference gauge data |
| 48 | +
|
| 49 | + Parameters |
| 50 | + ---------- |
| 51 | + qpefolder : str |
| 52 | + Main directory where the QPE data is stored, each model corresponding |
| 53 | + to a subfolder, as given by qpe_compute.py |
| 54 | + gaugepattern : str |
| 55 | + The pattern of gauge files that contain the gauge data. |
| 56 | + on CSCS: '/store/msrad/radar/radar_database/gauge/*.csv.gz' |
| 57 | + t0: datetime.datetime instance |
| 58 | + starting time of the time range, default is first timestep available |
| 59 | + t1 : datetime.datetime instance |
| 60 | + end time of the time range, default is last timestep available |
| 61 | + bounds10 : list of float |
| 62 | + list of precipitation bounds for which to compute scores separately |
| 63 | + at 10 min time resolution |
| 64 | + [0,2,10,100] will give scores in range [0-2], [2-10] and [10-100] |
| 65 | + bounds60 : list of float |
| 66 | + list of precipitation bounds for which to compute scores separately |
| 67 | + at hourly time resolution |
| 68 | + [0,1,10,100] will give scores in range [0-1], [1-10] and [10-100] |
| 69 | + list_models : list of str |
| 70 | + list of models to use in the evaluation, default is to use all |
| 71 | + subfolders (models) available in qpefolder |
| 72 | +
|
| 73 | +
|
61 | 74 | """ |
62 | | - |
63 | | - if type(qpefolder) == list: |
64 | | - logging.info('Getting all files from multiple qpe folders') |
65 | | - tmp = get_qpe_files_multiple_dirs(qpefolder, time_agg = 10, list_models = list_models) |
| 75 | + |
| 76 | + if type(qpefolder) is list: |
| 77 | + logging.info("Getting all files from multiple qpe folders") |
| 78 | + tmp = get_qpe_files_multiple_dirs( |
| 79 | + qpefolder, time_agg=10, list_models=list_models |
| 80 | + ) |
66 | 81 | else: |
67 | | - logging.info('Getting all files from qpe folder {:s}'.format(qpefolder)) |
68 | | - tmp = get_qpe_files(qpefolder, time_agg = 10, list_models = list_models) |
69 | | - |
| 82 | + logging.info("Getting all files from qpe folder {:s}".format(qpefolder)) |
| 83 | + tmp = get_qpe_files(qpefolder, time_agg=10, list_models=list_models) |
| 84 | + |
70 | 85 | # Get only timesteps where at least 2 files are available during 10 min period |
71 | 86 | qpe_files10 = copy.deepcopy(tmp) |
72 | 87 | for k in tmp.keys(): |
73 | 88 | for m in tmp[k].keys(): |
74 | 89 | if len(tmp[k][m]) < 2: |
75 | 90 | del qpe_files10[k][m] |
76 | | - |
| 91 | + |
77 | 92 | # number of models by timestep |
78 | 93 | nmodels = np.array([len(d) for d in qpe_files10.values()]) |
79 | | - |
80 | | - |
| 94 | + |
81 | 95 | # Get only timesteps where all qpe models are available |
82 | 96 | qpe_files10_filt = {} |
83 | 97 | for i, k in enumerate(qpe_files10.keys()): |
84 | 98 | if nmodels[i] == max(nmodels): |
85 | 99 | qpe_files10_filt[k] = qpe_files10[k] |
86 | | - |
| 100 | + |
87 | 101 | models = list(list(qpe_files10_filt.values())[0].keys()) |
88 | | - if list_models == None: |
| 102 | + if list_models is None: |
89 | 103 | list_models = nmodels |
90 | | - |
| 104 | + |
91 | 105 | tsteps = sorted(list(qpe_files10_filt.keys())) |
92 | | - |
93 | | - logging.info('Reading gauge data from pattern {:s}'.format(gaugepattern)) |
| 106 | + |
| 107 | + logging.info("Reading gauge data from pattern {:s}".format(gaugepattern)) |
94 | 108 | df = read_df(gaugepattern) |
95 | | - |
96 | | - |
97 | | - logging.info('Converting to pandas dataframe...') |
| 109 | + |
| 110 | + logging.info("Converting to pandas dataframe...") |
98 | 111 | t0 = timestamp_from_datetime(tsteps[0]) |
99 | 112 | t1 = timestamp_from_datetime(tsteps[-1]) |
100 | | - df = df[(df['TIMESTAMP'] >= t0) & (df['TIMESTAMP'] <= t1)].compute() |
101 | | - stations = np.unique(df['STATION']) |
102 | | - |
103 | | - logging.info('Getting lookup table') |
104 | | - lut = get_lookup('station_to_qpegrid') |
105 | | - |
| 113 | + df = df[(df["TIMESTAMP"] >= t0) & (df["TIMESTAMP"] <= t1)].compute() |
| 114 | + stations = np.unique(df["STATION"]) |
| 115 | + |
| 116 | + logging.info("Getting lookup table") |
| 117 | + lut = get_lookup("station_to_qpegrid") |
| 118 | + |
106 | 119 | # Initialize matrices of precip at stations |
107 | 120 | precip_qpe = {} |
108 | 121 | for m in models: |
109 | | - precip_qpe[m] = np.zeros((len(tsteps), len(stations))) |
| 122 | + precip_qpe[m] = np.zeros((len(tsteps), len(stations))) |
110 | 123 | precip_ref = np.zeros((len(tsteps), len(stations))) |
111 | | - |
112 | | - for i, tstep in enumerate(tsteps): # Loop on timesteps |
113 | | - logging.info('Reading timestep {:d}/{:d}'.format(i+1, len(tsteps))) |
| 124 | + |
| 125 | + for i, tstep in enumerate(tsteps): # Loop on timesteps |
| 126 | + logging.info("Reading timestep {:d}/{:d}".format(i + 1, len(tsteps))) |
114 | 127 | # Get reference precip |
115 | 128 | tstamp = timestamp_from_datetime(tstep) |
116 | | - measures_10 = df[df['TIMESTAMP'] == tstamp] |
117 | | - idx = np.searchsorted(stations, measures_10['STATION']) # idx of stations for this timestep |
118 | | - precip_ref[i, idx] = measures_10['RRE150Z0'] * 6 |
119 | | - |
| 129 | + measures_10 = df[df["TIMESTAMP"] == tstamp] |
| 130 | + idx = np.searchsorted( |
| 131 | + stations, measures_10["STATION"] |
| 132 | + ) # idx of stations for this timestep |
| 133 | + precip_ref[i, idx] = measures_10["RRE150Z0"] * 6 |
| 134 | + |
120 | 135 | # Get QPE precip |
121 | 136 | for m in models: |
122 | 137 | for f in qpe_files10_filt[tstep][m]: |
123 | 138 | data = read_cart(f) |
124 | | - for j,s in enumerate(stations): |
125 | | - precip_qpe[m][i,j] += data[lut[s]['00'][0], lut[s]['00'][1]] |
126 | | - |
| 139 | + for j, s in enumerate(stations): |
| 140 | + precip_qpe[m][i, j] += data[lut[s]["00"][0], lut[s]["00"][1]] |
| 141 | + |
127 | 142 | precip_qpe[m][i] /= len(qpe_files10_filt[tstep][m]) |
128 | | - |
129 | | - #Get avg over 10min peri |
| 143 | + |
| 144 | + # Get avg over 10min peri |
130 | 145 | scores10 = {} |
131 | 146 | # COmpute 10min scores |
132 | 147 | valid_ref = np.isfinite(precip_ref.ravel()) |
133 | 148 | for m in models: |
134 | | - scores10[m] = perfscores(precip_qpe[m].ravel()[valid_ref], |
135 | | - precip_ref.ravel()[valid_ref], |
136 | | - bounds10) |
137 | | - |
| 149 | + scores10[m] = perfscores( |
| 150 | + precip_qpe[m].ravel()[valid_ref], precip_ref.ravel()[valid_ref], bounds10 |
| 151 | + ) |
| 152 | + |
138 | 153 | # Hourly resolution |
139 | 154 | # get hour of tsteps |
140 | 155 | hours = np.array([nearest_time(t, 60) for t in tsteps]) |
141 | | - hours_u,cnt = np.unique(hours, return_counts = True) |
| 156 | + hours_u, cnt = np.unique(hours, return_counts=True) |
142 | 157 | precip_qpe60 = {} |
143 | 158 |
|
144 | 159 | for m in models: |
145 | 160 | data = precip_qpe[m] |
146 | | - precip_qpe60[m] = np.array([np.nanmean(data[hours == h], axis = 0) |
147 | | - for h in hours_u[cnt == 6]] ) |
148 | | - if 'CPC.CV' in list_models: |
149 | | - logging.info('Retrieving CPC.CV data. at hourly resolution...') |
150 | | - |
151 | | - precip_qpe60['CPC.CV'] = [] |
| 161 | + precip_qpe60[m] = np.array( |
| 162 | + [np.nanmean(data[hours == h], axis=0) for h in hours_u[cnt == 6]] |
| 163 | + ) |
| 164 | + if "CPC.CV" in list_models: |
| 165 | + logging.info("Retrieving CPC.CV data. at hourly resolution...") |
| 166 | + |
| 167 | + precip_qpe60["CPC.CV"] = [] |
152 | 168 | for h in hours_u: |
153 | | - precip_qpe60['CPC.CV'] .append(retrieve_CPCCV(h, stations)) |
154 | | - precip_qpe60['CPC.CV'] = np.array(precip_qpe60['CPC.CV']) |
155 | | - models.append('CPC.CV') |
156 | | - |
157 | | - precip_ref60 = np.array([np.nanmean(precip_ref[hours == h], axis = 0) |
158 | | - for h in hours_u[cnt == 6]] ) |
159 | | - |
| 169 | + precip_qpe60["CPC.CV"].append(retrieve_CPCCV(h, stations)) |
| 170 | + precip_qpe60["CPC.CV"] = np.array(precip_qpe60["CPC.CV"]) |
| 171 | + models.append("CPC.CV") |
| 172 | + |
| 173 | + precip_ref60 = np.array( |
| 174 | + [np.nanmean(precip_ref[hours == h], axis=0) for h in hours_u[cnt == 6]] |
| 175 | + ) |
| 176 | + |
160 | 177 | scores60 = {} |
161 | 178 | # COmpute 60min scores |
162 | 179 | valid_ref = np.isfinite(precip_ref60.ravel()) |
163 | 180 | for m in models: |
164 | | - scores60[m] = perfscores(precip_qpe60[m].ravel()[valid_ref], |
165 | | - precip_ref60.ravel()[valid_ref], |
166 | | - bounds60) |
167 | | - |
| 181 | + scores60[m] = perfscores( |
| 182 | + precip_qpe60[m].ravel()[valid_ref], |
| 183 | + precip_ref60.ravel()[valid_ref], |
| 184 | + bounds60, |
| 185 | + ) |
| 186 | + |
168 | 187 | # Make score plots |
169 | | - timerange = datetime.datetime.strftime(tsteps[0], '%Y%m%d%H%M') + '_' +\ |
170 | | - datetime.datetime.strftime(tsteps[-1], '%Y%m%d%H%M') |
171 | | - |
| 188 | + timerange = ( |
| 189 | + datetime.datetime.strftime(tsteps[0], "%Y%m%d%H%M") |
| 190 | + + "_" |
| 191 | + + datetime.datetime.strftime(tsteps[-1], "%Y%m%d%H%M") |
| 192 | + ) |
| 193 | + |
172 | 194 | # Save the data as parquet |
173 | 195 | for m in models: |
174 | | - df_precip = pd.DataFrame(precip_qpe[m],columns=stations, index=tsteps) |
175 | | - df_precip.to_csv(outputfolder+'/'+str(m)+'_qpe10min_'+timerange+'.csv', float_format='%.3f') |
176 | | - |
177 | | - title = datetime.datetime.strftime(tsteps[0], '%d %b %Y %H:%M') |
178 | | - title += ' - ' + datetime.datetime.strftime(tsteps[-1], '%d %b %Y %H:%M') |
179 | | - score_plot(scores10, title + ', agg = 10 min ', figsize = (13,8)) |
180 | | - plt.savefig(outputfolder + 'scores_agg10_' + timerange+ '.png', |
181 | | - bbox_inches='tight', |
182 | | - dpi = 300) |
183 | | - |
184 | | - score_plot(scores60, title + ', agg = 60 min ', figsize = (13,8)) |
185 | | - plt.savefig(outputfolder + 'scores_agg60_' + timerange+ '.png', |
186 | | - bbox_inches='tight', |
187 | | - dpi = 300) |
188 | | - |
| 196 | + df_precip = pd.DataFrame(precip_qpe[m], columns=stations, index=tsteps) |
| 197 | + df_precip.to_csv( |
| 198 | + outputfolder + "/" + str(m) + "_qpe10min_" + timerange + ".csv", |
| 199 | + float_format="%.3f", |
| 200 | + ) |
| 201 | + |
| 202 | + title = datetime.datetime.strftime(tsteps[0], "%d %b %Y %H:%M") |
| 203 | + title += " - " + datetime.datetime.strftime(tsteps[-1], "%d %b %Y %H:%M") |
| 204 | + score_plot(scores10, title + ", agg = 10 min ", figsize=(13, 8)) |
| 205 | + plt.savefig( |
| 206 | + outputfolder + "scores_agg10_" + timerange + ".png", |
| 207 | + bbox_inches="tight", |
| 208 | + dpi=300, |
| 209 | + ) |
| 210 | + |
| 211 | + score_plot(scores60, title + ", agg = 60 min ", figsize=(13, 8)) |
| 212 | + plt.savefig( |
| 213 | + outputfolder + "scores_agg60_" + timerange + ".png", |
| 214 | + bbox_inches="tight", |
| 215 | + dpi=300, |
| 216 | + ) |
| 217 | + |
189 | 218 | # Make scatterplots |
190 | | - qpe_scatterplot(precip_qpe, precip_ref, figsize = (10,8.5), |
191 | | - title_prefix = title + ', agg = 10 min ') |
192 | | - plt.savefig(outputfolder + 'scatterplots10_' + timerange+ '.png', |
193 | | - bbox_inches='tight', |
194 | | - dpi = 300) |
195 | | - qpe_scatterplot(precip_qpe60, precip_ref60, figsize = (10,8.5), |
196 | | - title_prefix = title + ', agg = 60 min ') |
197 | | - plt.savefig(outputfolder + 'scatterplots60_' + timerange+ '.png', |
198 | | - bbox_inches='tight', |
199 | | - dpi = 300) |
200 | | - |
201 | | - |
| 219 | + qpe_scatterplot( |
| 220 | + precip_qpe, |
| 221 | + precip_ref, |
| 222 | + figsize=(10, 8.5), |
| 223 | + title_prefix=title + ", agg = 10 min ", |
| 224 | + ) |
| 225 | + plt.savefig( |
| 226 | + outputfolder + "scatterplots10_" + timerange + ".png", |
| 227 | + bbox_inches="tight", |
| 228 | + dpi=300, |
| 229 | + ) |
| 230 | + qpe_scatterplot( |
| 231 | + precip_qpe60, |
| 232 | + precip_ref60, |
| 233 | + figsize=(10, 8.5), |
| 234 | + title_prefix=title + ", agg = 60 min ", |
| 235 | + ) |
| 236 | + plt.savefig( |
| 237 | + outputfolder + "scatterplots60_" + timerange + ".png", |
| 238 | + bbox_inches="tight", |
| 239 | + dpi=300, |
| 240 | + ) |
0 commit comments