From acc8d9395afaff6271729acc8f479a6dab086fc5 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Mon, 17 Jun 2024 07:48:58 -0400 Subject: [PATCH] refactor: np.where(cond) -> np.asarray(cond).nonzero() The docs for `np.where()` (https://numpy.org/doc/stable/reference/generated/numpy.where.html) suggest to prefer `nonzero()` over `where()` without `x` and `y` arguments. In the spirit of defensive programming I included `np.asarray(cond)` even where `cond` is already an array. This PR also fixes a bug I introduced in the model splitter in #2124: while E711 (https://www.flake8rules.com/rules/E711.html) dictates comparisons to `None` should use identity rather than equality, this rule should not be applied to NumPy array selection conditions as it will change the semantics: ```shell >>> a = np.array([None, None]) >>> a[a != None] array([], dtype=object) >>> a[a is not None] array([[None, None]], dtype=object) ``` Unrelatedly, mark `test_mt3d.py::test_mfnwt_keat_uzf()` slow, it should not be included in smoke tests. --- autotest/conftest.py | 4 +- autotest/test_lake_connections.py | 6 +- autotest/test_mt3d.py | 1 + autotest/test_sfr.py | 2 +- autotest/test_zonbud_utility.py | 4 +- flopy/discretization/grid.py | 6 +- flopy/discretization/structuredgrid.py | 4 +- flopy/export/netcdf.py | 2 +- flopy/export/utils.py | 2 +- flopy/export/vtk.py | 6 +- flopy/mf6/utils/lakpak_utils.py | 6 +- flopy/mf6/utils/model_splitter.py | 94 +++++++------ flopy/modflow/mfrch.py | 4 +- flopy/modflow/mfsfr2.py | 6 +- flopy/pest/params.py | 2 +- flopy/plot/plotutil.py | 6 +- flopy/utils/binaryfile.py | 30 ++-- flopy/utils/check.py | 2 +- flopy/utils/compare.py | 4 +- flopy/utils/cvfdutil.py | 2 +- flopy/utils/datafile.py | 8 +- flopy/utils/flopy_io.py | 2 +- flopy/utils/formattedfile.py | 2 +- flopy/utils/geometry.py | 4 +- flopy/utils/gridgen.py | 2 +- flopy/utils/lgrutil.py | 2 +- flopy/utils/observationfile.py | 4 +- flopy/utils/particletrackfile.py | 14 +- flopy/utils/rasters.py | 4 +- flopy/utils/sfroutputfile.py | 4 +- flopy/utils/swroutputfile.py | 4 +- flopy/utils/triangle.py | 2 +- flopy/utils/util_list.py | 10 +- flopy/utils/voronoi.py | 4 +- flopy/utils/zonbud.py | 184 ++++++++++++------------- 35 files changed, 235 insertions(+), 208 deletions(-) diff --git a/autotest/conftest.py b/autotest/conftest.py index 6f770ef54d..13e25a530a 100644 --- a/autotest/conftest.py +++ b/autotest/conftest.py @@ -1,12 +1,10 @@ import re from importlib import metadata -from io import BytesIO, StringIO from pathlib import Path from platform import system -from typing import List, Optional +from typing import List import matplotlib.pyplot as plt -import numpy as np import pytest from modflow_devtools.misc import is_in_ci diff --git a/autotest/test_lake_connections.py b/autotest/test_lake_connections.py index ae0690b55d..71bab3d002 100644 --- a/autotest/test_lake_connections.py +++ b/autotest/test_lake_connections.py @@ -216,7 +216,7 @@ def test_lake(function_tmpdir, example_data_path): # mm.plot_array(bot_tm) # determine a reasonable lake bottom - idx = np.where(lakes > -1) + idx = np.asarray(lakes > -1).nonzero() lak_bot = bot_tm[idx].max() + 2.0 # interpolate top elevations @@ -634,9 +634,9 @@ def test_embedded_lak_prudic_mixed(example_data_path): lake_map[0, :, :] = lakibd[:, :] - 1 lakebed_leakance = np.zeros(shape2d, dtype=object) - idx = np.where(lake_map[0, :, :] == 0) + idx = np.asarray(lake_map[0, :, :] == 0).nonzero() lakebed_leakance[idx] = "none" - idx = np.where(lake_map[0, :, :] == 1) + idx = np.asarray(lake_map[0, :, :] == 1).nonzero() lakebed_leakance[idx] = 1.0 lakebed_leakance = lakebed_leakance.tolist() diff --git a/autotest/test_mt3d.py b/autotest/test_mt3d.py index 27dfafff9a..313809cc25 100644 --- a/autotest/test_mt3d.py +++ b/autotest/test_mt3d.py @@ -287,6 +287,7 @@ def test_mf2000_zeroth(function_tmpdir, mf2kmt3d_model_path): assert success, f"{mt.name} did not run" +@pytest.mark.slow @flaky(max_runs=3) @requires_exe("mfnwt", "mt3dms") @excludes_platform( diff --git a/autotest/test_sfr.py b/autotest/test_sfr.py index 452e871231..c50dcffe50 100644 --- a/autotest/test_sfr.py +++ b/autotest/test_sfr.py @@ -236,7 +236,7 @@ def interpolate_to_reaches(sfr): sfr.get_slopes(minimum_slope=-100, maximum_slope=100) reach_inds = 29 outreach = sfr.reach_data.outreach[reach_inds] - out_inds = np.where(sfr.reach_data.reachID == outreach) + out_inds = np.asarray(sfr.reach_data.reachID == outreach).nonzero() assert ( sfr.reach_data.slope[reach_inds] == ( diff --git a/autotest/test_zonbud_utility.py b/autotest/test_zonbud_utility.py index 8b4f35b66d..86991e6b26 100644 --- a/autotest/test_zonbud_utility.py +++ b/autotest/test_zonbud_utility.py @@ -113,8 +113,8 @@ def test_compare2zonebudget(cbc_f, zon_f, zbud_f, rtol): zb_arr = zba[zba["totim"] == time] fp_arr = fpa[fpa["totim"] == time] for name in fp_arr["name"]: - r1 = np.where(zb_arr["name"] == name) - r2 = np.where(fp_arr["name"] == name) + r1 = np.asarray(zb_arr["name"] == name).nonzero() + r2 = np.asarray(fp_arr["name"] == name).nonzero() if r1[0].shape[0] < 1 or r2[0].shape[0] < 1: continue if r1[0].shape[0] != r2[0].shape[0]: diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 0a091320de..2681a671cd 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -455,15 +455,15 @@ def saturated_thickness(self, array, mask=None): bot = self.remove_confining_beds(bot) array = self.remove_confining_beds(array) - idx = np.where((array < top) & (array > bot)) + idx = np.asarray((array < top) & (array > bot)).nonzero() thickness[idx] = array[idx] - bot[idx] - idx = np.where(array <= bot) + idx = np.asarray(array <= bot).nonzero() thickness[idx] = 0.0 if mask is not None: if isinstance(mask, (float, int)): mask = [float(mask)] for mask_value in mask: - thickness[np.where(array == mask_value)] = np.nan + thickness[np.asarray(array == mask_value).nonzero()] = np.nan return thickness def saturated_thick(self, array, mask=None): diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index ba5a1a4302..13f15346ce 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -930,7 +930,7 @@ def intersect(self, x, y, z=None, local=False, forgive=False): "x, y point given is outside of the model area" ) else: - col = np.where(xcomp)[0][-1] + col = np.asarray(xcomp).nonzero()[0][-1] ycomp = y < ye if np.all(ycomp) or not np.any(ycomp): @@ -941,7 +941,7 @@ def intersect(self, x, y, z=None, local=False, forgive=False): "x, y point given is outside of the model area" ) else: - row = np.where(ycomp)[0][-1] + row = np.asarray(ycomp).nonzero()[0][-1] if np.any(np.isnan([row, col])): row = col = np.nan if z is not None: diff --git a/flopy/export/netcdf.py b/flopy/export/netcdf.py index cfc2907d23..eda44e1172 100644 --- a/flopy/export/netcdf.py +++ b/flopy/export/netcdf.py @@ -632,7 +632,7 @@ def difference( d_data[np.isnan(d_data)] = FILLVALUE if mask_zero_diff: - d_data[np.where(d_data == 0.0)] = FILLVALUE + d_data[np.asarray(d_data == 0.0).nonzero()] = FILLVALUE var = new_net.create_variable( vname, attrs, s_var.dtype, dimensions=s_var.dimensions diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 5754134360..d0b5773983 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -203,7 +203,7 @@ def _add_output_nc_variable( logger.log(f"creating array for {var_name}") for mask_val in mask_vals: - array[np.where(array == mask_val)] = np.nan + array[np.asarray(array == mask_val).nonzero()] = np.nan mx, mn = np.nanmax(array), np.nanmin(array) array[np.isnan(array)] = netcdf.FILLVALUE diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 16aa5dcc0d..e25c94695e 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -574,7 +574,9 @@ def _build_hfbs(self, pkg): pts = [] for v in v1: - ix = np.where((v2.T[0] == v[0]) & (v2.T[1] == v[1])) + ix = np.asarray( + (v2.T[0] == v[0]) & (v2.T[1] == v[1]) + ).nonzero() if len(ix[0]) > 0 and len(pts) < 2: pts.append(v2[ix[0][0]]) @@ -652,7 +654,7 @@ def _build_point_scalar_array(self, array): ps_array[pt] = array[value["idx"][ix]] else: ps_graph = self._point_scalar_numpy_graph.copy() - idxs = np.where(np.isnan(array)) + idxs = np.asarray(np.isnan(array)).nonzero() not_graphed = np.isin(ps_graph, idxs[0]) ps_graph[not_graphed] = -1 ps_array = np.where(ps_graph >= 0, array[ps_graph], np.nan) diff --git a/flopy/mf6/utils/lakpak_utils.py b/flopy/mf6/utils/lakpak_utils.py index 9dc3282936..8b1e1b25b1 100644 --- a/flopy/mf6/utils/lakpak_utils.py +++ b/flopy/mf6/utils/lakpak_utils.py @@ -125,7 +125,7 @@ def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=None): unique = np.unique(lake_map) # exclude lakes with lake numbers less than 0 - idx = np.where(unique > -1) + idx = np.asarray(unique > -1).nonzero() unique = unique[idx] dx, dy = None, None @@ -199,7 +199,9 @@ def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=None): # reset idomain for lake if iconn > 0: - idx = np.where((lake_map == lake_number) & (idomain > 0)) + idx = np.asarray( + (lake_map == lake_number) & (idomain > 0) + ).nonzero() idomain[idx] = 0 return idomain, connection_dict, connectiondata diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index 7313877618..b16ba93b6b 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -316,7 +316,7 @@ def load_node_mapping(self, sim, filename): for mkey in models: ncpl = self._new_ncpl[mkey] array = np.full((ncpl,), -1, dtype=int) - onode = np.where(model_array == mkey)[0] + onode = np.asarray(model_array == mkey).nonzero()[0] nnode = split_array[onode] array[nnode] = onode grid_info[mkey] = (array,) @@ -413,7 +413,7 @@ def optimize_splitting_mask(self, nparts): membership = np.array(membership, dtype=int) if laks: for lak in laks: - idx = np.where(lak_array == lak)[0] + idx = np.asarray(lak_array == lak).nonzero()[0] mnum = np.unique(membership[idx])[0] membership[idx] = mnum @@ -429,7 +429,7 @@ def optimize_splitting_mask(self, nparts): ev = np.equal(mnums1, mnums2) if np.all(ev): continue - idx = np.where(~ev)[0] + idx = np.asarray(~ev).nonzero()[0] mnum_to = mnums1[idx] adj_nodes = nodes2[idx] membership[adj_nodes] = mnum_to @@ -471,7 +471,7 @@ def reconstruct_array(self, arrays): array = array.ravel() ncpl = self._new_ncpl[mkey] mapping = self._grid_info[mkey][-1] - old_nodes = np.where(mapping != -1) + old_nodes = np.asarray(mapping != -1).nonzero() new_nodes = mapping[old_nodes] old_nodes = np.tile(old_nodes, (nlay, 1)) @@ -645,7 +645,7 @@ def _remap_nodes(self, array): bad_keys = [] for mkey in mkeys: count = 0 - mask = np.where(array == mkey) + mask = np.asarray(array == mkey).nonzero() for arr in idomain: check = arr[mask] count += np.count_nonzero(check) @@ -670,7 +670,7 @@ def _remap_nodes(self, array): if self._modelgrid.grid_type == "structured": a = array.reshape(self._modelgrid.nrow, self._modelgrid.ncol) for m in np.unique(a): - cells = np.where(a == m) + cells = np.asarray(a == m).nonzero() rmin, rmax = np.min(cells[0]), np.max(cells[0]) cmin, cmax = np.min(cells[1]), np.max(cells[1]) cellids = list(zip([0] * len(cells[0]), cells[0], cells[1])) @@ -702,7 +702,7 @@ def _remap_nodes(self, array): xverts, yverts = None, None for m in np.unique(array): - cells = np.where(array == m)[0] + cells = np.asarray(array == m).nonzero()[0] mapping = np.zeros( ( len( @@ -718,9 +718,9 @@ def _remap_nodes(self, array): if xverts is not None: mxv = xverts[cells] myv = yverts[cells] - xmidx = np.where(mxv == np.nanmin(mxv))[0] + xmidx = np.asarray(mxv == np.nanmin(mxv)).nonzero()[0] myv = myv[xmidx] - ymidx = np.where(myv == np.nanmin(myv))[0] + ymidx = np.asarray(myv == np.nanmin(myv)).nonzero()[0] self._offsets[m] = { "xorigin": np.nanmin(mxv[xmidx[0]]), @@ -736,11 +736,11 @@ def _remap_nodes(self, array): new_ncpl[m] *= i for mdl in np.unique(array): - mnodes = np.where(array == mdl)[0] + mnodes = np.asarray(array == mdl).nonzero()[0] mg_info = grid_info[mdl] if mg_info is not None: mapping = mg_info[-1] - new_nodes = np.where(mapping != -1)[0] + new_nodes = np.asarray(mapping != -1).nonzero()[0] old_nodes = mapping[new_nodes] for ix, nnode in enumerate(new_nodes): self._node_map[old_nodes[ix]] = (mdl, nnode) @@ -1163,7 +1163,7 @@ def _remap_array(self, item, mfarray, mapped_data, **kwargs): new_ncpl = self._new_ncpl[mkey] new_array = np.zeros(new_ncpl * nlay, dtype=dtype) mapping = self._grid_info[mkey][-1] - new_nodes = np.where(mapping != -1) + new_nodes = np.asarray(mapping != -1).nonzero() old_nodes = mapping[new_nodes] old_nodes = np.tile(old_nodes, (nlay, 1)) @@ -1263,7 +1263,7 @@ def _remap_mflist( new_model, new_node = self._get_new_model_new_node(nodes) for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if self._pkg_mover and transient: mvr_remap = { idx[i]: (model.name, i) for i in range(len(idx)) @@ -1363,7 +1363,7 @@ def _remap_uzf(self, package, mapped_data): name = package.filename self._uzf_remaps[name] = {} for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if len(idx) == 0: new_recarray = None else: @@ -1401,7 +1401,9 @@ def _remap_uzf(self, package, mapped_data): spd = {} for per, recarray in perioddata.items(): - idx = np.where(np.isin(recarray.ifno, uzf_nodes)) + idx = np.asarray( + np.isin(recarray.ifno, uzf_nodes) + ).nonzero() new_period = recarray[idx] new_period["ifno"] = [ uzf_remap[i] for i in new_period["ifno"] @@ -1547,7 +1549,7 @@ def _remap_lak(self, package, mapped_data): new_model, new_node = self._get_new_model_new_node(nodes) for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if len(idx) == 0: new_recarray = None else: @@ -1586,7 +1588,9 @@ def _remap_lak(self, package, mapped_data): if meta[0] == mkey: mapnos.append(lak) - idxs = np.where(np.isin(outlets.lakein, mapnos))[0] + idxs = np.asarray( + np.isin(outlets.lakein, mapnos) + ).nonzero()[0] if len(idxs) == 0: new_outlets = None else: @@ -1680,7 +1684,7 @@ def _remap_sfr(self, package, mapped_data): new_model, new_node = self._get_new_model_new_node(nodes) for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if len(idx) == 0: new_recarray = None continue @@ -1709,7 +1713,9 @@ def _remap_sfr(self, package, mapped_data): ) # now let's remap connection data and tag external exchanges - idx = np.where(np.isin(connectiondata.ifno, old_rno))[0] + idx = np.asarray( + np.isin(connectiondata.ifno, old_rno) + ).nonzero()[0] new_connectiondata = connectiondata[idx] ncons = [] for ix, rec in enumerate(new_connectiondata): @@ -1776,8 +1782,12 @@ def _remap_sfr(self, package, mapped_data): if m0 != m1: div_mover_ix.append(ix) - idx = np.where(np.isin(diversions.ifno, old_rno))[0] - idx = np.where(~np.isin(idx, div_mover_ix))[0] + idx = np.asarray( + np.isin(diversions.ifno, old_rno) + ).nonzero()[0] + idx = np.asarray( + ~np.isin(idx, div_mover_ix) + ).nonzero()[0] new_diversions = diversions[idx] new_rno = [ @@ -1802,23 +1812,25 @@ def _remap_sfr(self, package, mapped_data): # now we can do the stress period data spd = {} for kper, recarray in perioddata.items(): - idx = np.where(np.isin(recarray.ifno, old_rno))[0] + idx = np.asarray( + np.isin(recarray.ifno, old_rno) + ).nonzero()[0] new_spd = recarray[idx] if diversions is not None: - external_divs = np.where( + external_divs = np.asarray( np.isin(new_spd.idv, list(div_mvr_conn.keys())) - )[0] + ).nonzero()[0] if len(external_divs) > 0: for ix in external_divs: rec = recarray[ix] idv = recarray["idv"] div_mvr_conn[idv].append(rec["divflow"]) - idx = np.where( + idx = np.asarray( ~np.isin( new_spd.idv, list(div_mvr_conn.keys()) ) - )[0] + ).nonzero()[0] new_spd = new_spd[idx] @@ -1931,7 +1943,7 @@ def _remap_maw(self, package, mapped_data): maw_remaps = {} for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] new_connectiondata = connectiondata[idx] if len(new_connectiondata) == 0: continue @@ -1965,7 +1977,9 @@ def _remap_maw(self, package, mapped_data): spd = {} for per, recarray in perioddata.items(): - idx = np.where(np.isin(recarray.ifno, maw_wellnos))[0] + idx = np.asarray( + np.isin(recarray.ifno, maw_wellnos) + ).nonzero()[0] if len(idx) > 0: new_recarray = recarray[idx] new_wellno = [ @@ -2030,7 +2044,7 @@ def _remap_csub(self, package, mapped_data): ninterbeds = None for mkey, model in self._model_dict.items(): - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if len(idx) == 0: new_packagedata = None else: @@ -2052,7 +2066,7 @@ def _remap_csub(self, package, mapped_data): layers, nodes = self._cellid_to_layer_node(recarray.cellid) new_model, new_node = self._get_new_model_new_node(nodes) - idx = np.where(new_model == mkey)[0] + idx = np.asarray(new_model == mkey).nonzero()[0] if len(idx) == 0: continue @@ -2158,7 +2172,7 @@ def _remap_hfb(self, package, mapped_data): raise AssertionError("Models cannot be split along faults") for mkey, model in self._model_dict.items(): - idx = np.where(new_model1 == mkey)[0] + idx = np.asarray(new_model1 == mkey).nonzero()[0] if len(idx) == 0: new_recarray = None else: @@ -2262,7 +2276,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): dtype=object, ) for mkey, model in self._model_dict.items(): - idx = np.where(new_model1 == mkey) + idx = np.asarray(new_model1 == mkey).nonzero() tmp_cellid = self._new_node_to_cellid( model, new_node1, layers1, idx ) @@ -2297,7 +2311,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): ) for idt in set(idtype): remaps = remapper[idt] - idx = np.where(idtype == idt) + idx = np.asarray(idtype == idt).nonzero() new_cellid1[idx] = [ ( remaps[i][-1] + 1 @@ -2364,7 +2378,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): dtype=object, ) for mkey, model in self._model_dict.items(): - idx = np.where(new_model1 == mkey) + idx = np.asarray(new_model1 == mkey).nonzero() idx = [ ix for ix, i in enumerate(recarray.id[idx]) @@ -2399,7 +2413,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): new_model1[mm_idx] = tmp_models cellid2 = recarray.id2 - conv_idx = np.where((cellid2 is not None))[0] + conv_idx = np.asarray(cellid2 != None).nonzero()[0] # noqa: E711 if len(conv_idx) > 0: # do stuff # need to trap layers... if pkg_type is None: @@ -2454,9 +2468,9 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): (len(new_node2),), None, dtype=object ) for mkey, model in self._model_dict.items(): - idx = np.where(new_model2 == mkey) + idx = np.asarray(new_model2 == mkey).nonzero() tmp_node = new_node2[idx] - cidx = np.where((tmp_node is not None)) + cidx = np.asarray((tmp_node != None)).nonzero() # noqa: E711 tmp_cellid = model.modelgrid.get_lrc( tmp_node[cidx].to_list() ) @@ -2501,7 +2515,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): if idt is None: continue remaps = remapper[idt] - idx = np.where(idtype == idt) + idx = np.asarray(idtype == idt).nonzero() new_cellid2[idx] = [ ( remaps[i][-1] + 1 @@ -2536,7 +2550,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): new_model1[idx] = mkey # now we remap the continuous data!!!! - idx = np.where(new_model1 == mkey)[0] + idx = np.asarray(new_model1 == mkey).nonzero()[0] if len(idx) == 0: continue @@ -2682,7 +2696,7 @@ def _remap_adv_tag(self, mkey, recarray, item, mapper): if meta[0] == mkey: mapnos.append(lak) - idxs = np.where(np.isin(recarray[item], mapnos))[0] + idxs = np.asarray(np.isin(recarray[item], mapnos)).nonzero()[0] if len(idxs) == 0: new_recarray = None else: diff --git a/flopy/modflow/mfrch.py b/flopy/modflow/mfrch.py index 4803dcfdc4..a9a22e394f 100644 --- a/flopy/modflow/mfrch.py +++ b/flopy/modflow/mfrch.py @@ -235,8 +235,8 @@ def check( if Tmean != 0: R_T = period_means / Tmean - lessthan = np.where(R_T < RTmin)[0] - greaterthan = np.where(R_T > RTmax)[0] + lessthan = np.asarray(R_T < RTmin).nonzero()[0] + greaterthan = np.asarray(R_T > RTmax).nonzero()[0] if len(lessthan) > 0: txt = ( diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index e21915e29e..62112fb51a 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -1511,7 +1511,7 @@ def plot_path(self, start_seg=None, end_seg=0, plot_segment_lines=True): # slice the path path = np.array(self.paths[start_seg]) - endidx = np.where(path == end_seg)[0] + endidx = np.asarray(path == end_seg).nonzero()[0] endidx = endidx if len(endidx) > 0 else None path = path[: np.squeeze(endidx)] path = [s for s in path if s > 0] # skip lakes for now @@ -1523,7 +1523,7 @@ def plot_path(self, start_seg=None, end_seg=0, plot_segment_lines=True): dist = np.cumsum(tmp.rchlen.values) * to_miles.get(mfunits, 1.0) # segment starts - starts = dist[np.where(tmp.ireach.values == 1)[0]] + starts = dist[np.asarray(tmp.ireach.values == 1).nonzero()[0]] ax = plt.subplots(figsize=(11, 8.5))[-1] ax.plot(dist, tops, label="Model top") @@ -2411,7 +2411,7 @@ def routing(self): # max node with * a tolerance # 1.25 * hyp is greater than distance of two diagonally adjacent nodes # where one is 1.5x larger than the other - breaks = np.where(dist > hyp * 1.25) + breaks = np.asarray(dist > hyp * 1.25).nonzero() breaks_reach_data = rd[breaks] segments_with_breaks = set(breaks_reach_data.iseg) if len(breaks) > 0: diff --git a/flopy/pest/params.py b/flopy/pest/params.py index 7dc2535f9e..ece259d54f 100644 --- a/flopy/pest/params.py +++ b/flopy/pest/params.py @@ -74,7 +74,7 @@ def zonearray2params( plist = [] for i, iz in enumerate(parzones): span = {} - span["idx"] = np.where(zonearray == iz) + span["idx"] = np.asarray(zonearray == iz).nonzero() parname = f"{partype}_{iz}" startvalue = parvals[i] p = Params( diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 2a55346d10..4efab156a6 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -1647,7 +1647,7 @@ def line_intersect_grid(ptsin, xgrid, ygrid): numb = (x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3) denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1) ua = np.ones(denom.shape, dtype=denom.dtype) * np.nan - idx = np.where(denom != 0.0) + idx = np.asarray(denom != 0.0).nonzero() ua[idx] = numa[idx] / denom[idx] del numa del numb @@ -2231,7 +2231,7 @@ def advanced_package_bc_helper(pkg, modelgrid, kper): idx = np.array([list(i) for i in mflist["cellid"]], dtype=int).T else: iuzfbnd = pkg.iuzfbnd.array - idx = np.where(iuzfbnd != 0) + idx = np.asarray(iuzfbnd != 0).nonzero() idx = np.append([[0] * idx[-1].size], idx, axis=0) elif pkg.package_type in ("lak", "maw"): if pkg.parent.version == "mf6": @@ -2239,7 +2239,7 @@ def advanced_package_bc_helper(pkg, modelgrid, kper): idx = np.array([list(i) for i in mflist["cellid"]], dtype=int).T else: lakarr = pkg.lakarr.array[kper] - idx = np.where(lakarr != 0) + idx = np.asarray(lakarr != 0).nonzero() idx = np.array(idx) else: raise NotImplementedError( diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index 157cb17404..94074f4275 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -598,7 +598,7 @@ def get_ts(self, idx): # Find the time index and then put value into result in the # correct location. - itim = np.where(result[:, 0] == header["totim"])[0] + itim = np.asarray(result[:, 0] == header["totim"]).nonzero()[0] result[itim, istat] = binaryread(self.file, self.realtype) istat += 1 return result @@ -887,7 +887,9 @@ def _get_data_array(self, totim=0.0): """ if totim >= 0.0: - keyindices = np.where(self.recordarray["totim"] == totim)[0] + keyindices = np.asarray( + self.recordarray["totim"] == totim + ).nonzero()[0] if len(keyindices) == 0: msg = f"totim value ({totim}) not found in file..." raise Exception(msg) @@ -1604,7 +1606,9 @@ def get_indices(self, text=None): # check and make sure that text is in file if text is not None: text16 = self._find_text(text) - select_indices = np.where(self.recordarray["text"] == text16) + select_indices = np.asarray( + self.recordarray["text"] == text16 + ).nonzero() if isinstance(select_indices, tuple): select_indices = select_indices[0] else: @@ -1869,7 +1873,7 @@ def get_ts(self, idx, text=None, times=None): for vv in v: field = vv.dtype.names[2] - dix = np.where(np.isin(vv["node"], ndx))[0] + dix = np.asarray(np.isin(vv["node"], ndx)).nonzero()[0] if len(dix) > 0: result[itim, 1:] = vv[field][dix] @@ -2176,7 +2180,9 @@ def get_residual(self, totim, scaled=False): residual = np.zeros((nlay, nrow, ncol), dtype=float) if scaled: inflow = np.zeros((nlay, nrow, ncol), dtype=float) - select_indices = np.where(self.recordarray["totim"] == totim)[0] + select_indices = np.asarray( + self.recordarray["totim"] == totim + ).nonzero()[0] for i in select_indices: text = self.recordarray[i]["text"].decode() @@ -2187,9 +2193,9 @@ def get_residual(self, totim, scaled=False): residual -= flow[:, :, :] residual[:, :, 1:] += flow[:, :, :-1] if scaled: - idx = np.where(flow < 0.0) + idx = np.asarray(flow < 0.0).nonzero() inflow[idx] -= flow[idx] - idx = np.where(flow > 0.0) + idx = np.asarray(flow > 0.0).nonzero() l, r, c = idx idx = (l, r, c + 1) inflow[idx] += flow[idx] @@ -2197,9 +2203,9 @@ def get_residual(self, totim, scaled=False): residual -= flow[:, :, :] residual[:, 1:, :] += flow[:, :-1, :] if scaled: - idx = np.where(flow < 0.0) + idx = np.asarray(flow < 0.0).nonzero() inflow[idx] -= flow[idx] - idx = np.where(flow > 0.0) + idx = np.asarray(flow > 0.0).nonzero() l, r, c = idx idx = (l, r + 1, c) inflow[idx] += flow[idx] @@ -2207,16 +2213,16 @@ def get_residual(self, totim, scaled=False): residual -= flow[:, :, :] residual[1:, :, :] += flow[:-1, :, :] if scaled: - idx = np.where(flow < 0.0) + idx = np.asarray(flow < 0.0).nonzero() inflow[idx] -= flow[idx] - idx = np.where(flow > 0.0) + idx = np.asarray(flow > 0.0).nonzero() l, r, c = idx idx = (l + 1, r, c) inflow[idx] += flow[idx] else: residual += flow if scaled: - idx = np.where(flow > 0.0) + idx = np.asarray(flow > 0.0).nonzero() inflow[idx] += flow[idx] if scaled: diff --git a/flopy/utils/check.py b/flopy/utils/check.py index 662643d60b..432dc52214 100644 --- a/flopy/utils/check.py +++ b/flopy/utils/check.py @@ -507,7 +507,7 @@ def values(self, a, criteria, error_name="", error_type="Warning"): True value in criteria. """ if np.any(criteria): - inds = np.where(criteria) + inds = np.asarray(criteria).nonzero() v = a[inds] # works with structured or unstructured pn = [self.package.name] * len(v) en = [error_name] * len(v) diff --git a/flopy/utils/compare.py b/flopy/utils/compare.py index 3e4100fba8..2e5f3fd6b9 100644 --- a/flopy/utils/compare.py +++ b/flopy/utils/compare.py @@ -40,7 +40,7 @@ def _diffmax(v1, v2): diff = abs(v1 - v2) diffmax = diff.max() - return diffmax, np.where(diff == diffmax) + return diffmax, np.asarray(diff == diffmax).nonzero() def _difftol(v1, v2, tol): @@ -75,7 +75,7 @@ def _difftol(v1, v2, tol): raise Exception(err) diff = abs(v1 - v2) - return diff.max(), np.where(diff > tol) + return diff.max(), np.asarray(diff > tol).nonzero() def compare_budget( diff --git a/flopy/utils/cvfdutil.py b/flopy/utils/cvfdutil.py index 3eb33d6e88..fea3b56027 100644 --- a/flopy/utils/cvfdutil.py +++ b/flopy/utils/cvfdutil.py @@ -324,7 +324,7 @@ def gridlist_to_verts(gridlist): vertdict = {} icell = 0 for sg in gridlist: - ilays, irows, icols = np.where(sg.idomain > 0) + ilays, irows, icols = np.asarray(sg.idomain > 0).nonzero() for _, i, j in zip(ilays, irows, icols): v = sg.get_cell_vertices(i, j) vertdict[icell] = v + [v[0]] diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 51244a4b42..ca922c516e 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -463,7 +463,9 @@ def _get_data_array(self, totim=0): """ if totim >= 0.0: - keyindices = np.where(self.recordarray["totim"] == totim)[0] + keyindices = np.asarray( + self.recordarray["totim"] == totim + ).nonzero()[0] if len(keyindices) == 0: msg = f"totim value ({totim}) not found in file..." raise Exception(msg) @@ -545,10 +547,10 @@ def get_data(self, kstpkper=None, idx=None, totim=None, mflay=None): if kstpkper is not None: kstp1 = kstpkper[0] + 1 kper1 = kstpkper[1] + 1 - idx = np.where( + idx = np.asarray( (self.recordarray["kstp"] == kstp1) & (self.recordarray["kper"] == kper1) - ) + ).nonzero() if idx[0].shape[0] == 0: raise Exception( f"get_data() error: kstpkper not found:{kstpkper}" diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py index c3f94b62f0..403f3e7a97 100644 --- a/flopy/utils/flopy_io.py +++ b/flopy/utils/flopy_io.py @@ -305,7 +305,7 @@ def flux_to_wel(cbc_file, text, precision="single", model=None, verbose=False): arr = arr[0] print(arr.max(), arr.min(), arr.sum()) # masked where zero - arr[np.where(arr == 0.0)] = np.nan + arr[np.asarray(arr == 0.0).nonzero()] = np.nan m4d[iper + 1] = arr iper += 1 diff --git a/flopy/utils/formattedfile.py b/flopy/utils/formattedfile.py index 02903ca219..ca9436ed21 100644 --- a/flopy/utils/formattedfile.py +++ b/flopy/utils/formattedfile.py @@ -303,7 +303,7 @@ def get_ts(self, idx): # Find the time index and then put value into result in the # correct location. - itim = np.where(result[:, 0] == header["totim"])[0] + itim = np.asarray(result[:, 0] == header["totim"]).nonzero()[0] result[itim, istat] = self._read_val(j) istat += 1 return result diff --git a/flopy/utils/geometry.py b/flopy/utils/geometry.py index dac27c1441..578040d559 100644 --- a/flopy/utils/geometry.py +++ b/flopy/utils/geometry.py @@ -872,9 +872,9 @@ def point_in_polygon(xc, yc, polygon): yc - polygon[i][1] ) / (polygon[j][1] - polygon[i][1]) - comp = np.where( + comp = np.asarray( ((polygon[i][1] > yc) ^ (polygon[j][1] > yc)) & (xc < tmp) - ) + ).nonzero() j = i if len(comp[0]) > 0: diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py index 27e9b1b344..437f3a0be5 100644 --- a/flopy/utils/gridgen.py +++ b/flopy/utils/gridgen.py @@ -733,7 +733,7 @@ def plot( shapename = os.path.join(self.model_ws, "qtgrid") xmin, xmax, ymin, ymax = shapefile_extents(shapename) - idx = np.where(self.qtra.layer == layer)[0] + idx = np.asarray(self.qtra.layer == layer).nonzero()[0] pc = plot_shapefile( shapename, diff --git a/flopy/utils/lgrutil.py b/flopy/utils/lgrutil.py index b17839483c..043b1a9990 100644 --- a/flopy/utils/lgrutil.py +++ b/flopy/utils/lgrutil.py @@ -161,7 +161,7 @@ def __init__( # idomain assert idomainp.shape == (nlayp, nrowp, ncolp) self.idomain = idomainp - idxl, idxr, idxc = np.where(idomainp == 0) + idxl, idxr, idxc = np.asarray(idomainp == 0).nonzero() assert idxl.shape[0] > 1, "no zero values found in idomain" # child cells per parent and child cells per parent layer diff --git a/flopy/utils/observationfile.py b/flopy/utils/observationfile.py index 8d85aaf3d7..3e5308b3aa 100644 --- a/flopy/utils/observationfile.py +++ b/flopy/utils/observationfile.py @@ -104,7 +104,7 @@ def get_data(self, idx=None, obsname=None, totim=None): i0 = 0 i1 = self.data.shape[0] if totim is not None: - idx = np.where(self.data["totim"] == totim)[0][0] + idx = np.asarray(self.data["totim"] == totim).nonzero()[0][0] i0 = idx i1 = idx + 1 elif idx is not None: @@ -183,7 +183,7 @@ def get_dataframe( i0 = 0 i1 = self.data.shape[0] if totim is not None: - idx = np.where(self.data["totim"] == totim)[0][0] + idx = np.asarray(self.data["totim"] == totim).nonzero()[0][0] i0 = idx i1 = idx + 1 elif idx is not None: diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py index 75ff3972da..50c004db8e 100644 --- a/flopy/utils/particletrackfile.py +++ b/flopy/utils/particletrackfile.py @@ -106,16 +106,16 @@ def get_data( """ data = self._data[list(self.outdtype.names)] if minimal else self._data idx = ( - np.where(data["particleid"] == partid)[0] + np.asarray(data["particleid"] == partid).nonzero()[0] if totim is None else ( - np.where( + np.asarray( (data["time"] >= totim) & (data["particleid"] == partid) - )[0] + ).nonzero()[0] if ge - else np.where( + else np.asarray( (data["time"] <= totim) & (data["particleid"] == partid) - )[0] + ).nonzero()[0] ) ) @@ -145,9 +145,9 @@ def get_alldata(self, totim=None, ge=True, minimal=False): data = self._data[list(self.outdtype.names)] if minimal else self._data if totim is not None: idx = ( - np.where(data["time"] >= totim)[0] + np.asarray(data["time"] >= totim).nonzero()[0] if ge - else np.where(data["time"] <= totim)[0] + else np.asarray(data["time"] <= totim).nonzero()[0] ) if len(idx) > 0: data = data[idx] diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index 72e32eabff..f8f965edb8 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -256,7 +256,7 @@ def sample_point(self, *point, band=1): dist = np.sqrt(xt + yt) # 3: find indices of minimum distance - md = np.where(dist == np.nanmin(dist)) + md = np.asarray(dist == np.nanmin(dist)).nonzero() # 4: sample the array and average if necessary vals = [] @@ -534,7 +534,7 @@ def crop(self, polygon, invert=False): xt = (pt[0] - xc) ** 2 yt = (pt[1] - yc) ** 2 hypot = np.sqrt(xt + yt) - ind = np.where(hypot == np.min(hypot)) + ind = np.asarray(hypot == np.min(hypot)).nonzero() yind.append(ind[0][0]) xind.append(ind[1][0]) diff --git a/flopy/utils/sfroutputfile.py b/flopy/utils/sfroutputfile.py index 7754d4da8f..ba23a9b5b4 100644 --- a/flopy/utils/sfroutputfile.py +++ b/flopy/utils/sfroutputfile.py @@ -147,7 +147,9 @@ def get_nstrm(df): Number of SFR cells """ - wherereach1 = np.where((df.segment == 1) & (df.reach == 1))[0] + wherereach1 = np.asarray( + (df.segment == 1) & (df.reach == 1) + ).nonzero()[0] if len(wherereach1) == 1: return len(df) elif len(wherereach1) > 1: diff --git a/flopy/utils/swroutputfile.py b/flopy/utils/swroutputfile.py index 991571ebab..eff2488688 100644 --- a/flopy/utils/swroutputfile.py +++ b/flopy/utils/swroutputfile.py @@ -231,11 +231,11 @@ def get_data(self, idx=None, kswrkstpkper=None, totim=None): kper1 = kswrkstpkper[2] totim1 = self._recordarray[ - np.where( + np.asarray( (self._recordarray["kswr"] == kswr1) & (self._recordarray["kstp"] == kstp1) & (self._recordarray["kper"] == kper1) - ) + ).nonzero() ]["totim"][0] elif totim is not None: totim1 = totim diff --git a/flopy/utils/triangle.py b/flopy/utils/triangle.py index a0d86d18f0..d001bd73ca 100644 --- a/flopy/utils/triangle.py +++ b/flopy/utils/triangle.py @@ -309,7 +309,7 @@ def plot_boundary(self, ibm, ax=None, **kwargs): """ if ax is None: ax = plt.gca() - idx = np.where(self.edge["boundary_marker"] == ibm)[0] + idx = np.asarray(self.edge["boundary_marker"] == ibm).nonzero()[0] for i in idx: iv1 = self.edge["endpoint1"][i] iv2 = self.edge["endpoint2"][i] diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index 382dbc65cf..8616cf11df 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -810,15 +810,15 @@ def check_kij(self): data = self[kper] if data is not None: k = data["k"] - k_idx = np.where(np.logical_or(k < 0, k >= nl)) + k_idx = np.asarray(np.logical_or(k < 0, k >= nl)).nonzero() if k_idx[0].shape[0] > 0: out_idx.extend(list(k_idx[0])) i = data["i"] - i_idx = np.where(np.logical_or(i < 0, i >= nr)) + i_idx = np.asarray(np.logical_or(i < 0, i >= nr)).nonzero() if i_idx[0].shape[0] > 0: out_idx.extend(list(i_idx[0])) j = data["j"] - j_idx = np.where(np.logical_or(j < 0, j >= nc)) + j_idx = np.asarray(np.logical_or(j < 0, j >= nc)).nonzero() if j_idx[0].shape[0]: out_idx.extend(list(j_idx[0])) @@ -887,7 +887,9 @@ def attribute_by_kper(self, attr, function=np.mean, idx_val=None): kper_data = self.__data[kper] if idx_val is not None: kper_data = kper_data[ - np.where(kper_data[idx_val[0]] == idx_val[1]) + np.asarray( + kper_data[idx_val[0]] == idx_val[1] + ).nonzero() ] v = function(kper_data[attr]) values.append(v) diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index 15982d205a..d3f52566ae 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -153,11 +153,11 @@ def tri2vor(tri, **kwargs): polygon = [(x, y) for x, y in tri._polygons[ipolygon]] vor_vert_notindomain = point_in_polygon(xc, yc, polygon) vor_vert_notindomain = vor_vert_notindomain.flatten() - idx = np.where(vor_vert_notindomain == True) + idx = np.asarray(vor_vert_notindomain == True).nonzero() vor_vert_indomain[idx] = False idx_vertindex = -1 * np.ones((nvertices), int) - idx_filtered = np.where(vor_vert_indomain == True) + idx_filtered = np.asarray(vor_vert_indomain == True).nonzero() nvalid_vertices = len(idx_filtered[0]) # renumber valid vertices consecutively idx_vertindex[idx_filtered] = np.arange(nvalid_vertices) diff --git a/flopy/utils/zonbud.py b/flopy/utils/zonbud.py index 37593dfa5f..39949a1b0a 100644 --- a/flopy/utils/zonbud.py +++ b/flopy/utils/zonbud.py @@ -548,18 +548,18 @@ def _update_budget_recordarray( try: if kstpkper is not None: for rn, cn, flux in zip(rownames, colnames, fluxes): - rowidx = np.where( + rowidx = np.asarray( (self._budget["time_step"] == kstpkper[0]) & (self._budget["stress_period"] == kstpkper[1]) & (self._budget["name"] == rn) - ) + ).nonzero() self._budget[cn][rowidx] += flux elif totim is not None: for rn, cn, flux in zip(rownames, colnames, fluxes): - rowidx = np.where( + rowidx = np.asarray( (self._budget["totim"] == totim) & (self._budget["name"] == rn) - ) + ).nonzero() self._budget[cn][rowidx] += flux except Exception as e: @@ -592,9 +592,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): # ZONE 4 TO 3 IS THE NEGATIVE OF FLOW FROM 3 TO 4. # 1ST, CALCULATE FLOW BETWEEN NODE J,I,K AND J-1,I,K - k, i, j = np.where( + k, i, j = np.asarray( self.izone[:, :, 1:] > self.izone[:, :, :-1] - ) + ).nonzero() # Adjust column values to account for the starting position of "nz" j += 1 @@ -613,9 +613,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): # Don't include CH to CH flow (can occur if CHTOCH option is used) # Create an iterable tuple of (from zone, to zone, flux) # Then group tuple by (from_zone, to_zone) and sum the flux values - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, i, jl] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nzl[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim @@ -625,18 +625,18 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): # Don't include CH to CH flow (can occur if CHTOCH option is used) # Create an iterable tuple of (from zone, to zone, flux) # Then group tuple by (from_zone, to_zone) and sum the flux values - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, i, jl] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nzl[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # FLOW BETWEEN NODE J,I,K AND J+1,I,K - k, i, j = np.where( + k, i, j = np.asarray( self.izone[:, :, :-1] > self.izone[:, :, 1:] - ) + ).nonzero() # Define the zone from which flow is coming nz = self.izone[k, i, j] @@ -652,9 +652,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): # Don't include CH to CH flow (can occur if CHTOCH option is used) # Create an iterable tuple of (from zone, to zone, flux) # Then group tuple by (from_zone, to_zone) and sum the flux values - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, i, jr] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nzr[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim @@ -664,24 +664,24 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): # Don't include CH to CH flow (can occur if CHTOCH option is used) # Create an iterable tuple of (from zone, to zone, flux) # Then group tuple by (from_zone, to_zone) and sum the flux values - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, i, jr] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nzr[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # CALCULATE FLOW TO CONSTANT-HEAD CELLS IN THIS DIRECTION - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = k[j > 0], i[j > 0], j[j > 0] jl = j - 1 nzl = self.izone[k, i, jl] nz = self.izone[k, i, j] q = data[k, i, jl] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, i, jl] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzl[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -689,9 +689,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, i, jl] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzl[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi[tzi != 0]] @@ -699,7 +699,7 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = ( k[j < self.ncol - 1], i[j < self.ncol - 1], @@ -709,9 +709,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): jr = j + 1 nzr = self.izone[k, i, jr] q = data[k, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, i, jr] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzr[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -719,9 +719,9 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, i, jr] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzr[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi] @@ -732,7 +732,6 @@ def _accumulate_flow_frf(self, recname, ich, kstpkper, totim): except Exception as e: print(e) raise - return def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): """ @@ -756,64 +755,64 @@ def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): # "FLOW FRONT FACE" # CALCULATE FLOW BETWEEN NODE J,I,K AND J,I-1,K - k, i, j = np.where( + k, i, j = np.asarray( self.izone[:, 1:, :] < self.izone[:, :-1, :] - ) + ).nonzero() i += 1 ia = i - 1 nza = self.izone[k, ia, j] nz = self.izone[k, i, j] q = data[k, ia, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, ia, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nza[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, ia, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nza[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # CALCULATE FLOW BETWEEN NODE J,I,K AND J,I+1,K. - k, i, j = np.where( + k, i, j = np.asarray( self.izone[:, :-1, :] < self.izone[:, 1:, :] - ) + ).nonzero() nz = self.izone[k, i, j] ib = i + 1 nzb = self.izone[k, ib, j] q = data[k, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, ib, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nzb[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, ib, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # CALCULATE FLOW TO CONSTANT-HEAD CELLS IN THIS DIRECTION - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = k[i > 0], i[i > 0], j[i > 0] ia = i - 1 nza = self.izone[k, ia, j] nz = self.izone[k, i, j] q = data[k, ia, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, ia, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nza[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -821,9 +820,9 @@ def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, ia, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nza[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi] @@ -831,7 +830,7 @@ def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = ( k[i < self.nrow - 1], i[i < self.nrow - 1], @@ -841,9 +840,9 @@ def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): ib = i + 1 nzb = self.izone[k, ib, j] q = data[k, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[k, ib, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -851,9 +850,9 @@ def _accumulate_flow_fff(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[k, ib, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi] @@ -888,64 +887,64 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): # "FLOW LOWER FACE" # CALCULATE FLOW BETWEEN NODE J,I,K AND J,I,K-1 - k, i, j = np.where( + k, i, j = np.asarray( self.izone[1:, :, :] < self.izone[:-1, :, :] - ) + ).nonzero() k += 1 ka = k - 1 nza = self.izone[ka, i, j] nz = self.izone[k, i, j] q = data[ka, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[ka, i, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nza[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[ka, i, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nza[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # CALCULATE FLOW BETWEEN NODE J,I,K AND J,I,K+1 - k, i, j = np.where( + k, i, j = np.asarray( self.izone[:-1, :, :] < self.izone[1:, :, :] - ) + ).nonzero() nz = self.izone[k, i, j] kb = k + 1 nzb = self.izone[kb, i, j] q = data[k, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[kb, i, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nz[idx], nzb[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[kb, i, j] != 1)) - ) + ).nonzero() fzi, tzi, fi = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) self._update_budget_fromfaceflow( fzi, tzi, np.abs(fi), kstpkper, totim ) # CALCULATE FLOW TO CONSTANT-HEAD CELLS IN THIS DIRECTION - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = k[k > 0], i[k > 0], j[k > 0] ka = k - 1 nza = self.izone[ka, i, j] nz = self.izone[k, i, j] q = data[ka, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[ka, i, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nza[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -953,9 +952,9 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[ka, i, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nza[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi] @@ -963,7 +962,7 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - k, i, j = np.where(ich == 1) + k, i, j = np.asarray(ich == 1).nonzero() k, i, j = ( k[k < self.nlay - 1], i[k < self.nlay - 1], @@ -973,9 +972,9 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): kb = k + 1 nzb = self.izone[kb, i, j] q = data[k, i, j] - idx = np.where( + idx = np.asarray( (q > 0) & ((ich[k, i, j] != 1) | (ich[kb, i, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) fz = ["FROM_CONSTANT_HEAD"] * len(tzi) tz = [self._zonenamedict[z] for z in tzi] @@ -983,9 +982,9 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): fz, tz, np.abs(f), kstpkper, totim ) - idx = np.where( + idx = np.asarray( (q < 0) & ((ich[k, i, j] != 1) | (ich[kb, i, j] != 1)) - ) + ).nonzero() fzi, tzi, f = sum_flux_tuples(nzb[idx], nz[idx], q[idx]) fz = ["TO_CONSTANT_HEAD"] * len(fzi) tz = [self._zonenamedict[z] for z in tzi] @@ -996,7 +995,6 @@ def _accumulate_flow_flf(self, recname, ich, kstpkper, totim): except Exception as e: print(e) raise - return def _accumulate_flow_ssst(self, recname, kstpkper, totim): # NOT AN INTERNAL FLOW TERM, SO MUST BE A SOURCE TERM OR STORAGE @@ -1049,9 +1047,9 @@ def _accumulate_flow_ssst(self, recname, kstpkper, totim): # 1-LAYER ARRAY THAT DEFINES LAYER 1 qin = np.ma.zeros(self.cbc_shape, self.float_type) qout = np.ma.zeros(self.cbc_shape, self.float_type) - r, c = np.where(data > 0) + r, c = np.asarray(data > 0).nonzero() qin[0, r, c] = data[r, c] - r, c = np.where(data < 0) + r, c = np.asarray(data < 0).nonzero() qout[0, r, c] = data[r, c] else: # Should not happen @@ -1103,16 +1101,16 @@ def _compute_mass_balance(self, kstpkper, totim): innames = [n for n in recnames if n.startswith("FROM_")] outnames = [n for n in recnames if n.startswith("TO_")] if kstpkper is not None: - rowidx = np.where( + rowidx = np.asarray( (self._budget["time_step"] == kstpkper[0]) & (self._budget["stress_period"] == kstpkper[1]) - & np.isin(self._budget["name"], innames) - ) + & np.in1d(self._budget["name"], innames) + ).nonzero() elif totim is not None: - rowidx = np.where( + rowidx = np.asarray( (self._budget["totim"] == totim) - & np.isin(self._budget["name"], innames) - ) + & np.in1d(self._budget["name"], innames) + ).nonzero() a = _numpyvoid2numeric( self._budget[list(self._zonenamedict.values())][rowidx] ) @@ -1125,16 +1123,16 @@ def _compute_mass_balance(self, kstpkper, totim): # Compute outflows if kstpkper is not None: - rowidx = np.where( + rowidx = np.asarray( (self._budget["time_step"] == kstpkper[0]) & (self._budget["stress_period"] == kstpkper[1]) - & np.isin(self._budget["name"], outnames) - ) + & np.in1d(self._budget["name"], outnames) + ).nonzero() elif totim is not None: - rowidx = np.where( + rowidx = np.asarray( (self._budget["totim"] == totim) - & np.isin(self._budget["name"], outnames) - ) + & np.in1d(self._budget["name"], outnames) + ).nonzero() a = _numpyvoid2numeric( self._budget[list(self._zonenamedict.values())][rowidx] ) @@ -2462,7 +2460,7 @@ def _get_budget(recarray, zonenamedict, names=None, zones=None, net=False): if "totim" in recarray.dtype.names: standard_fields.insert(0, "totim") select_fields = standard_fields + list(zonenamedict.values()) - select_records = np.where(recarray["name"] == recarray["name"]) + select_records = np.asarray(recarray["name"] == recarray["name"]).nonzero() if zones is not None: for idx, z in enumerate(zones): if isinstance(z, int): @@ -2945,10 +2943,10 @@ def _pivot_recarray(recarray): pvt_rec = np.recarray((1,), dtype=dtype) n = 0 for kstp, kper in kstp_kper: - idxs = np.where( + idxs = np.asarray( (recarray["time_step"] == kstp) & (recarray["stress_period"] == kper) - ) + ).nonzero() if len(idxs) == 0: pass else: @@ -3008,7 +3006,7 @@ def _volumetric_flux(recarray, modeltime, extrapolate_kper=False): perlen = modeltime.perlen totim = np.add.accumulate(perlen) for per in range(nper): - idx = np.where(recarray["kper"] == per)[0] + idx = np.asarray(recarray["kper"] == per).nonzero()[0] if len(idx) == 0: continue @@ -3019,7 +3017,7 @@ def _volumetric_flux(recarray, modeltime, extrapolate_kper=False): if zone == 0: continue - zix = np.where(temp["zone"] == zone)[0] + zix = np.asarray(temp["zone"] == zone).nonzero()[0] if len(zix) == 0: raise Exception @@ -3052,9 +3050,9 @@ def _volumetric_flux(recarray, modeltime, extrapolate_kper=False): totim = modeltime.totim for ix, nstp in enumerate(modeltime.nstp): for stp in range(nstp): - idx = np.where( + idx = np.asarray( (recarray["kper"] == ix) & (recarray["kstp"] == stp) - ) + ).nonzero() if len(idx[0]) == 0: continue elif n == 0: