From a57b3a674948ea3239925258e167a259e7a1e55b Mon Sep 17 00:00:00 2001 From: apryet Date: Thu, 4 Apr 2024 14:19:30 -0400 Subject: [PATCH 1/4] Specifying raster file to zonal_stats in setup_basic_stress_data() lead to proj issues. Bypassed when source raster is loaded with rasterio and affine transform specified --- mfsetup/bcs.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/mfsetup/bcs.py b/mfsetup/bcs.py index 7869c2a8..ba2b1983 100644 --- a/mfsetup/bcs.py +++ b/mfsetup/bcs.py @@ -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]) From 1d02f714e9013a8d44d2346b8882fa76a3c17925 Mon Sep 17 00:00:00 2001 From: apryet Date: Tue, 25 Jun 2024 17:38:14 -0400 Subject: [PATCH 2/4] added simulate_high_k_lakes condition to setup_array in sourcedata --- mfsetup/sourcedata.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mfsetup/sourcedata.py b/mfsetup/sourcedata.py index 5ed64c79..acfa0672 100644 --- a/mfsetup/sourcedata.py +++ b/mfsetup/sourcedata.py @@ -1381,7 +1381,7 @@ def setup_array(model, package, var, data=None, # special handling of some variables # (for lakes) simulate_high_k_lakes = model.cfg['high_k_lakes']['simulate_high_k_lakes'] - if var == 'botm': + if var == 'botm' and simulate_high_k_lakes: # only execute this code if building the model (not loading) if not model._load: bathy = model.lake_bathymetry @@ -1605,7 +1605,9 @@ def setup_array(model, package, var, data=None, # write the top array again, because top was filled # with botm array above - if var == 'botm': + # NOTE: that tends to corrupt the top array, + # is that only applicable to simulate_high_k_lakes ?) + if var == 'botm' and simulate_high_k_lakes: top_filepath = model.setup_external_filepaths(package, 'top', model.cfg[package]['top_filename_fmt'])[0] save_array(top_filepath, top, From cc70d27abf43a5380bb0f047f9b7793875393fca Mon Sep 17 00:00:00 2001 From: "Leaf, Andrew T" Date: Thu, 25 Jul 2024 11:18:55 -0500 Subject: [PATCH 3/4] fix(sourcedata.py): when setting up the bottom array(s), only make edits to the model top if the lake package is being built, and a lake bathymetry raster was supplied (to be subtracted from the initial model top, which presumably represents a typical lake surface elevation, so that the model top is set at the lake bottom. --- mfsetup/sourcedata.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/mfsetup/sourcedata.py b/mfsetup/sourcedata.py index acfa0672..4d1fbd78 100644 --- a/mfsetup/sourcedata.py +++ b/mfsetup/sourcedata.py @@ -1381,9 +1381,16 @@ def setup_array(model, package, var, data=None, # special handling of some variables # (for lakes) simulate_high_k_lakes = model.cfg['high_k_lakes']['simulate_high_k_lakes'] - if var == 'botm' and simulate_high_k_lakes: - # only execute this code if building the model (not loading) - if not model._load: + apply_lake_bathymetry_to_model_top = False + if var == 'botm': + 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) + apply_lake_bathymetry_to_model_top = True bathy = model.lake_bathymetry # save a copy of original top elevations @@ -1417,7 +1424,6 @@ 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] @@ -1427,7 +1433,11 @@ def setup_array(model, package, var, data=None, 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: @@ -1605,9 +1615,9 @@ def setup_array(model, package, var, data=None, # write the top array again, because top was filled # with botm array above - # NOTE: that tends to corrupt the top array, + # NOTE: that tends to corrupt the top array, # is that only applicable to simulate_high_k_lakes ?) - if var == 'botm' and simulate_high_k_lakes: + if var == 'botm' and apply_lake_bathymetry_to_model_top: top_filepath = model.setup_external_filepaths(package, 'top', model.cfg[package]['top_filename_fmt'])[0] save_array(top_filepath, top, From 0be02d45ed58b12b878949f7f4defc52ebda87b3 Mon Sep 17 00:00:00 2001 From: "Leaf, Andrew T" Date: Mon, 23 Sep 2024 10:30:45 -0500 Subject: [PATCH 4/4] refactor(discretization.fill_cells_vertically): only return bottom array, since top shouldn't be modified; refactor(sourcedata.setup_array): move resetting of model top to branch where the top modified by supplied lake bathymetry (since this should be the only situation where this is needed) --- mfsetup/discretization.py | 6 ++-- mfsetup/sourcedata.py | 51 ++++++++++++---------------- mfsetup/tests/test_discretization.py | 10 ++++-- 3 files changed, 31 insertions(+), 36 deletions(-) diff --git a/mfsetup/discretization.py b/mfsetup/discretization.py index 7b95bede..666c5a55 100644 --- a/mfsetup/discretization.py +++ b/mfsetup/discretization.py @@ -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, @@ -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 diff --git a/mfsetup/sourcedata.py b/mfsetup/sourcedata.py index 4d1fbd78..6727bde2 100644 --- a/mfsetup/sourcedata.py +++ b/mfsetup/sourcedata.py @@ -1381,7 +1381,6 @@ def setup_array(model, package, var, data=None, # special handling of some variables # (for lakes) simulate_high_k_lakes = model.cfg['high_k_lakes']['simulate_high_k_lakes'] - apply_lake_bathymetry_to_model_top = False if var == 'botm': if 'lak' in model.cfg['model']['packages'] and not model._load and\ np.sum(model.lake_bathymetry) != 0: @@ -1390,7 +1389,6 @@ def setup_array(model, package, var, data=None, # (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) - apply_lake_bathymetry_to_model_top = True bathy = model.lake_bathymetry # save a copy of original top elevations @@ -1427,6 +1425,25 @@ def setup_array(model, package, var, data=None, 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') @@ -1502,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 @@ -1613,21 +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 - # NOTE: that tends to corrupt the top array, - # is that only applicable to simulate_high_k_lakes ?) - if var == 'botm' and apply_lake_bathymetry_to_model_top: - 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: diff --git a/mfsetup/tests/test_discretization.py b/mfsetup/tests/test_discretization.py index 73d55446..fb6fc359 100644 --- a/mfsetup/tests/test_discretization.py +++ b/mfsetup/tests/test_discretization.py @@ -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): @@ -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,