Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

suggested fixes #82

Merged
merged 4 commits into from
Sep 23, 2024
Merged

suggested fixes #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a57b3a6
Specifying raster file to zonal_stats in setup_basic_stress_data() le…
apryet Apr 4, 2024
1d02f71
added simulate_high_k_lakes condition to setup_array in sourcedata
apryet Jun 25, 2024
cc70d27
fix(sourcedata.py): when setting up the bottom array(s), only make ed…
aleaf Jul 25, 2024
0be02d4
refactor(discretization.fill_cells_vertically): only return bottom ar…
aleaf Sep 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion mfsetup/bcs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,11 @@ def setup_basic_stress_data(model, shapefile=None, csvfile=None,
if meta['transform'][0] > m.modelgrid.delr[0]:
all_touched = True
stat = entry['stat']
results = zonal_stats(polygons, filename, stats=stat,
# load raster and specify affine transform to avoid issues with proj
with rasterio.open(filename) as src:
affine = src.transform
array = src.read(1)
results = zonal_stats(polygons, array, affine=affine, stats=stat,
all_touched=all_touched)
#values = np.ones((m.nrow * m.ncol), dtype=float) * np.nan
#values[cells_with_bc] = np.array([r[stat] for r in results])
Expand Down
6 changes: 3 additions & 3 deletions mfsetup/discretization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,7 @@ def fill_cells_vertically(top, botm):
filled += np.nanmin(all_surfaces, axis=0) # botm[-1]
# append the model bottom elevations
filled = np.append(filled, [np.nanmin(all_surfaces, axis=0)], axis=0)
return filled[0].copy(), filled[1:].copy()
return filled[1:].copy()


def fix_model_layer_conflicts(top_array, botm_array,
Expand Down Expand Up @@ -748,6 +748,6 @@ def voxels_to_layers(voxel_array, z_edges, model_top=None, model_botm=None, no_d

# finally, fill any remaining nans with next layer elevation (going upward)
# might still have nans in areas where there are no voxel values, but model top and botm values
top, botm = fill_cells_vertically(layers[0], layers[1:])
layers = np.vstack([np.reshape(top, (1, *top.shape)), botm])
botm = fill_cells_vertically(layers[0], layers[1:])
layers = np.vstack([np.reshape(layers[0], (1, *layers[0].shape)), botm])
return layers
63 changes: 33 additions & 30 deletions mfsetup/sourcedata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1382,8 +1382,13 @@ def setup_array(model, package, var, data=None,
# (for lakes)
simulate_high_k_lakes = model.cfg['high_k_lakes']['simulate_high_k_lakes']
if var == 'botm':
# only execute this code if building the model (not loading)
if not model._load:
if 'lak' in model.cfg['model']['packages'] and not model._load and\
np.sum(model.lake_bathymetry) != 0:
# only execute this code if building the model (not loading)
# and if lake bathymetry was supplied
# (bathymetry is loaded from the configuration file input
# when the model.lake_bathymetry property is called for the first time;
# see MFsetupMixin._set_lake_bathymetry)
bathy = model.lake_bathymetry

# save a copy of original top elevations
Expand Down Expand Up @@ -1417,17 +1422,39 @@ def setup_array(model, package, var, data=None,
shutil.copy(model.cfg['intermediate_data']['top'][0],
original_top_file)
top = model.load_array(original_top_file)
lake_botm_elevations = top[bathy != 0] - bathy[bathy != 0]
if model.version == 'mf6':
# reset the model top to the lake bottom
top[bathy != 0] -= bathy[bathy != 0]

# update the top in the model
# todo: refactor this code to be consolidated with _set_idomain and setup_dis
model._setup_array('dis', 'top',
data={0: top},
datatype='array2d', resample_method='linear',
write_fmt='%.2f', dtype=float)
if hasattr(model, 'dis') and model.dis is not None:
if model.version == 'mf6':
model.dis.top = model.cfg['dis']['griddata']['top']
else:
model.dis.top = model.cfg['dis']['top'][0]

# write the top array again
top_filepath = model.setup_external_filepaths(package, 'top',
model.cfg[package]['top_filename_fmt'])[0]
save_array(top_filepath, top,
nodata=write_nodata,
fmt=write_fmt)
# if loading the model; use the model top that was just loaded in
else:
top_filename = model.cfg['dis']['griddata'].get('top')
if top_filename is None:
model.setup_external_filepaths('dis', 'top',
model.cfg['dis']['top_filename_fmt'])
top = model.load_array(model.cfg['dis']['griddata']['top'][0]['filename'])
if model.version == 'mf6':
top_filename = model.cfg['dis']['griddata']['top'][0]['filename']
else:
top_filename = model.cfg['dis']['top'][0]
top = model.load_array(top_filename)

# fill missing layers if any
if len(data) < model.nlay:
Expand Down Expand Up @@ -1492,19 +1519,8 @@ def setup_array(model, package, var, data=None,
raise Exception('Model layers less than {} {} thickness'.format(min_thickness,
model.length_units))
# fill any nan values that are above or below active cells to avoid cell thickness errors
top, botm = fill_cells_vertically(top, botm)
# the top may have been modified by fill_cells_vertically
# update the top in the model
# todo: refactor this code to be consolidated with _set_idomain and setup_dis
model._setup_array('dis', 'top',
data={0: top},
datatype='array2d', resample_method='linear',
write_fmt='%.2f', dtype=float)
if hasattr(model, 'dis') and model.dis is not None:
if model.version == 'mf6':
model.dis.top = model.cfg['dis']['griddata']['top']
else:
model.dis.top = model.cfg['dis']['top'][0]
botm = fill_cells_vertically(top, botm)

data = {i: arr for i, arr in enumerate(botm)}

# special case of LGR models
Expand Down Expand Up @@ -1603,19 +1619,6 @@ def setup_array(model, package, var, data=None,
dst = model.cfg['intermediate_data'][var][i]
shutil.copy(src, dst)

# write the top array again, because top was filled
# with botm array above
if var == 'botm':
top_filepath = model.setup_external_filepaths(package, 'top',
model.cfg[package]['top_filename_fmt'])[0]
save_array(top_filepath, top,
nodata=write_nodata,
fmt=write_fmt)
if model.version == 'mf6':
src = filepaths[i]['filename']
dst = model.cfg['intermediate_data'][var][i]
shutil.copy(src, dst)


def get_source_data_file_ext(cfg_data, package, var):
if 'filenames' in cfg_data:
Expand Down
10 changes: 7 additions & 3 deletions mfsetup/tests/test_discretization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,12 +118,16 @@ def test_fill_na(all_layers):
botm[-2] = 2
botm[-2, 2, 2] = np.nan

top, botm = fill_cells_vertically(top, botm)
botm = fill_cells_vertically(top, botm)
filled = all_layers.copy()
filled[0] = top
filled[1:] = botm
assert filled[:, 2, 2].tolist() == [10., 10., 8., 8., 8., 5., 5., 5., 5, 1.]
assert filled[:, 0, 0].tolist() == [8] * 8 + [2, 1]
# only compare layers 1:
# because fill_cells_vertically doesn't modify the model top
# (assuming that in most cases, there will be value
# for the model top elevation at each cell location)
assert filled[1:, 0, 0].tolist() == [8.] * 7 + [2., 1.]


def test_make_ibound(all_layers):
Expand All @@ -134,7 +138,7 @@ def test_make_ibound(all_layers):
botm[-2] = 2
botm[-2, 2, 2] = np.nan
botm[:, 3, 3] = np.arange(1, 10)[::-1] # column of cells all eq. to min thickness
filled_top, filled_botm = fill_cells_vertically(top, botm)
filled_botm = fill_cells_vertically(top, botm)
ibound = make_ibound(top, botm, nodata=nodata,
minimum_layer_thickness=1,
drop_thin_cells=True,
Expand Down
Loading