-
Notifications
You must be signed in to change notification settings - Fork 37
/
std_of_meandiurnalcycle.py
executable file
·212 lines (178 loc) · 6.82 KB
/
std_of_meandiurnalcycle.py
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
#!/usr/bin/env python
# Globally average the variance of CMIP5/AMIP mean-diurnal-cycle precipitation from a given month and range of years.
# Square the standard deviatios before averaging over area, and then take the square root. This is the correct order of
# operations to get the fraction of total variance, per Covey and Gehne
# 2016 (https://e-reports-ext.llnl.gov/pdf/823104.pdf).
# This has the PMP Parser "wrapped" around it, so it's executed with both input and output parameters specified in the
# Unix command line.
# Charles Doutriaux September 2017
# Curt Covey (from ./old_meandiurnalcycWrappedInOut.py) July 2017
from __future__ import print_function
import collections
import glob
import json
import multiprocessing as mp
import os
import cdms2
import cdp
import cdutil
import genutil
import pcmdi_metrics
from pcmdi_metrics import resources
from pcmdi_metrics.diurnal.common import (
INPUT,
P,
monthname_d,
populateStringConstructor,
)
def main():
def compute(param):
template = populateStringConstructor(args.filename_template, args)
template.variable = param.varname
template.month = param.monthname
fnameRoot = param.fileName
reverted = template.reverse(os.path.basename(fnameRoot))
model = reverted["model"]
print("Specifying latitude / longitude domain of interest ...")
datanameID = "diurnalmean" # Short ID name of output data
latrange = (param.args.lat1, param.args.lat2)
lonrange = (param.args.lon1, param.args.lon2)
region = cdutil.region.domain(latitude=latrange, longitude=lonrange)
if param.args.region_name == "":
region_name = "{:g}_{:g}&{:g}_{:g}".format(*(latrange + lonrange))
else:
region_name = param.args.region_name
print("Reading %s ..." % fnameRoot)
try:
f = cdms2.open(fnameRoot)
x = f(datanameID, region)
units = x.units
print(" Shape =", x.shape)
print("Finding standard deviation over first dimension (time of day) ...")
x = genutil.statistics.std(x)
print(" Shape =", x.shape)
print("Finding r.m.s. average over 2nd-3rd dimensions (area) ...")
x = x * x
x = cdutil.averager(x, axis="xy")
x = cdms2.MV2.sqrt(x)
print(
"For %8s in %s, average variance of hourly values = (%5.2f %s)^2"
% (model, monthname, x, units)
)
f.close()
except Exception as err:
print("Failed model %s with error" % (err))
x = 1.0e20
return model, region, {region_name: float(x)}
P.add_argument(
"-j",
"--outnamejson",
type=str,
dest="outnamejson",
default="pr_%(month)_%(firstyear)-%(lastyear)_std_of_meandiurnalcyc.json",
help="Output name for jsons",
)
P.add_argument("--lat1", type=float, default=-50.0, help="First latitude")
P.add_argument("--lat2", type=float, default=50.0, help="Last latitude")
P.add_argument("--lon1", type=float, default=0.0, help="First longitude")
P.add_argument("--lon2", type=float, default=360.0, help="Last longitude")
P.add_argument(
"--region_name",
type=str,
default="TRMM",
help="name for the region of interest",
)
P.add_argument(
"-t",
"--filename_template",
default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_avg.nc",
)
P.add_argument("--model", default="*")
P.add_argument(
"--cmec",
dest="cmec",
action="store_true",
default=False,
help="Use to save metrics in CMEC JSON format",
)
P.add_argument(
"--no_cmec",
dest="cmec",
action="store_false",
default=False,
help="Use to disable saving metrics in CMEC JSON format",
)
args = P.get_parameter()
month = args.month
monthname = monthname_d[month]
startyear = args.firstyear # noqa: F841
finalyear = args.lastyear # noqa: F841
cmec = args.cmec
template = populateStringConstructor(args.filename_template, args)
template.month = monthname
print("TEMPLATE NAME:", template())
print("Specifying latitude / longitude domain of interest ...")
# TRMM (observed) domain:
latrange = (args.lat1, args.lat2)
lonrange = (args.lon1, args.lon2)
region = cdutil.region.domain(latitude=latrange, longitude=lonrange)
# Amazon basin:
# latrange = (-15.0, -5.0)
# lonrange = (285.0, 295.0)
print("Preparing to write output to JSON file ...")
os.makedirs(args.results_dir, exist_ok=True)
jsonFile = populateStringConstructor(args.outnamejson, args)
jsonFile.month = monthname
jsonname = os.path.join(os.path.abspath(args.results_dir), jsonFile())
if not os.path.exists(jsonname) or args.append is False:
print("Initializing dictionary of statistical results ...")
stats_dic = {}
metrics_dictionary = collections.OrderedDict()
else:
with open(jsonname) as f:
metrics_dictionary = json.load(f)
print("LOADE WITH KEYS:", list(metrics_dictionary.keys()))
stats_dic = metrics_dictionary["RESULTS"]
OUT = pcmdi_metrics.io.base.Base(os.path.abspath(args.results_dir), jsonFile())
egg_pth = resources.resource_path()
disclaimer = open(os.path.join(egg_pth, "disclaimer.txt")).read()
metrics_dictionary["DISCLAIMER"] = disclaimer
metrics_dictionary["REFERENCE"] = (
"The statistics in this file are based on Trenberth, Zhang & Gehne, "
"J Hydromet. 2017"
)
files = glob.glob(os.path.join(args.modpath, template()))
print(files)
params = [INPUT(args, name, template) for name in files]
print("PARAMS:", params)
results = cdp.cdp_run.multiprocess(compute, params, num_workers=args.num_workers)
for r in results:
m, region, res = r
if r[0] not in stats_dic:
stats_dic[m] = res
else:
stats_dic[m].update(res)
print("Writing output to JSON file ...")
metrics_dictionary["RESULTS"] = stats_dic
print("KEYS AT END:", list(metrics_dictionary.keys()))
rgmsk = metrics_dictionary.get("RegionalMasking", {})
print("REG MASK:", rgmsk)
nm = list(res.keys())[0]
region.id = nm
rgmsk[nm] = {"id": nm, "domain": region}
metrics_dictionary["RegionalMasking"] = rgmsk
OUT.write(
metrics_dictionary,
json_structure=["model", "domain"],
indent=4,
separators=(",", ": "),
)
if cmec:
print("Writing cmec file")
OUT.write_cmec(indent=4, separators=(",", ": "))
print("done")
# Good practice to place contents of script under this check
if __name__ == "__main__":
# Related to script being installed as executable
mp.freeze_support()
main()