Skip to content

Commit 6cdf2a2

Browse files
lee1043lee1043
authored
Merge pull request #931 from PCMDI/mean_clim_loca2_analysis
PMP mean climate regional application
2 parents 26326e5 + 847a1a2 commit 6cdf2a2

File tree

6 files changed

+899
-530
lines changed

6 files changed

+899
-530
lines changed

‎doc/jupyter/Demo/Demo_1b_mean_climate.ipynb

Lines changed: 739 additions & 439 deletions
Large diffs are not rendered by default.

‎doc/obs_info_dictionary.json

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -100,12 +100,12 @@
100100
},
101101
"GPCP-2-3": {
102102
"CMIP_CMOR_TABLE": "Amon",
103-
"MD5sum": "3792901034585d3d495722f10a0dfecb",
103+
"MD5sum": "036e14d0124f4d00921d58e3c1c3e9e1",
104104
"RefTrackingDate": "Wed Aug 4 12:22:10 2021",
105-
"filename": "pr_mon_GPCP-2-3_PCMDI_gn.200301-201812.AC.v20210804.nc",
105+
"filename": "pr_monC_GPCP-2-3_PCMDI_gn_197901-201907.AC.v20230516.nc",
106106
"period": "200301-201812",
107107
"shape": "(12, 72, 144)",
108-
"template": "pr/GPCP-2-3/v20210804/pr_mon_GPCP-2-3_PCMDI_gn.200301-201812.AC.v20210804.nc"
108+
"template": "NOAA-NCEI/GPCP-2-3/monC/pr/gn/v20230516/pr_monC_GPCP-2-3_PCMDI_gn_197901-201907.AC.v20230516.nc"
109109
},
110110
"TRMM-3B43v-7": {
111111
"CMIP_CMOR_TABLE": "Amon",

‎pcmdi_metrics/mean_climate/lib/compute_metrics.py

Lines changed: 54 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
import pcmdi_metrics
44

55

6-
def compute_metrics(Var, dm, do, debug=False):
6+
def compute_metrics(Var, dm, do, debug=False, time_dim_sync=False):
77
# Var is sometimes sent with level associated
88
var = Var.split("_")[0]
99
# Did we send data? Or do we just want the info?
@@ -35,44 +35,43 @@ def compute_metrics(Var, dm, do, debug=False):
3535

3636
# unify time and time bounds between observation and model
3737
if debug:
38-
print("before time and time bounds unifying")
3938
print("dm.time: ", dm["time"])
4039
print("do.time: ", do["time"])
4140

42-
# Below is temporary...
43-
dm["time"] = do["time"]
44-
dm[dm.time.attrs["bounds"]] = do[do.time.attrs["bounds"]]
41+
if time_dim_sync:
42+
# Below is temporary...
43+
dm["time"] = do["time"]
44+
dm[dm.time.encoding["bounds"]] = do[do.time.attrs["bounds"]]
4545

46-
if debug:
47-
print("after time and time bounds unifying")
48-
print("dm.time: ", dm["time"])
49-
print("do.time: ", do["time"])
46+
if debug:
47+
print("time and time bounds synced")
48+
print("dm.time: ", dm["time"])
49+
print("do.time: ", do["time"])
5050

51-
# if debug:
52-
# dm.to_netcdf('dm.nc')
53-
# do.to_netcdf('do.nc')
51+
dm.to_netcdf("dm.nc")
52+
do.to_netcdf("do.nc")
5453

5554
metrics_dictionary = OrderedDict()
5655

57-
# SET CONDITIONAL ON INPUT VARIABLE
58-
if var == "pr":
59-
conv = 86400.0
60-
else:
61-
conv = 1.0
62-
6356
if var in ["hus"]:
6457
sig_digits = ".5f"
6558
else:
6659
sig_digits = ".3f"
6760

6861
# CALCULATE ANNUAL CYCLE SPACE-TIME RMS, CORRELATIONS and STD
6962
print("compute_metrics-CALCULATE ANNUAL CYCLE SPACE-TIME RMS, CORRELATIONS and STD")
63+
7064
print("compute_metrics, rms_xyt")
7165
rms_xyt = pcmdi_metrics.mean_climate.lib.rms_xyt(dm, do, var)
66+
print("compute_metrics, rms_xyt:", rms_xyt)
67+
7268
print("compute_metrics, stdObs_xyt")
7369
stdObs_xyt = pcmdi_metrics.mean_climate.lib.std_xyt(do, var)
70+
print("compute_metrics, stdObs_xyt:", stdObs_xyt)
71+
7472
print("compute_metrics, std_xyt")
7573
std_xyt = pcmdi_metrics.mean_climate.lib.std_xyt(dm, var)
74+
print("compute_metrics, std_xyt:", std_xyt)
7675

7776
# CALCULATE ANNUAL MEANS
7877
print("compute_metrics-CALCULATE ANNUAL MEANS")
@@ -81,19 +80,23 @@ def compute_metrics(Var, dm, do, debug=False):
8180
# CALCULATE ANNUAL MEAN BIAS
8281
print("compute_metrics-CALCULATE ANNUAL MEAN BIAS")
8382
bias_xy = pcmdi_metrics.mean_climate.lib.bias_xy(dm_am, do_am, var)
83+
print("compute_metrics-CALCULATE ANNUAL MEAN BIAS, bias_xy:", bias_xy)
8484

8585
# CALCULATE MEAN ABSOLUTE ERROR
8686
print("compute_metrics-CALCULATE MSE")
8787
mae_xy = pcmdi_metrics.mean_climate.lib.meanabs_xy(dm_am, do_am, var)
88+
print("compute_metrics-CALCULATE MSE, mae_xy:", mae_xy)
8889

8990
# CALCULATE ANNUAL MEAN RMS (centered and uncentered)
9091
print("compute_metrics-CALCULATE MEAN RMS")
9192
rms_xy = pcmdi_metrics.mean_climate.lib.rms_xy(dm_am, do_am, var)
9293
rmsc_xy = pcmdi_metrics.mean_climate.lib.rmsc_xy(dm_am, do_am, var)
94+
print("compute_metrics-CALCULATE MEAN RMS: rms_xy, rmsc_xy: ", rms_xy, rmsc_xy)
9395

9496
# CALCULATE ANNUAL MEAN CORRELATION
9597
print("compute_metrics-CALCULATE MEAN CORR")
9698
cor_xy = pcmdi_metrics.mean_climate.lib.cor_xy(dm_am, do_am, var)
99+
print("compute_metrics-CALCULATE MEAN CORR: cor_xy:", cor_xy)
97100

98101
# CALCULATE ANNUAL OBS and MOD STD
99102
print("compute_metrics-CALCULATE ANNUAL OBS AND MOD STD")
@@ -155,25 +158,23 @@ def compute_metrics(Var, dm, do, debug=False):
155158
):
156159
metrics_dictionary[stat] = OrderedDict()
157160

158-
metrics_dictionary["mean-obs_xy"]["ann"] = format(meanObs_xy * conv, sig_digits)
159-
metrics_dictionary["mean_xy"]["ann"] = format(mean_xy * conv, sig_digits)
160-
metrics_dictionary["std-obs_xy"]["ann"] = format(stdObs_xy * conv, sig_digits)
161-
metrics_dictionary["std_xy"]["ann"] = format(std_xy * conv, sig_digits)
162-
metrics_dictionary["std-obs_xyt"]["ann"] = format(stdObs_xyt * conv, sig_digits)
163-
metrics_dictionary["std_xyt"]["ann"] = format(std_xyt * conv, sig_digits)
164-
metrics_dictionary["std-obs_xy_devzm"]["ann"] = format(
165-
stdObs_xy_devzm * conv, sig_digits
166-
)
167-
metrics_dictionary["std_xy_devzm"]["ann"] = format(std_xy_devzm * conv, sig_digits)
168-
metrics_dictionary["rms_xyt"]["ann"] = format(rms_xyt * conv, sig_digits)
169-
metrics_dictionary["rms_xy"]["ann"] = format(rms_xy * conv, sig_digits)
170-
metrics_dictionary["rmsc_xy"]["ann"] = format(rmsc_xy * conv, sig_digits)
161+
metrics_dictionary["mean-obs_xy"]["ann"] = format(meanObs_xy, sig_digits)
162+
metrics_dictionary["mean_xy"]["ann"] = format(mean_xy, sig_digits)
163+
metrics_dictionary["std-obs_xy"]["ann"] = format(stdObs_xy, sig_digits)
164+
metrics_dictionary["std_xy"]["ann"] = format(std_xy, sig_digits)
165+
metrics_dictionary["std-obs_xyt"]["ann"] = format(stdObs_xyt, sig_digits)
166+
metrics_dictionary["std_xyt"]["ann"] = format(std_xyt, sig_digits)
167+
metrics_dictionary["std-obs_xy_devzm"]["ann"] = format(stdObs_xy_devzm, sig_digits)
168+
metrics_dictionary["std_xy_devzm"]["ann"] = format(std_xy_devzm, sig_digits)
169+
metrics_dictionary["rms_xyt"]["ann"] = format(rms_xyt, sig_digits)
170+
metrics_dictionary["rms_xy"]["ann"] = format(rms_xy, sig_digits)
171+
metrics_dictionary["rmsc_xy"]["ann"] = format(rmsc_xy, sig_digits)
171172
metrics_dictionary["cor_xy"]["ann"] = format(cor_xy, sig_digits)
172-
metrics_dictionary["bias_xy"]["ann"] = format(bias_xy * conv, sig_digits)
173-
metrics_dictionary["mae_xy"]["ann"] = format(mae_xy * conv, sig_digits)
173+
metrics_dictionary["bias_xy"]["ann"] = format(bias_xy, sig_digits)
174+
metrics_dictionary["mae_xy"]["ann"] = format(mae_xy, sig_digits)
174175
# ZONAL MEAN CONTRIBUTIONS
175-
metrics_dictionary["rms_y"]["ann"] = format(rms_y * conv, sig_digits)
176-
metrics_dictionary["rms_devzm"]["ann"] = format(rms_xy_devzm * conv, sig_digits)
176+
metrics_dictionary["rms_y"]["ann"] = format(rms_y, sig_digits)
177+
metrics_dictionary["rms_devzm"]["ann"] = format(rms_xy_devzm, sig_digits)
177178

178179
# CALCULATE SEASONAL MEANS
179180
for sea in ["djf", "mam", "jja", "son"]:
@@ -195,17 +196,15 @@ def compute_metrics(Var, dm, do, debug=False):
195196
meanObs_xy_sea = pcmdi_metrics.mean_climate.lib.mean_xy(do_sea, var)
196197
mean_xy_sea = pcmdi_metrics.mean_climate.lib.mean_xy(dm_sea, var)
197198

198-
metrics_dictionary["bias_xy"][sea] = format(bias_sea * conv, sig_digits)
199-
metrics_dictionary["rms_xy"][sea] = format(rms_sea * conv, sig_digits)
200-
metrics_dictionary["rmsc_xy"][sea] = format(rmsc_sea * conv, sig_digits)
199+
metrics_dictionary["bias_xy"][sea] = format(bias_sea, sig_digits)
200+
metrics_dictionary["rms_xy"][sea] = format(rms_sea, sig_digits)
201+
metrics_dictionary["rmsc_xy"][sea] = format(rmsc_sea, sig_digits)
201202
metrics_dictionary["cor_xy"][sea] = format(cor_sea, ".2f")
202-
metrics_dictionary["mae_xy"][sea] = format(mae_sea * conv, sig_digits)
203-
metrics_dictionary["std-obs_xy"][sea] = format(stdObs_xy_sea * conv, sig_digits)
204-
metrics_dictionary["std_xy"][sea] = format(std_xy_sea * conv, sig_digits)
205-
metrics_dictionary["mean-obs_xy"][sea] = format(
206-
meanObs_xy_sea * conv, sig_digits
207-
)
208-
metrics_dictionary["mean_xy"][sea] = format(mean_xy_sea * conv, sig_digits)
203+
metrics_dictionary["mae_xy"][sea] = format(mae_sea, sig_digits)
204+
metrics_dictionary["std-obs_xy"][sea] = format(stdObs_xy_sea, sig_digits)
205+
metrics_dictionary["std_xy"][sea] = format(std_xy_sea, sig_digits)
206+
metrics_dictionary["mean-obs_xy"][sea] = format(meanObs_xy_sea, sig_digits)
207+
metrics_dictionary["mean_xy"][sea] = format(mean_xy_sea, sig_digits)
209208

210209
rms_mo_l = []
211210
rmsc_mo_l = []
@@ -251,15 +250,15 @@ def compute_metrics(Var, dm, do, debug=False):
251250
meanObs_xy_mo = pcmdi_metrics.mean_climate.lib.mean_xy(do_mo, var)
252251
mean_xy_mo = pcmdi_metrics.mean_climate.lib.mean_xy(dm_mo, var)
253252

254-
rms_mo_l.append(format(rms_mo * conv, sig_digits))
255-
rmsc_mo_l.append(format(rmsc_mo * conv, sig_digits))
253+
rms_mo_l.append(format(rms_mo, sig_digits))
254+
rmsc_mo_l.append(format(rmsc_mo, sig_digits))
256255
cor_mo_l.append(format(cor_mo, ".2f"))
257-
mae_mo_l.append(format(mae_mo * conv, sig_digits))
258-
bias_mo_l.append(format(bias_mo * conv, sig_digits))
259-
stdObs_xy_mo_l.append(format(stdObs_xy_mo * conv, sig_digits))
260-
std_xy_mo_l.append(format(std_xy_mo * conv, sig_digits))
261-
meanObs_xy_mo_l.append(format(meanObs_xy_mo * conv, sig_digits))
262-
mean_xy_mo_l.append(format(mean_xy_mo * conv, sig_digits))
256+
mae_mo_l.append(format(mae_mo, sig_digits))
257+
bias_mo_l.append(format(bias_mo, sig_digits))
258+
stdObs_xy_mo_l.append(format(stdObs_xy_mo, sig_digits))
259+
std_xy_mo_l.append(format(std_xy_mo, sig_digits))
260+
meanObs_xy_mo_l.append(format(meanObs_xy_mo, sig_digits))
261+
mean_xy_mo_l.append(format(mean_xy_mo, sig_digits))
263262

264263
metrics_dictionary["bias_xy"]["CalendarMonths"] = bias_mo_l
265264
metrics_dictionary["rms_xy"]["CalendarMonths"] = rms_mo_l
@@ -276,8 +275,8 @@ def compute_metrics(Var, dm, do, debug=False):
276275

277276
# ZONAL AND SEASONAL MEAN CONTRIBUTIONS
278277
# metrics_dictionary['rms_y'][sea] = format(
279-
# rms_y * conv,
278+
# rms_y,
280279
# sig_digits)
281280
# metrics_dictionary['rms_devzm'][sea] = format(
282-
# rms_xy_devzm * conv,
281+
# rms_xy_devzm,
283282
# sig_digits)

‎pcmdi_metrics/mean_climate/lib/load_and_regrid.py

Lines changed: 34 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,14 @@ def load_and_regrid(
3737
data_path, data_var=varname_in_file, decode_times=decode_times
3838
) # NOTE: decode_times=False will be removed once obs4MIP written using xcdat
3939

40+
# SET CONDITIONAL ON INPUT VARIABLE
41+
if varname == "pr":
42+
print("Adjust units for pr")
43+
if ds[varname_in_file].units == "kg m-2 s-1":
44+
ds[varname_in_file] = ds[varname_in_file] * 86400
45+
print("pr units adjusted to [mm d-1] from [kg m-2 s-1] by 86400 multiplied")
46+
47+
"""
4048
# calendar quality check
4149
if "calendar" in list(ds.time.attrs.keys()):
4250
if debug:
@@ -45,15 +53,18 @@ def load_and_regrid(
4553
if debug:
4654
print("ds.calendar:", ds.calendar)
4755
if ds.calendar != ds.time.attrs["calendar"]:
48-
print(
49-
'[WARNING]: calendar info mismatch. ds.time.attrs["calendar"] is adjusted to ds.calendar'
50-
)
51-
ds.time.attrs["calendar"] = ds.calendar
56+
ds.time.encoding["calendar"] = ds.calendar
57+
print('[WARNING]: calendar info mismatch. ds.time.attrs["calendar"] is adjusted to ds.calendar, ', ds.calendar)
5258
else:
5359
if "calendar" in ds.attrs.keys():
5460
ds.time.attrs["calendar"] = ds.calendar
61+
print('[WARNING]: calendar info not found for time axis. ds.time.attrs["calendar"] is adjusted to ds.calendar, ', ds.calendar)
62+
else:
63+
ds.time.attrs["calendar"] = 'standard'
64+
print('[WARNING]: calendar info not found for time axis. ds.time.attrs["calendar"] is adjusted to standard')
65+
"""
5566

56-
# time bound check -- add proper time bound info if cdms-generated annual cycle is loaded
67+
# time bound check #1 -- add proper time bound info if cdms-generated annual cycle is loaded
5768
if isinstance(
5869
ds.time.values[0], np.float64
5970
): # and "units" not in list(ds.time.attrs.keys()):
@@ -62,6 +73,13 @@ def load_and_regrid(
6273
if debug:
6374
print("decode_time done")
6475

76+
# time bound check #2 -- add time bounds itself it it is missing
77+
if "bounds" in list(ds.time.attrs.keys()):
78+
time_bnds_key = ds.time.attrs["bounds"]
79+
if time_bnds_key not in list(ds.keys()):
80+
ds = ds.bounds.add_missing_bounds(["T"])
81+
print("[WARNING]: bounds.add_missing_bounds conducted for T axis")
82+
6583
# level - extract a specific level if needed
6684
if level is not None:
6785
if isinstance(level, int) or isinstance(level, float):
@@ -114,8 +132,9 @@ def load_and_regrid(
114132

115133
if varname != varname_in_file:
116134
ds_regridded[varname] = ds_regridded[varname_in_file]
135+
ds_regridded = ds_regridded.drop_vars([varname_in_file])
117136

118-
# preserve units
137+
# preserve units in regridded dataset
119138
try:
120139
units = ds[varname].units
121140
except Exception as e:
@@ -125,6 +144,15 @@ def load_and_regrid(
125144

126145
ds_regridded[varname] = ds_regridded[varname].assign_attrs({"units": units})
127146

147+
ds_regridded[varname] = ds_regridded[varname].assign_attrs({"units": units})
148+
149+
# time bound check #3 -- preserve time bnds in regridded dataset
150+
if "bounds" in list(ds_regridded.time.attrs.keys()):
151+
time_bnds_key = ds_regridded.time.attrs["bounds"]
152+
if time_bnds_key not in list(ds_regridded.keys()):
153+
ds_regridded = ds_regridded.bounds.add_missing_bounds(["T"])
154+
print("[WARNING]: bounds.add_missing_bounds conducted for T axis")
155+
128156
if debug:
129157
print("ds_regridded:", ds_regridded)
130158
return ds_regridded

‎pcmdi_metrics/mean_climate/lib/mean_climate_metrics_to_json.py

Lines changed: 22 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -22,12 +22,28 @@ def mean_climate_metrics_to_json(
2222
models_in_dict = list(json_dict["RESULTS"].keys())
2323
for m in models_in_dict:
2424
if m == model:
25-
for ref in list(json_dict["RESULTS"][m].keys()):
26-
if ref != "units":
27-
runs_in_model_dict = list(json_dict["RESULTS"][m][ref].keys())
28-
for r in runs_in_model_dict:
29-
if (r != run) and (run is not None):
30-
del json_dict["RESULTS"][m][ref][r]
25+
ref_list = list(json_dict["RESULTS"][m].keys())
26+
if debug:
27+
print(
28+
"debug:: mean_climate_metrics_to_json:: m, ref_list: ",
29+
m,
30+
ref_list,
31+
)
32+
if "units" in ref_list:
33+
ref_list.remove("units")
34+
for ref in ref_list:
35+
if debug:
36+
print(
37+
'debug:: mean_climate_metrics_to_json:: m, ref, json_dict["RESULTS"][m][ref]: ',
38+
m,
39+
ref,
40+
json_dict["RESULTS"][m][ref],
41+
)
42+
runs_in_model_dict = list(json_dict["RESULTS"][m][ref].keys())
43+
for r in runs_in_model_dict:
44+
if (r != run) and (run is not None):
45+
del json_dict["RESULTS"][m][ref][r]
46+
3147
else:
3248
del json_dict["RESULTS"][m]
3349
# Write selected dict to JSON

0 commit comments

Comments
 (0)