Skip to content

Commit eb662f1

Browse files
author
Joe Hamman
committed
Changes to allow for routing on the same grid as the target domain.
i.e. if REMAP = False and AGGREGATE = False.
1 parent 9b6ecb6 commit eb662f1

File tree

5 files changed

+112
-60
lines changed

5 files changed

+112
-60
lines changed

‎make_parameters.py‎

Lines changed: 30 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -64,7 +64,6 @@ def main():
6464

6565
outlets = OrderedDict(sorted(results.items(), key=lambda t: t[0]))
6666
else:
67-
print(outlets)
6867
for name, outlet in outlets.iteritems():
6968
outlet = gen_uh_run(uh_box, fdr_data, fdr_vatts, dom_data, outlet,
7069
config_dict, directories)
@@ -167,7 +166,7 @@ def gen_uh_init(config_file):
167166

168167
# ---------------------------------------------------------------- #
169168
# Check latitude order, flip if necessary.
170-
if fdr_data[fdr_lat][-1] > fdr_data[fdr_lat][0]:
169+
if fdr_data[fdr_lat][-1] < fdr_data[fdr_lat][0]:
171170
log.debug('Inputs came in upside down, flipping everything now.')
172171
var_list = fdr_data.keys()
173172
var_list.remove(fdr_lon)
@@ -210,9 +209,9 @@ def gen_uh_init(config_file):
210209
# ---------------------------------------------------------------- #
211210
# If remap is False, domain coordinates needs to be in the fdr coordinates
212211
# We can move the unit hydrographs to the domain grid later
213-
if not options['REMAP']:
214-
log.error('RVIC parameter generation requires REMAP option to be True')
215-
log.error('In theory, this is possible but it has not been tested.')
212+
if options['AGGREGATE'] and not options['REMAP']:
213+
log.error('RVIC parameter generation requires REMAP option to be True'
214+
' if AGGREGATE is True')
216215
raise ValueError('Invalid option')
217216
# ---------------------------------------------------------------- #
218217

@@ -257,7 +256,7 @@ def gen_uh_init(config_file):
257256
# ---------------------------------------------------------------- #
258257

259258
log.info('Finished with gen_uh_init')
260-
log.info('------------------------------------------------------------\n')
259+
log.info('-------------------------------------------------------------\n')
261260

262261
return (uh_box, fdr_data, fdr_vatts, dom_data, outlets,
263262
config_dict, directories)
@@ -351,17 +350,38 @@ def gen_uh_run(uh_box, fdr_data, fdr_vatts, dom_data, outlet, config_dict,
351350
clean_file(temp_file_2)
352351
# -------------------------------------------------------- #
353352

353+
# -------------------------------------------------------- #
354+
# Set the domain index offest to zero and use the remapped fraction
355+
# as the final fractions
356+
y0, x0 = 0, 0
357+
final_fracs = final_data['fraction']
358+
# -------------------------------------------------------- #
354359
else:
360+
# -------------------------------------------------------- #
361+
# Put the agg data back onto the original grid
355362
final_data = agg_data
363+
final_fracs = np.zeros_like(fdr_data['velocity'],
364+
dtype=np.float64)
365+
x0 = final_data['dom_x_min']
366+
x1 = final_data['dom_x_max']
367+
y0 = final_data['dom_y_min']
368+
y1 = final_data['dom_y_max']
369+
final_fracs[y0:y1, x0:x1] = final_data['fraction']
370+
# -------------------------------------------------------- #
356371
# ------------------------------------------------------------ #
357372

358373
# ------------------------------------------------------------ #
359-
# Add to adjust fractions Structure
360-
y, x = np.nonzero((final_data['fraction'] > 0.0) *
374+
# Add to "adjust fractions structure"
375+
y, x = np.nonzero((final_fracs > 0.0) *
361376
(dom_data[config_dict['DOMAIN']['LAND_MASK_VAR']] == 1))
377+
yi = y - y0
378+
xi = x - x0
379+
380+
# From final data
381+
outlet.fractions = final_data['fraction'][yi, xi]
382+
outlet.unit_hydrograph = final_data['unit_hydrograph'][:, yi, xi]
362383

363-
outlet.fractions = final_data['fraction'][y, x]
364-
outlet.unit_hydrograph = final_data['unit_hydrograph'][:, y, x]
384+
# From domain data
365385
outlet.time = np.arange(final_data['unit_hydrograph'].shape[0])
366386
outlet.lon_source = dom_data[config_dict['DOMAIN']['LONGITUDE_VAR']][y, x]
367387
outlet.lat_source = dom_data[config_dict['DOMAIN']['LATITUDE_VAR']][y, x]

‎rvic/aggregate.py‎

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,9 @@ def make_agg_pairs(pour_points, dom_data, fdr_data, config_dict):
8383

8484
# ---------------------------------------------------------------- #
8585
# Sort based on outlet total source area pour_point.source_area
86-
outlets = OrderedDict(sorted(outlets.items(), key=lambda t: t[1].upstream_area, reverse=True))
86+
outlets = OrderedDict(sorted(outlets.items(),
87+
key=lambda t: t[1].upstream_area,
88+
reverse=True))
8789
# ---------------------------------------------------------------- #
8890

8991
# ---------------------------------------------------------------- #
@@ -104,7 +106,7 @@ def make_agg_pairs(pour_points, dom_data, fdr_data, config_dict):
104106
log.info('NUMBER OF POUR POINTS AGGREGATED: %i' % pp_count)
105107
log.info('EFFECIENCY OF: %.2f %%' % (100.*pp_count / num))
106108
log.info('UNASSIGNED POUR POINTS: %i \n' % (num-pp_count))
107-
log.info('---------------------------------------------------------------\n')
109+
log.info('-------------------------------------------------------------\n')
108110

109111
# ---------------------------------------------------------------- #
110112
return outlets
@@ -143,7 +145,8 @@ def aggregate(in_data, agg_data, res=0, pad=0, maskandnorm=False):
143145
lons = np.arange((lon_min-res*pad), (lon_max+res*(pad+1)), res)
144146

145147
fraction = np.zeros((len(lats), len(lons)), dtype=np.float64)
146-
unit_hydrograph = np.zeros((tshape, len(lats), len(lons)), dtype=np.float64)
148+
unit_hydrograph = np.zeros((tshape, len(lats), len(lons)),
149+
dtype=np.float64)
147150
log.debug('fraction shape %s' % str(fraction.shape))
148151
log.debug('unit_hydrograph shape %s' % str(unit_hydrograph.shape))
149152
# ---------------------------------------------------------------- #
@@ -156,15 +159,17 @@ def aggregate(in_data, agg_data, res=0, pad=0, maskandnorm=False):
156159
ilon_min_ind = find_nearest(lons, np.min(in_data['lon']))
157160
ilon_max_ind = find_nearest(lons, np.max(in_data['lon']))+1
158161

159-
log.debug('in_data fraction shape: %s' % str(in_data['fraction'].shape))
162+
log.debug('in_data fraction shape: %s',
163+
str(in_data['fraction'].shape))
160164

161165
if agg_data:
162166
alat_min_ind = find_nearest(lats, np.max(agg_data['lat']))
163167
alat_max_ind = find_nearest(lats, np.min(agg_data['lat']))+1
164168
alon_min_ind = find_nearest(lons, np.min(agg_data['lon']))
165169
alon_max_ind = find_nearest(lons, np.max(agg_data['lon']))+1
166170

167-
log.debug('agg_data fraction shape: %s' % str(agg_data['fraction'].shape))
171+
log.debug('agg_data fraction shape: %s',
172+
str(agg_data['fraction'].shape))
168173

169174
# ---------------------------------------------------------------- #
170175

@@ -186,7 +191,8 @@ def aggregate(in_data, agg_data, res=0, pad=0, maskandnorm=False):
186191
yv, xv = np.nonzero(fraction > 0.0)
187192
unit_hydrograph[:, yv, xv] /= unit_hydrograph[:, yv, xv].sum(axis=0)
188193

189-
# Mask the unit_hydrograph and make sure they sum to 1 at each grid cell
194+
# Mask the unit_hydrograph and make sure they sum to 1 at each grid
195+
# cell
190196
ym, xm = np.nonzero(fraction <= 0.0)
191197
unit_hydrograph[:, ym, xm] = FILLVALUE_F
192198

@@ -195,7 +201,6 @@ def aggregate(in_data, agg_data, res=0, pad=0, maskandnorm=False):
195201

196202
# ---------------------------------------------------------------- #
197203
# Put all the data into agg_data variable and return to main
198-
199204
agg_data['timesteps'] = in_data['timesteps']
200205
agg_data['unit_hydrograph_dt'] = in_data['unit_hydrograph_dt']
201206
agg_data['lon'] = lons

‎rvic/convert.py‎

Lines changed: 16 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
log = logging.getLogger(LOG_NAME)
1515
# -------------------------------------------------------------------- #
1616

17+
1718
# -------------------------------------------------------------------- #
1819
# read station file
1920
def read_station_file(file_name, dom_data, config_dict):
@@ -38,17 +39,17 @@ def read_station_file(file_name, dom_data, config_dict):
3839
x = int(x)-1
3940

4041
# make sure files exist
41-
log.info('On station: %s, active: %s' %(name, active))
42+
log.info('On station: %s, active: %s', name, active)
43+
uhs2_file = os.path.join(config_dict['UHS_FILES']['ROUT_DIR'],
44+
name+'.uh_s2')
4245
if active == '1':
4346
if os.path.isfile(uhs_file):
4447
active = True
45-
elif os.path.isfile(os.path.join(config_dict['UHS_FILES']['ROUT_DIR'],
46-
name+'.uh_s2')):
47-
uhs_file = os.path.join(config_dict['UHS_FILES']['ROUT_DIR'],
48-
name+'.uh_s2')
48+
elif os.path.isfile(uhs2_file):
49+
active = True
4950
else:
50-
raise ValueError('missing uhs_file: (%s or %s)' %(uhs_file, os.path.join(config_dict['UHS_FILES']['ROUT_DIR'],
51-
name+'.uh_s2')))
51+
raise ValueError('missing uhs_file: (%s or %s)', uhs_file,
52+
uhs2_file)
5253
else:
5354
active = False
5455

@@ -62,11 +63,12 @@ def read_station_file(file_name, dom_data, config_dict):
6263
outlets[i].lon = dom_data[config_dict['DOMAIN']['LONGITUDE_VAR']][y, x]
6364
outlets[i].lat = dom_data[config_dict['DOMAIN']['LATITUDE_VAR']][y, x]
6465
else:
65-
log.info('%s not active... skipping' %name)
66+
log.info('%s not active... skipping', name)
6667
f.close()
6768
return outlets
6869
# -------------------------------------------------------------------- #
6970

71+
7072
# -------------------------------------------------------------------- #
7173
# Read uhs files
7274
def read_uhs_files(outlets, dom_data, config_dict):
@@ -81,17 +83,17 @@ def read_uhs_files(outlets, dom_data, config_dict):
8183
"""
8284
if config_dict['UHS_FILES']['ROUT_PROGRAM'] == 'C':
8385
for cell_id, outlet in outlets.iteritems():
84-
log.info('Reading outlet %i: %s' %(cell_id, outlet.name))
86+
log.info('Reading outlet %i: %s', cell_id, outlet.name)
8587
log.debug(outlet.uhs_file)
8688
f = open(outlet.uhs_file, 'r')
8789
num_sources = int(f.readline())
88-
log.debug('Number of sources in file: %i' %num_sources)
90+
log.debug('Number of sources in file: %i', num_sources)
8991
# setup some empty arrays
9092
outlets[cell_id].lon_source = np.empty(num_sources)
9193
outlets[cell_id].lat_source = np.empty(num_sources)
9294
outlets[cell_id].fractions = np.empty(num_sources)
93-
outlets[cell_id].x_source = np.empty(num_sources, dtype=int)
94-
outlets[cell_id].y_source = np.empty(num_sources, dtype=int)
95+
outlets[cell_id].x_source = np.empty(num_sources, dtype=np.int16)
96+
outlets[cell_id].y_source = np.empty(num_sources, dtype=np.int16)
9597

9698
uh = []
9799

@@ -107,7 +109,7 @@ def read_uhs_files(outlets, dom_data, config_dict):
107109
outlets[cell_id].fractions[j] = float(fracs)
108110
outlets[cell_id].x_source[j] = int(x)
109111
outlets[cell_id].y_source[j] = int(y)
110-
line = re.sub(' +',' ',f.readline())
112+
line = re.sub(' +', ' ', f.readline())
111113
uh.append(map(float, line.split()))
112114

113115
outlets[cell_id].unit_hydrograph = np.rot90(np.array(uh,
@@ -159,6 +161,7 @@ def read_uhs_files(outlets, dom_data, config_dict):
159161
return outlets
160162
# -------------------------------------------------------------------- #
161163

164+
162165
# -------------------------------------------------------------------- #
163166
# Adjust Domain Bounds
164167
def move_domain(dom_data, new_dom_data, outlets):

‎rvic/make_uh.py‎

Lines changed: 21 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -40,9 +40,9 @@ def rout(pour_point, uh_box, fdr_data, fdr_atts, rout_dict):
4040
# Find Basin Dims and ID
4141
basin_id = fdr_data[rout_dict['BASIN_ID_VAR']][pour_point.routy, pour_point.routx]
4242

43-
log.info('Input Latitude: %f' % pour_point.lat)
44-
log.info('Input Longitude: %f' % pour_point.lon)
45-
log.info('Global Basid ID: %i' % basin_id)
43+
log.info('Input Latitude: %f', pour_point.lat)
44+
log.info('Input Longitude: %f', pour_point.lon)
45+
log.info('Global Basid ID: %i', basin_id)
4646

4747
y_inds, x_inds = np.nonzero(fdr_data[rout_dict['BASIN_ID_VAR']] == basin_id)
4848
y = np.arange(len(fdr_data[rout_dict['LATITUDE_VAR']]))
@@ -65,7 +65,7 @@ def rout(pour_point, uh_box, fdr_data, fdr_atts, rout_dict):
6565
basin['velocity'] = fdr_data['velocity'][y_min:y_max, x_min:x_max]
6666
basin['diffusion'] = fdr_data['diffusion'][y_min:y_max, x_min:x_max]
6767

68-
log.debug('Grid cells in subset: %i' % basin['velocity'].size)
68+
log.debug('Grid cells in subset: %i', basin['velocity'].size)
6969

7070
pour_point.basiny, pour_point.basinx = latlon2yx(plats=pour_point.lat,
7171
plons=pour_point.lon,
@@ -75,8 +75,14 @@ def rout(pour_point, uh_box, fdr_data, fdr_atts, rout_dict):
7575

7676
# ---------------------------------------------------------------- #
7777
# Create the rout_data Dictionary
78-
rout_data = {'lat': basin['lat'], 'lon': basin['lon']}
79-
78+
rout_data = {'lat': basin['lat'], 'lon': basin['lon'],
79+
'dom_x_min': x_min, 'dom_y_min': y_min,
80+
'dom_x_max': x_max, 'dom_y_max': y_max}
81+
log.debug('Clipping indicies:')
82+
log.debug('dom_x_min: %s', rout_data['dom_x_min'])
83+
log.debug('dom_x_max: %s', rout_data['dom_x_max'])
84+
log.debug('dom_y_min: %s', rout_data['dom_y_min'])
85+
log.debug('dom_y_max: %s', rout_data['dom_y_max'])
8086
# ---------------------------------------------------------------- #
8187
# Determine low direction syntax
8288
if 'VIC' in fdr_atts[rout_dict['FLOW_DIRECTION_VAR']]:
@@ -132,12 +138,12 @@ def rout(pour_point, uh_box, fdr_data, fdr_atts, rout_dict):
132138

133139
# ---------------------------------------------------------------- #
134140
# Agregate to output timestep
135-
rout_data['unit_hydrograph'], rout_data['timesteps'] = adjust_uh_timestep(uh_s, t_uh,
136-
input_interval,
137-
rout_dict['OUTPUT_INTERVAL'],
138-
catchment['x_inds'],
139-
catchment['y_inds'])
140-
print ">>>>>>>>unit_hydrograph", rout_data['unit_hydrograph'].sum()
141+
rout_data['unit_hydrograph'], \
142+
rout_data['timesteps'] = adjust_uh_timestep(uh_s, t_uh,
143+
input_interval,
144+
rout_dict['OUTPUT_INTERVAL'],
145+
catchment['x_inds'],
146+
catchment['y_inds'])
141147
# ---------------------------------------------------------------- #
142148
return rout_data
143149
# -------------------------------------------------------------------- #
@@ -150,7 +156,7 @@ def find_ts(uh_t):
150156
Determines the (input_interval) based on the timestep given in uhfile
151157
"""
152158
input_interval = uh_t[1]-uh_t[0]
153-
log.debug('Input Timestep = %i seconds' % input_interval)
159+
log.debug('Input Timestep = %i seconds', input_interval)
154160
return input_interval
155161
# -------------------------------------------------------------------- #
156162

@@ -170,8 +176,8 @@ def read_direction(fdr, dy, dx):
170176

171177
for (y, x), d in np.ndenumerate(fdr):
172178
try:
173-
to_y[y, x] = y+dy[d]
174-
to_x[y, x] = x+dx[d]
179+
to_y[y, x] = y - dy[d]
180+
to_x[y, x] = x + dx[d]
175181
except KeyError as e:
176182
if (d == 0) or (d == -9999):
177183
to_y[y, x] = -9999

‎rvic/utilities.py‎

Lines changed: 33 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -275,41 +275,59 @@ def read_domain(domain_dict):
275275

276276
# ---------------------------------------------------------------- #
277277
# Create the cell_ids variable
278-
dom_data['cell_ids'] = np.arange(dom_data[domain_dict['LAND_MASK_VAR']].size).reshape(dom_data[domain_dict['LAND_MASK_VAR']].shape)
278+
dom_mask = domain_dict['LAND_MASK_VAR']
279+
temp = np.arange(dom_data[dom_mask].size)
280+
dom_data['cell_ids'] = temp.reshape(dom_data[dom_mask].shape)
279281
# ---------------------------------------------------------------- #
280282

281283
# ---------------------------------------------------------------- #
282284
# Make sure the longitude / latitude vars are 2d
283-
dom_data['cord_lons'] = dom_data[domain_dict['LONGITUDE_VAR']]
284-
dom_data['cord_lats'] = dom_data[domain_dict['LATITUDE_VAR']]
285-
if dom_data[domain_dict['LONGITUDE_VAR']].ndim == 1:
286-
dom_data[domain_dict['LONGITUDE_VAR']], \
287-
dom_data[domain_dict['LATITUDE_VAR']] = np.meshgrid(dom_data[domain_dict['LONGITUDE_VAR']],
288-
dom_data[domain_dict['LATITUDE_VAR']])
285+
dom_lat = domain_dict['LATITUDE_VAR']
286+
dom_lon = domain_dict['LONGITUDE_VAR']
287+
if dom_data[dom_lon].ndim == 1:
288+
# ------------------------------------------------------------- #
289+
# Check latitude order, flip if necessary.
290+
if dom_data[dom_lat][-1] < dom_data[dom_lat][0]:
291+
log.debug('Inputs came in upside down, flipping everything now.')
292+
var_list = dom_data.keys()
293+
var_list.remove(dom_lon)
294+
for var in var_list:
295+
dom_data[var] = np.flipud(dom_data[var])
296+
# ------------------------------------------------------------ #
289297

298+
# ------------------------------------------------------------- #
299+
# Make 2d coordinate vars
300+
dom_data[dom_lon], dom_data[dom_lat] = np.meshgrid(dom_data[dom_lon],
301+
dom_data[dom_lat])
302+
# ------------------------------------------------------------- #
303+
dom_data['cord_lons'] = dom_data[dom_lon]
304+
dom_data['cord_lats'] = dom_data[dom_lat]
290305
# ---------------------------------------------------------------- #
291306

292307
# ---------------------------------------------------------------- #
293308
# Make sure the area is in m2
294-
area_units = dom_vatts[domain_dict['AREA_VAR']]['units']
309+
dom_area = domain_dict['AREA_VAR']
310+
area_units = dom_vatts[dom_area]['units']
295311

296312
if area_units in ["rad2", "radians2", "radian2", "radian^2", "rad^2",
297-
"radians^2", "rads^2", "radians squared", "square-radians"]:
298-
dom_data[domain_dict['AREA_VAR']] = dom_data[domain_dict['AREA_VAR']]*EARTHRADIUS*EARTHRADIUS
313+
"radians^2", "rads^2", "radians squared",
314+
"square-radians"]:
315+
dom_data[dom_area] = dom_data[dom_area]*EARTHRADIUS*EARTHRADIUS
299316
elif area_units in ["m2", "m^2", "meters^2", "meters2", "square-meters",
300317
"meters squared"]:
301-
dom_data[domain_dict['AREA_VAR']] = dom_data[domain_dict['AREA_VAR']]
318+
dom_data[dom_area] = dom_data[dom_area]
302319
elif area_units in ["km2", "km^2", "kilometers^2", "kilometers2",
303320
"square-kilometers", "kilometers squared"]:
304-
dom_data[domain_dict['AREA_VAR']] = dom_data[domain_dict['AREA_VAR']]*METERSPERKM*METERSPERKM
321+
dom_data[dom_area] = dom_data[dom_area]*METERSPERKM*METERSPERKM
305322
elif area_units in ["mi2", "mi^2", "miles^2", "miles", "square-miles",
306323
"miles squared"]:
307-
dom_data[domain_dict['AREA_VAR']] = dom_data[domain_dict['AREA_VAR']]*METERSPERMILE*METERSPERMILE
324+
dom_data[dom_area] = dom_data[dom_area]*METERSPERMILE*METERSPERMILE
308325
elif area_units in ["acres", "ac", "ac."]:
309-
dom_data[domain_dict['AREA_VAR']] = dom_data[domain_dict['AREA_VAR']]*METERS2PERACRE
326+
dom_data[dom_area] = dom_data[dom_area]*METERS2PERACRE
310327
else:
311328
log.warning("WARNING: UNKNOWN AREA units (%s), ASSUMING THEY ARE IN "
312-
"SQUARE METERS" % dom_data[domain_dict['AREA_VAR']]['units'])
329+
"SQUARE METERS",
330+
dom_data[domain_dict['AREA_VAR']]['units'])
313331
# ---------------------------------------------------------------- #
314332
return dom_data, dom_vatts, dom_gatts
315333
# -------------------------------------------------------------------- #

0 commit comments

Comments
 (0)