-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathdemo.py
215 lines (184 loc) · 8.97 KB
/
demo.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
213
214
215
import ee
from gridwxcomp.prep_metadata import prep_metadata
from gridwxcomp.ee_download import download_grid_data
from gridwxcomp.plot import daily_comparison, monthly_comparison, station_bar_plot
from gridwxcomp.calc_bias_ratios import calc_bias_ratios
from gridwxcomp.spatial import make_grid, interpolate
from gridwxcomp.util import reproject_crs_for_bounds
# name of the dataset comparison is being made with
gridded_dataset_name = 'conus404'
# local path for config file
config_path = f'gridwxcomp/example_data/gridwxcomp_config_{gridded_dataset_name}.ini'
# Path to station metadata file with lat/long coords
station_meta_path = 'gridwxcomp/example_data/Station_Data.txt'
# initialize earth engine
ee.Authenticate()
ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')
# local path for prep_metadata output
gridwxcomp_input = f'{gridded_dataset_name}_gridwxcomp_metadata.csv'
# Directory that bias ratio/interpolation outputs will be saved to
output_dir = f'bias_outputs_{gridded_dataset_name}'
# Defining the parameters for the interpolation
search_radius = 15 * 100 * 1000
params = {
'power': 4,
'smoothing': 0.0,
'max_points': 25,
'min_points': 0,
'radius': search_radius,
'nodata': -999}
# in meters, since interpolation will happen in Lambert Conformal Conic
lcc_interpolation_resolution = 1000
'''
prep_metadata
The purpose of this module is to open the file provided at 'station_meta_path' and verify
it has everything needed in order to proceed with acquiring gridded data and comparing it against
the observed station data
'''
prep_metadata(
station_meta_path,
config_path,
gridded_dataset_name,
out_path=gridwxcomp_input)
'''
define_projection_parameters
Functions in gridwxcomp.spatial and gridwxcomp.interpgdal are expecting projection information in a dictionary
This might seem cumbersome but if this is fleshed out a little more we could have gridwxcomp function with any two
EPSG projections, although the first is required to be the same projection as the input data coordinates
This should eventually be made its own function under prep_metadata.py, where it asks for:
- two projection names (arbitrary, user will refer to them later in other functions)
- two projection resolutions
- optionally, bounds for the first projection, if none provided uses spatial.get_subgrid_bounds(),
which should also be reclassified under prep_metadata.py
'''
projection_dict = {
'wgs84': { # WGS84 is EPSG:4326 and is in decimal degrees
'bounds': {
'xmin': -111.8708332996666428,
'xmax': -108.6208332996662733,
'ymin': 38.0874999999668162,
'ymax': 40.5874999999666741},
'resolution': 0.1,
'crs_id': 'EPSG:4326'},
'lcc': { # Lambert Conformal Conic is ESRI:102004 and is in meters
'bounds': {},
'resolution': 1000,
'crs_id': 'ESRI:102004'}
}
# User has option to automatically pull bounds, see TODO in function over assumptions
# projection_dict['wgs84']['bounds'] = get_subgrid_bounds(gridwxcomp_input, config_path, buffer=25)
projection_dict['lcc']['bounds'] = reproject_crs_for_bounds(
projection_dict['wgs84']['bounds'], projection_dict['lcc']['resolution'],
projection_dict['wgs84']['crs_id'], projection_dict['lcc']['crs_id'], 0)
'''
download_grid_data
The purpose of this module is to make calls to the earth engine API to
export timeseries for the gridded dataset at the station locations
contained within the station metadata file generated by prep_metadata.
The data is saved locally as it is downloaded, and the metadata file
generated by prep_metadata is amended with paths to the downloaded files.
'''
download_grid_data(
gridwxcomp_input,
config_path,
local_folder=None,
force_download=False)
'''
plotting
The purpose of this module is to generate bokeh plots at both the daily and monthly timesteps to visualize
differences/similarities between the station and gridded datasets. Both methods will generate a subplot
for every variable which is shared between the two datasets. Each timeseries plot is accompanied by
a scatterplot featuring a LSR line that is forced through zero, giving a sense of the overall bias and
correlation between the two datasets
The 'monthly_comparison()' method will generate monthly averages to get a better sense of the overall bias
by reducing the noise of daily observations (Ex. observations from Jan 1st - 31st 2020 become one average).
This generates one plot per station included in the input file.
The 'daily_comparison()' will generate 12 plots per station, one for each month of the year. These months are
selected and plotted together, which allows the user to visualize if one particular year has had a substantial
deviation for other years in the record
These plots are mainly diagnostic and are not used as part of the later steps.
The plotting module also contains a function to visualize the results from 'calc_bias_ratios' or 'interpolate'
in the form of bar plots.
'''
monthly_comparison(
gridwxcomp_input,
config_path,
dataset_name=gridded_dataset_name)
daily_comparison(
gridwxcomp_input,
config_path,
dataset_name=gridded_dataset_name)
for var in ['prcp', 'etr', 'eto']: # Iterate over vars in list. Valid entries found in calc_bias_ratios.py VAR_LIST
ratio_filepath = f'{output_dir}/{var}_summary_comp_all_yrs.csv' # path to bias ratios output file
interpolation_out_path = (f'{var}_invdistnn_p{params["power"]}_' # directory for interpolation outputs
f's{params["smoothing"]}_maxpoints{params["max_points"]}_radius{params["radius"]}')
'''
calc_bias_ratios
The purpose of this module is to calculate the monthly bias
(either the long-term-mean [default] or mean-of-annual)
between the gridded dataset and the station observations.
IMPORTANT: The bias factors are either
(station - gridded) [Temperature] or (station / gridded) [Everything else].
This is the reverse relationship from how we normally think of modeled / observed, so a temperature bias
of -2 means that the gridded dataset is hotter than observed, not colder. This reversal was done to make
applying bias a matter of addition/multiplication.
This module requires the configuration .ini file (stored under the example_data folder)
to be set up so that it can correct interpret the columns within the station and gridded data files.
The acceptable parameters for 'comparison_var' are coded in ACCEPTABLE_VAR_LIST within calc_bias_ratios.py
'''
calc_bias_ratios(
gridwxcomp_input,
config_path,
output_dir,
method='long_term_mean',
grid_id_name='GRID_ID',
comparison_var=var,
grid_id=None,
day_limit=10,
years='all',
comp=True)
'''
station_bar_plot
Generates boxplots of the bias ratios to visualize overall performance.
Requires the outputs of calc_bias_ratios as an input.
'''
station_bar_plot(ratio_filepath, bar_plot_layer='grow_mean')
'''
make grid
This will generate a fishnet grid to Make fishnet grid (vector file of polygon geometry)
based on bounding coordinates. Each cell in the grid will be assigned
a unique numerical identifier (property specified by grid_id_name) and the grid will be in
the WGS84 coordinate system.
This step is not required if the desired output is just the interpolation surfaces. It becomes necessary
zonal statistics are going to be generated by spatial.interpolate.
'''
make_grid(
ratio_filepath,
projection_dict['wgs84']['resolution'],
bounds=projection_dict['wgs84']['bounds'],
overwrite=False,
grid_id_name='GRID_ID')
'''
interpolate
Contains several methods to create a 2-D interpolated surface based on the lat/long points and values
generated by calc_bias_ratios.py. Provides options allow for up/down-scaling the resolution of the
resampling grid and to select from multiple interpolation methods.
Interpolation will occur in lambert conformal conic, however the end product
Interploated surfaces are saved as GeoTIFF rasters and will be in the WGS84 CRS.
The default behavior will also calculate zonal statistics (which requires having run spatial.make_grid).
To disable this you must specifically pass z_stats=False
'''
interpolate(
ratio_filepath,
projection_dict,
proj_name='lcc',
layer='all',
out=interpolation_out_path,
scale_factor=1,
function='invdistnn',
params=params,
buffer=5,
z_stats=False,
res_plot=False,
grid_id_name='GRID_ID',
options=None)