-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathasop_duration_example.py
163 lines (136 loc) · 6.94 KB
/
asop_duration_example.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
import asop_duration as asop
import numpy as np
"""
Example use of ASoP Duration package to compute
and plot diagnostics of the histogram for
preciptiaiton duration in a dataset.
This example uses TRMM 3B42v7A and CMORPH v1.0
data at a common grid of N96.
Written by Amulya Chevuturi
(C) The author 2020
"""
def get_dictionary(key):
"""
The Coherence package relies on a dataset dictionary,
for which currently the user must specify most values.
This function shows how to build a dictionary. The dictionary to
be returned to the main program is selected through the use of
a "key", here either "TRMM" or "CMORPH"
Arguments:
* key:
A string used to select the correct dataset
Returns:
* asop_dict:
A dictionary containing required and optional
parameters for the ASoP Coherence package.
Required dictionary keys:
infile - Path to input file (must be readable by Iris)
name - Name for model in plot files (no spaces)
legend_name - Name for model in legends and titles on plots (spaces allowed)
dt - Timestep of input data (in seconds)
dx - Longitudinal grid spacing at equator (in km)
dy - Latitudinal grid spacing (in km)
constraint - standard_name of data to read from netCDF file (e.g., precipitation flux)
scale_factor - Multiplier necessary to convert precipitation to units of mm/day
region - Region of data to read [minlat, maxlat, minlon, maxlon]
box_size - Length of sub-regions (square boxes) to consider for correlation analysis
as a function of physical distance
(in km, recommended value is > 6.5*dx).
color - Name of color to use on line graphs (must be recognised by matplotlib).
region_size - Length of sub-regions (square boxes) to consider for correlation analysis
as a function of model gridpoints (including lag correlations, see below)
(in units of native gridpoints, odd integers strongly recommended).
lag_length - Maximum lag to consider for correlation analysis as a function of model
gridpoints, for constructing distance vs. lag correlation diagrams.
Correlations for lags from 0 (coincidence) until lag_length will be computed.
autocorr_length - The maximum autocorrelation to analyse (in seconds).
Optional dictionary keys, which are useful mainly if analysing the same
dataset on more than one grid / temporal sampling interval:
grid_type - A string describing the grid, used in output filenames
(recommend no spaces).
time_type - A string describing the temporal sampling, used in output filenames
(recommend no spaces).
grid_desc - A string describing the grid, used in plot titles (can contain spaces).
time_desc - A string describing the temporal sampling, used in plot titles
(can contain spaces).
"""
asop_dict = {}
if key == 'TRMM':
asop_dict['infile'] = '/gws/nopw/j04/klingaman/amulya/data/obs/TRMM/TRMM_3B42_V7_0.25deg-DAILY_DJF_1998-2018_brazil_n48.nc'
asop_dict['name'] = 'TRMM'
asop_dict['dt'] = 10800
asop_dict['dx'] = 27
asop_dict['dy'] = 27
asop_dict['constraint'] = 'precipitation_flux'
asop_dict['scale_factor'] = 1.0
asop_dict['legend_name'] = 'TRMM'
asop_dict['region'] = [-10,5,288,310]
asop_dict['box_size'] = 1680
asop_dict['color'] = 'red'
asop_dict['region_size'] = 7
asop_dict['lag_length'] = 6
asop_dict['grid_type'] = 'N48'
asop_dict['time_type'] = 'Day'
asop_dict['grid_desc'] = 'N48'
asop_dict['time_desc'] = 'Day'
asop_dict['autocorr_length'] = 60*60*24
elif key == 'CMORPH':
asop_dict['infile'] = '/gws/nopw/j04/klingaman/amulya/data/obs/CMORPH/CMORPH_V1.0_0.25deg-DAILY_DJF_1998-2017_brazil_n48.nc'
asop_dict['name'] = 'CMORPH'
asop_dict['dt'] = 10800
asop_dict['dx'] = 27
asop_dict['dy'] = 27
asop_dict['constraint'] = 'precipitation_flux'
asop_dict['scale_factor'] = 1.0
asop_dict['legend_name'] = 'CMORPH'
asop_dict['region'] = [-10,5,288,310]
asop_dict['box_size'] = 1680
asop_dict['color'] = 'blue'
asop_dict['region_size'] = 7
asop_dict['lag_length'] = 6
asop_dict['grid_type'] = 'N48'
asop_dict['time_type'] = 'Day'
asop_dict['grid_desc'] = 'N48'
asop_dict['time_desc'] = 'Day'
asop_dict['autocorr_length'] = 60*60*24
return(asop_dict)
if __name__ == '__main__':
datasets = ('TRMM','CMORPH')
n_datasets = len(datasets)
# Allocate memory for multi-model fields
max_box_distance,max_timesteps,max_boxes = asop.parameters()
all_distance_correlations = np.zeros((n_datasets,max_box_distance))
all_distance_ranges = np.zeros((n_datasets,3,max_box_distance))
all_distance_max = np.zeros((n_datasets),dtype=np.int)
all_time_correlations = np.zeros((n_datasets,max_timesteps))
all_time_max = np.zeros((n_datasets),dtype=np.int)
all_dt = np.zeros((n_datasets),dtype=np.int)
all_colors = []
all_legend_names = []
hist_1d = [0]*n_datasets
for i in xrange(n_datasets):
print '--> '+datasets[i]
asop_dict = get_dictionary(datasets[i])
# Read precipitation data
precip = asop.read_precip(asop_dict).data
#convert the iris cube to a numpy array
#precip = prp.data
# Define the precipitation bins as prbins
# Define the duration length bins as cobins
# Note that on plots, the first and last edges will be
# replaced by < and > signs, respectively.
prbins=[0,1,2,4,6,9,12,16,20,25,30,40,60,90,130,180,2e20]
cobins=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,2e20]
thbins=[0,0.1]
# Compute 2D precipitation-duration histograms
hist_2d = asop.compute_twod_histogram(precip,prbins,cobins)
# Plot 2D histograms
asop.plot_twod_histogram(hist_2d,asop_dict,prbins,cobins,title=True,colorbar=True)
# Compute 2D precipitation-duration histogram for a threshold bin
hist_1d[i] = asop.compute_oned_histogram(precip,thbins,cobins)
# Save color and legend information
all_colors.append(asop_dict['color'])
all_legend_names.append(asop_dict['legend_name'])
# Plot 1D histogram
asop.plot_oned_histogram(hist_1d,asop_dict,thbins,cobins,all_colors,all_legend_names,title=True,legend=True)