/
main.ncl
486 lines (417 loc) · 18.2 KB
/
main.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
; #############################################################################
; MAIN SCRIPT FOR PERFORMANCE METRICS
; Authors: Mattia Righi (DLR, Germany) and Franziska Winterstein (DLR, Germany)
; ESMVal project
; #############################################################################
;
; Description
; Calculates and (optionally) plots annual/seasonal cycles, zonal means,
; lat-lon fields and time-lat-lon fields from input 2D/3D monthly data.
; The calculated fields can be plotted as difference w.r.t. a given
; reference dataset. It also calculates grading and taylor metrics.
; Input data have to be regridded to a common grid in the preprocessor.
;
; Required diag_script_info attributes
; plot_type: cycle (time), zonal (plev, lat), latlon (lat, lon) or
; cycle_latlon (time, lat, lon)
; time_avg: type of time average (see time_operations in
; diag_scripts/shared/statistics.ncl)
; region: selected region (see select_region in
; diag_scripts/shared/latlon.ncl)
;
; Optional diag_script_info attributes
; styleset (for cycle): as in diag_scripts/shared/plot/style.ncl functions
; plot_stddev (for cycle): plot standard deviation
; legend_outside (for cycle): save legend in a separate file
; t_test (for zonal and latlon): calculate t-test in difference plots
; (default: False)
; conf_level (for zonal and latlon): confidence level for the t-test
; (default: False)
; projection: map projection for lat-lon plots (default:
; CylindricalEquidistant)
; plot_diff: draw difference plots (default: False)
; calc_grading: calculate grading (default: False)
; stippling: use stippling to mark stat. significant differences (default:
; False = mask out non-significant differences in gray)
; show_global_avg: diplay global avaerage as right string on lat-lon plots
; (default: False)
; metric: grading metric (if calc_grading is True)
; normalization: metric normalization (for RMSD and BIAS metrics)
; abs_levs: (min, max, spacing) contour levels for absolute plot
; diff_levs: (min, max, spacing) contour levels for difference plot
; zonal_cmap (for zonal): color table (default: "amwg_blueyellowred")
; zonal_ymin (for zonal): minimum pressure on the plots (default: 5. hPa)
; latlon_cmap (for latlon): color table (default: "amwg_blueyellowred")
; plot_units: plotting units (if different from standard CMOR units)
; add_tropopause: add an optional tropopause outline to the zonal plots
; res_*: any resource as applied to ncl gsn_csm_press_hgt plots
;
; Required variable_info attributes:
; reference_dataset: reference dataset to compare with (usualy observations)
;
; Optional variable_info attributes:
; alternative_dataset: a second dataset to compare with
;
; Caveats
;
; Modification history
; 20221209-winterstein_franziska: added optional tropopause and plotting
; 20220609-bock_lisa: added calculation of multi-model mean and median for
; each project
; 20211014-bock_lisa: added sorting by project
; 20200506-gier_bettina: implemented support for multiple occurence of
; models with different experiments and ensembles
; 20190405-righi_mattia: added provenance logging
; 20190315-hassler_birgit: extended to smpi metric
; 20180503-righi_mattia: completely rewritten and modularized
; 20171215-righi_mattia: merged with perfmetrics_grading and
; permetrics_taylor.ncl
; 20171124-righi_mattia: completely revised to adapt it to the new backend
; (level selection, regridding and masking now done
; by the python preprocessor)
; 20161220-lauer_axel: added option to set map projection for lat-lon plots
; (diag_script_info@projection)
; added option to choose how to plot t-test results:
; stippling or masking out in gray (lat-lon plots only)
; 20161019-lauer_axel: changed plotting of t-test results:
; now stippling significant grid cells (old version:
; masking out non-significant values in gray)
; 20160628-righi_mattia: moving ref_model specification from cfg- files to
; recipe file
; 20160628-senftleben_daniel: added regridding for irregular grids
; (ESMF_regrid).
; 20151027-lauer_axel: moved call to 'write_references' to the beginning
; of the code.
; 20151013-righi_mattia: fixed t-test mask in lat-lon difference plots.
; 20150325-lauer_axel: modified reference tags used for acknowledgements
; (projects, observations, etc.).
; 20150119-gottschaldt_klaus-dirk: removed "grid", "region" from req_atts
; (for T2Ms vmrco).
; 20150113-gottschaldt_klaus-dirk: reconciled generalised regridding with
; T1* & T0*
; 20140905-righi_mattia: consistent regridding and missing values mask.
; 20140701-gottschaldt_klaus-dirk: Adapted for T1M.
; 20140630-gottschaldt_klaus-dirk: Adapted for T0Ms.
; 20131203-winterstein_franziska: written.
;
; #############################################################################
load "$diag_scripts/../interface_scripts/interface.ncl"
load "$diag_scripts/shared/latlon.ncl"
load "$diag_scripts/shared/statistics.ncl"
load "$diag_scripts/shared/regridding.ncl"
load "$diag_scripts/shared/ensemble.ncl"
load "$diag_scripts/shared/scaling.ncl"
load "$diag_scripts/shared/plot/style.ncl"
begin
enter_msg(DIAG_SCRIPT, "")
vars = metadata_att_as_array(variable_info, "short_name")
; Get variables and datasets
var0 = variable_info[0]@short_name
info_items = select_metadata_by_name(input_file_info, var0)
nDatasets = ListCount(info_items)
if (dimsizes(vars) .gt. 1) then
var1 = variable_info[1]@short_name
info_items_1 = select_metadata_by_name(input_file_info, var1)
cnt = ListCount(info_items_1)
do ii = 0, cnt - 1
ListAppend(info_items, info_items_1[ii])
end do
nDatasets = ListCount(info_items)
end if
; Check required diag_script_info attributes
exit_if_missing_atts(diag_script_info, (/"plot_type", "time_avg", "region"/))
; Define region
region = select_region(diag_script_info@region)
; Store required attributes
ptype = diag_script_info@plot_type
if (all(ptype.ne. \
(/"cycle", "zonal", "latlon", "cycle_latlon", "cycle_zonal"/))) then
error_msg("f", DIAG_SCRIPT, "", "plot_type " + ptype + " is not a " + \
"supported plot_type in this diagnostic")
end if
; Check for plot-type specific settings
if (ptype.eq."cycle") then
exit_if_missing_atts(diag_script_info, \
(/"legend_outside", "styleset", "plot_stddev"/))
end if
if ((ptype.eq."zonal" .or. ptype.eq."latlon") .and. \
diag_script_info@region.ne."global") then
error_msg("f", DIAG_SCRIPT, "", "plot_type " + ptype + \
" implemented only for region 'global'")
end if
; Set default values for non-required diag_script_info attributes
set_default_att(diag_script_info, "projection", "CylindricalEquidistant")
set_default_att(diag_script_info, "plot_diff", False)
set_default_att(diag_script_info, "calc_grading", False)
set_default_att(diag_script_info, "stippling", False)
set_default_att(diag_script_info, "t_test", False)
set_default_att(diag_script_info, "show_global_avg", False)
set_default_att(diag_script_info, "zonal_ymin", 5.)
set_default_att(diag_script_info, "zonal_cmap", "amwg_blueyellowred")
set_default_att(diag_script_info, "latlon_cmap", "amwg_blueyellowred")
set_default_att(diag_script_info, "add_tropopause", False)
; Check consistency of diff plots settings
if (diag_script_info@t_test .and. .not.diag_script_info@plot_diff) then
error_msg("f", DIAG_SCRIPT, "", "plot_diff must be True to apply t-test")
end if
if (diag_script_info@t_test .and. .not.diag_script_info@conf_level) then
error_msg("f", DIAG_SCRIPT, "", \
"conf_level must be specified to apply t-test")
end if
; Check metric
if (diag_script_info@calc_grading) then
exit_if_missing_atts(diag_script_info, (/"metric", "normalization"/))
if (dimsizes(diag_script_info@metric).ne.\
dimsizes(diag_script_info@normalization)) then
error_msg("f", DIAG_SCRIPT, "", "normalization must be " + \
"provided for each requested metric")
end if
end if
; Set dataset names
datasetnames = metadata_att_as_array(info_items, "dataset")
; Save projects
projectnames = metadata_att_as_array(info_items, "project")
; Extend model name by ensemble/experiment on multiple occurence
if (.not. (dimsizes(datasetnames) .eq. \
count_unique_values(datasetnames))) then
new_datasetnames = datasetnames
experiments = metadata_att_as_array(info_items, "exp")
ensembles = metadata_att_as_array(info_items, "ensemble")
do imod = 0, dimsizes(datasetnames) - 1
do jmod = 0, dimsizes(datasetnames) - 1
if imod.eq.jmod then
continue
else
if (datasetnames(imod) .eq. datasetnames(jmod)) then
if (experiments(imod) .ne. experiments(jmod)) then
new_datasetnames(imod) = new_datasetnames(imod) + " " + \
experiments(imod)
break
end if
end if
end if
end do
end do
do imod = 0, dimsizes(datasetnames) - 1
do jmod = 0, dimsizes(datasetnames) - 1
if imod.eq.jmod then
continue
else
if (datasetnames(imod) .eq. datasetnames(jmod)) then
if (ensembles(imod) .ne. ensembles(jmod)) then
new_datasetnames(imod) = new_datasetnames(imod) + " " + \
ensembles(imod)
break
end if
end if
end if
end do
end do
datasetnames = new_datasetnames
delete(new_datasetnames)
end if
; Save list of preproc files for provenance in collect.ncl
preproc_files = metadata_att_as_array(info_items, "filename")
; Check for reference dataset definition
if (variable_info[0]@reference_dataset.eq."None") then
error_msg("f", DIAG_SCRIPT, "", "no reference dataset is specified")
end if
; Set index of the reference (and alternative) dataset
l_altern = False
nobs = 1
ref_ind = ind(datasetnames.eq.variable_info[0]@reference_dataset)
ref_inds = ref_ind
if (isatt(variable_info[0], "alternative_dataset")) then
l_altern = True
nobs = 2
alt_ind = ind(datasetnames.eq.variable_info[0]@alternative_dataset)
ref_inds := (/ref_ind, alt_ind/)
end if
; Create output plot directory
plot_dir = config_user_info@plot_dir
system("mkdir -p " + plot_dir)
; Plot file type
file_type = config_user_info@output_file_type
if (ismissing(file_type)) then
file_type = "ps"
end if
; Grading settings
if (diag_script_info@calc_grading) then
; Define grading arrays
nmetrics = dimsizes(diag_script_info@metric)
ncdf_dir = new(nmetrics, string)
nModels = dimsizes(datasetnames) - nobs
grading = new((/nmetrics, 1, nModels, nobs/), float)
grading!0 = "metric"
grading!1 = "diagnostics" ; dummy coord. to facilitate appending
grading!2 = "models"
grading!3 = "reference"
grading&diagnostics = \
variable_info[0]@diagnostic + "-" + diag_script_info@region
grading&models = remove_index(datasetnames, ref_inds)
if (isdim(grading, "reference")) then
grading&reference = datasetnames(ref_inds)
end if
grading@projects = str_join(remove_index(projectnames, ref_inds), " ")
; Special case Taylor
if (any(diag_script_info@metric.eq."taylor")) then
nModels = dimsizes(datasetnames) - 1 ; always 1 reference dataset
taylor = new((/1, nModels, 2/), float)
taylor!0 = "diagnostics" ; dummy coord. to facilitate appending
taylor!1 = "models"
taylor!2 = "statistic"
taylor&diagnostics = \
variable_info[0]@diagnostic + "-" + diag_script_info@region
taylor&statistic = (/"stddev_ratio", "correlation"/)
taylor&models = remove_index(datasetnames, ref_ind)
end if
; Special case SMPI
if (any(diag_script_info@metric.eq."SMPI")) then
nModels = dimsizes(datasetnames) - 1 ; always 1 reference model
smpi = new((/diag_script_info@smpi_n_bootstrap + 1, nModels/), float)
smpi!0 = "bootstrap_member"
smpi!1 = "models"
smpi&bootstrap_member = ispan(0, diag_script_info@smpi_n_bootstrap, 1)
smpi&models = remove_index(datasetnames, ref_ind)
end if
; Define grading filename
do met = 0, nmetrics - 1
ncdf_dir(met) = config_user_info@work_dir + \
diag_script_info@metric(met) + "_" + \
variable_info[0]@diagnostic + "-" + diag_script_info@region + ".nc"
end do
end if
; Load plot-type-specific script
loadscript("$diag_scripts/perfmetrics/" + ptype + ".ncl")
end
begin
; Call plot-type-specific script
perfmetrics_ptype_script()
; Finalize grading calculations
if (diag_script_info@calc_grading) then
do met = 0, nmetrics - 1
if (diag_script_info@metric(met).eq."RMSD") then
metric = grading(met, :, :, :)
metric@title = diag_script_info@metric(met) + " metric"
metric@long_name = \
"Grading table of metric " + diag_script_info@metric(met)
metric@var = "grade"
; Normalization
do iobs = 0, nobs - 1
metric(:, :, iobs) = \
normalize_metric(metric(:, :, iobs), \
diag_script_info@normalization(met))
end do
; Reduce dimensionality if no alternative dataset
if (.not.l_altern) then
metric := metric(:, :, 0)
delete(metric@reference)
end if
; Provenance information
statistics = (/"rmsd"/)
authors = (/"winterstein_franziska", "righi_mattia", \
"eyring_veronika"/)
references = (/"righi15gmd", "gleckler08jgr"/)
elseif (diag_script_info@metric(met).eq."BIAS") then
metric = grading(met, :, :, :)
metric@title = diag_script_info@metric(met) + " metric"
metric@long_name = \
"Grading table of metric " + diag_script_info@metric(met)
metric@var = "grade"
; Normalization
do iobs = 0, nobs - 1
metric(:, :, iobs) = \
normalize_metric(metric(:, :, iobs), \
diag_script_info@normalization(met))
end do
; Reduce dimensionality if no alternative dataset
if (.not.l_altern) then
metric := metric(:, :, 0)
delete(metric@reference)
end if
; Provenance information
statistics = (/"diff"/)
authors = (/"winterstein_franziska", "righi_mattia", \
"eyring_veronika"/)
references = (/"righi15gmd", "gleckler08jgr"/)
elseif (diag_script_info@metric(met).eq."taylor") then
metric = taylor
metric@title = diag_script_info@metric(met) + " metric"
metric@long_name = \
"Grading table of metric " + diag_script_info@metric(met)
metric@var = "grade"
; Provenance information
statistics = (/"rmsd", "corr"/)
authors = (/"winterstein_franziska", "righi_mattia", \
"eyring_veronika"/)
references = (/"righi15gmd", "gleckler08jgr"/)
elseif (diag_script_info@metric(met).eq."SMPI") then
metric = smpi
metric@title = "metrics"
metric@long_name = "1 variable's Performance Index for " + \
"the Single Model Performance Index"
metric@var = "performance_index"
metric@invar = var0
metric@ensemble_name = diag_script_info@normalization(met)
; Normalization
ens_idx = new(dimsizes(metric&models), integer)
atts = True
atts@project = diag_script_info@normalization(met)
info = select_metadata_by_atts(input_file_info, atts)
delete(atts)
do ii = 0, dimsizes(ens_idx) - 1
if (dimsizes(info).ne.0) then
ens_idx(ii) = ii
end if
end do
if (all(ismissing(ens_idx))) then
error_msg("f", DIAG_SCRIPT, "", "No datasets for the selected " + \
"normalization (" + diag_script_info@normalization(met) + \
") found")
end if
ens_idx := ens_idx(ind(.not.ismissing(ens_idx)))
do iboot = 0, dimsizes(metric&bootstrap_member)-1
metric(iboot, :) = metric(iboot, :) / avg(metric(iboot, ens_idx))
end do
; Provenance information
statistics = "smpi"
authors = (/"gier_bettina", "hassler_birgit"/)
references = (/"rk2008bams"/)
else
error_msg("f", DIAG_SCRIPT, "", "unrecognized metric " + \
diag_script_info@metric(met))
end if
; Common attributes
metric@metric = diag_script_info@metric(met)
metric@diag_script = (/DIAG_SCRIPT/)
metric@region = region@name
metric@num_preproc_files = dimsizes(preproc_files)
do imod = 0, metric@num_preproc_files - 1
num_preproc = "preproc_file_" + imod
metric@$num_preproc$ = preproc_files(imod)
end do
metric@ncdf_dir = ncdf_dir(met)
ncdf_outfile = ncdf_write(metric, metric@ncdf_dir)
; Call provenance logger
log_provenance(ncdf_outfile, "n/a", "n/a", statistics, \
diag_script_info@region, "other", authors, references, \
preproc_files)
delete([/statistics, authors, references/])
; Write results of temporary grading list
temp_dir = config_user_info@work_dir + "/" + \
diag_script_info@metric(met) + ".nc"
if (fileexists(temp_dir)) then
system("rm -f " + temp_dir)
end if
ncdf_char = tochar(ncdf_dir(met))
temp_list = new((/1, dimsizes(ncdf_char)/), character)
temp_list(0, :) = ncdf_char
delete(ncdf_char)
; Create new file and add list
temp = addfile(temp_dir, "c")
temp->temp_list = temp_list
delete([/metric, temp_dir, temp_list/])
end do
end if
leave_msg(DIAG_SCRIPT, "")
end