Skip to content

Commit

Permalink
update(Grid): refactor thick and saturated_thick: (#1768)
Browse files Browse the repository at this point in the history
* update(Grid): refactor thick and saturated_thick:

* Deprecate thick and saturated_thick
* Added cell_thickness and saturated_thickness to Grid 
* changed calls to thick and saturated_thick to cell_thickness and saturated_thickness

* refactor sat_thk to saturated_thickness
  • Loading branch information
jlarsen-usgs committed Apr 15, 2023
1 parent baa322a commit 945be7c
Show file tree
Hide file tree
Showing 13 changed files with 149 additions and 100 deletions.
74 changes: 37 additions & 37 deletions autotest/test_grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -909,84 +909,84 @@ def unstructured_grid():
return GridCases().unstructured_small()


def test_structured_thick(structured_grid):
thick = structured_grid.thick
assert np.allclose(thick, 5.0), "thicknesses != 5."
def test_structured_thickness(structured_grid):
thickness = structured_grid.cell_thickness
assert np.allclose(thickness, 5.0), "thicknesses != 5."

sat_thick = structured_grid.saturated_thick(structured_grid.botm + 10.0)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
sat_thick = structured_grid.saturated_thickness(structured_grid.botm + 10.0)
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = structured_grid.saturated_thick(structured_grid.botm + 5.0)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
sat_thick = structured_grid.saturated_thickness(structured_grid.botm + 5.0)
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = structured_grid.saturated_thick(structured_grid.botm + 2.5)
sat_thick = structured_grid.saturated_thickness(structured_grid.botm + 2.5)
assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"

sat_thick = structured_grid.saturated_thick(structured_grid.botm)
sat_thick = structured_grid.saturated_thickness(structured_grid.botm)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."

sat_thick = structured_grid.saturated_thick(structured_grid.botm - 100.0)
sat_thick = structured_grid.saturated_thickness(structured_grid.botm - 100.0)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."


def test_vertices_thick(vertex_grid):
thick = vertex_grid.thick
assert np.allclose(thick, 5.0), "thicknesses != 5."
def test_vertices_thickness(vertex_grid):
thickness = vertex_grid.cell_thickness
assert np.allclose(thickness, 5.0), "thicknesses != 5."

sat_thick = vertex_grid.saturated_thick(vertex_grid.botm + 10.0)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
sat_thick = vertex_grid.saturated_thickness(vertex_grid.botm + 10.0)
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = vertex_grid.saturated_thick(vertex_grid.botm + 5.0)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
sat_thick = vertex_grid.saturated_thickness(vertex_grid.botm + 5.0)
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = vertex_grid.saturated_thick(vertex_grid.botm + 2.5)
sat_thick = vertex_grid.saturated_thickness(vertex_grid.botm + 2.5)
assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"

sat_thick = vertex_grid.saturated_thick(vertex_grid.botm)
sat_thick = vertex_grid.saturated_thickness(vertex_grid.botm)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."

sat_thick = vertex_grid.saturated_thick(vertex_grid.botm - 100.0)
sat_thick = vertex_grid.saturated_thickness(vertex_grid.botm - 100.0)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."


def test_unstructured_thick(unstructured_grid):
thick = unstructured_grid.thick
assert np.allclose(thick, 5.0), "thicknesses != 5."
def test_unstructured_thickness(unstructured_grid):
thickness = unstructured_grid.cell_thickness
assert np.allclose(thickness, 5.0), "thicknesses != 5."

sat_thick = unstructured_grid.saturated_thick(
sat_thick = unstructured_grid.saturated_thickness(
unstructured_grid.botm + 10.0
)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = unstructured_grid.saturated_thick(unstructured_grid.botm + 5.0)
assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm + 5.0)
assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."

sat_thick = unstructured_grid.saturated_thick(unstructured_grid.botm + 2.5)
sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm + 2.5)
assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"

sat_thick = unstructured_grid.saturated_thick(unstructured_grid.botm)
sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."

sat_thick = unstructured_grid.saturated_thick(
sat_thick = unstructured_grid.saturated_thickness(
unstructured_grid.botm - 100.0
)
assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."


@parametrize_with_cases("grid", cases=GridCases, prefix="structured_cbd")
def test_structured_ncb_thick(grid):
thick = grid.thick
def test_structured_ncb_thickness(grid):
thickness = grid.cell_thickness

assert thick.shape[0] == grid.nlay + np.count_nonzero(
assert thickness.shape[0] == grid.nlay + np.count_nonzero(
grid.laycbd
), "grid thick attribute returns incorrect shape"
), "grid cell_thickness attribute returns incorrect shape"

thick = grid.remove_confining_beds(grid.thick)
thickness = grid.remove_confining_beds(grid.cell_thickness)
assert (
thick.shape == grid.shape
thickness.shape == grid.shape
), "quasi3d confining beds not properly removed"

sat_thick = grid.saturated_thick(grid.thick)
sat_thick = grid.saturated_thickness(grid.cell_thickness)
assert (
sat_thick.shape == grid.shape
), "saturated_thickness confining beds not removed"
Expand Down
8 changes: 4 additions & 4 deletions autotest/test_modflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -316,11 +316,11 @@ def test_load_twri_grid(example_data_path):
assert (
mg.shape == shape
), f"modelgrid shape {mg.shape} not equal to {shape}"
thick = mg.thick
thickness = mg.cell_thickness
shape = (5, 15, 15)
assert (
thick.shape == shape
), f"thickness shape {thick.shape} not equal to {shape}"
thickness.shape == shape
), f"cell_thickness shape {thickness.shape} not equal to {shape}"


def test_mg(function_tmpdir):
Expand All @@ -347,7 +347,7 @@ def test_mg(function_tmpdir):
botm=botm,
)
bas = ModflowBas(ms, ifrefm=True)
t = ms.modelgrid.thick
t = ms.modelgrid.cell_thickness

# test instantiation of an empty basic Structured Grid
mg = StructuredGrid(dis.delc.array, dis.delr.array)
Expand Down
2 changes: 1 addition & 1 deletion autotest/test_postprocessing.py
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,7 @@ def test_get_sat_thickness_gradients(function_tmpdir):
dz = np.array([np.nan, -0.9])
assert np.nansum(np.abs(dh / dz - grad[:, 1, 0])) < 1e-6

sat_thick = m.modelgrid.saturated_thick(hds, mask=nodata)
sat_thick = m.modelgrid.saturated_thickness(hds, mask=nodata)
assert (
np.abs(np.sum(sat_thick[:, 1, 1] - np.array([0.2, 1.0, 1.0]))) < 1e-6
), "failed saturated thickness comparison (grid.thick())"
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,7 @@
"metadata": {},
"source": [
"### Get the saturated thickness of a layer\n",
"`m.modelgrid.saturated_thick()` returns an nlay, nrow, ncol array of saturated thicknesses."
"`m.modelgrid.saturated_thickness()` returns an nlay, nrow, ncol array of saturated thicknesses."
]
},
{
Expand All @@ -399,7 +399,7 @@
}
],
"source": [
"st = m.modelgrid.saturated_thick(hds, mask=-9999.0)\n",
"st = m.modelgrid.saturated_thickness(hds, mask=-9999.0)\n",
"\n",
"plt.imshow(st[0])\n",
"plt.colorbar(label=\"Saturated thickness\")\n",
Expand Down Expand Up @@ -633,7 +633,7 @@
}
],
"source": [
"m.modelgrid.thick[:, r, c]"
"m.modelgrid.cell_thickness[:, r, c]"
]
},
{
Expand Down Expand Up @@ -674,7 +674,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
"version": "3.9.12"
}
},
"nbformat": 4,
Expand Down
13 changes: 3 additions & 10 deletions examples/Notebooks/flopy3_SEAWAT_henry_problem.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2023-02-22T02:45:17.582730Z",
"iopub.status.busy": "2023-02-22T02:45:17.582364Z",
Expand Down Expand Up @@ -474,13 +473,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Util2d:hk layer 1: resetting 'how' to external\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Util2d:hk layer 1: resetting 'how' to external\n",
" Elapsed run time: 8.265 Seconds\n",
"\n",
" Normal termination of SEAWAT\n"
Expand Down Expand Up @@ -518,7 +511,7 @@
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -532,7 +525,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
"version": "3.9.12"
}
},
"nbformat": 4,
Expand Down
20 changes: 9 additions & 11 deletions examples/Notebooks/flopy3_WatertableRecharge_example.ipynb

Large diffs are not rendered by default.

14 changes: 7 additions & 7 deletions examples/Notebooks/flopy3_demo_of_modelgrid_classes.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -1290,7 +1290,7 @@
" - `nnodes` : returns the number of cells in the model\n",
" - `top` : returns an array of the top elevation of the model or for unstructured grids the top elevation of all cells\n",
" - `botm` : retuns an array of the botm elevations of model cells\n",
" - `thick` : returns an array of model cell thickness\n",
" - `cell_thickness` : returns an array of model cell thickness\n",
" - `idomain` : retruns the idomain array associated with the model\n",
" - `xvertices` : returns an array of x-vertices for model cells\n",
" - `yvertices` : returns an array of y-vertices for model cells\n",
Expand Down Expand Up @@ -1409,7 +1409,7 @@
" nrows=1, ncols=3, figsize=(18, 6), subplot_kw={\"aspect\": \"equal\"}\n",
")\n",
"\n",
"arrays = [modelgrid.top, modelgrid.botm, modelgrid.thick]\n",
"arrays = [modelgrid.top, modelgrid.botm, modelgrid.cell_thickness]\n",
"\n",
"labels = [\"top\", \"botm\", \"thickness\"]\n",
"\n",
Expand Down Expand Up @@ -1450,7 +1450,7 @@
" - `get_coords()` : Method to convert model coordinates to real-world coordinates based on coordinate reference information\n",
" - `get_local_coords()` : Method to convert real-world coordinates to model coordinates based on coordinate reference information\n",
" - `intersect()` : Method to get the cellid (`StructuredGrid`=(row, column) OR `VertexGrid` & `UnstrucuturedGrid`=node number) from either model coordinates or from real-world coordinates\n",
" - `saturated_thick()` : Method to get the saturated thickness\n",
" - `saturated_thickness()` : Method to get the saturated thickness\n",
" - `write_shapefile()` : Method to write a shapefile of the grid with just the cellid attributes"
]
},
Expand Down Expand Up @@ -1715,7 +1715,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### saturated_thick()\n",
"#### saturated_thickness()\n",
"\n",
"Method to calculate the saturated thickness of each model cell from an array of elevations. Input parameters to this method include:\n",
"\n",
Expand Down Expand Up @@ -1802,10 +1802,10 @@
"head = hds.get_alldata()[0]\n",
"\n",
"# calculate saturated thickness\n",
"sat_thk = ml.modelgrid.saturated_thick(head, mask=[1e30])\n",
"sat_thk = ml.modelgrid.saturated_thickness(head, mask=[1e30])\n",
"\n",
"# get thick for comparison\n",
"thk = ml.modelgrid.thick\n",
"thk = ml.modelgrid.cell_thickness\n",
"\n",
"# plot thickness and saturated thickness for comparison\n",
"fig, ax = plt.subplots(\n",
Expand Down Expand Up @@ -1898,7 +1898,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.2"
"version": "3.9.12"
}
},
"nbformat": 4,
Expand Down
66 changes: 55 additions & 11 deletions flopy/discretization/grid.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import copy
import os
import warnings

import numpy as np

Expand Down Expand Up @@ -299,7 +300,7 @@ def laycbd(self):
return self._laycbd

@property
def thick(self):
def cell_thickness(self):
"""
Get the cell thickness for a structured, vertex, or unstructured grid.
Expand All @@ -309,7 +310,24 @@ def thick(self):
"""
return -np.diff(self.top_botm, axis=0).reshape(self._botm.shape)

def saturated_thick(self, array, mask=None):
@property
def thick(self):
"""
DEPRECATED method. thick will be removed in version 3.3.9
Get the cell thickness for a structured, vertex, or unstructured grid.
Returns
-------
thick : calculated thickness
"""
warnings.warn(
"thick has been replaced with cell_thickness and will be removed in version 3.3.9,",
DeprecationWarning,
)
return self.cell_thickness

def saturated_thickness(self, array, mask=None):
"""
Get the saturated thickness for a structured, vertex, or unstructured
grid. If the optional array is passed then thickness is returned
Expand All @@ -325,26 +343,52 @@ def saturated_thick(self, array, mask=None):
Returns
-------
thick : calculated saturated thickness
thickness : calculated saturated thickness
"""
thick = self.thick
top = self.top_botm[:-1].reshape(thick.shape)
bot = self.top_botm[1:].reshape(thick.shape)
thick = self.remove_confining_beds(thick)
thickness = self.cell_thickness
top = self.top_botm[:-1].reshape(thickness.shape)
bot = self.top_botm[1:].reshape(thickness.shape)
thickness = self.remove_confining_beds(thickness)
top = self.remove_confining_beds(top)
bot = self.remove_confining_beds(bot)
array = self.remove_confining_beds(array)

idx = np.where((array < top) & (array > bot))
thick[idx] = array[idx] - bot[idx]
thickness[idx] = array[idx] - bot[idx]
idx = np.where(array <= bot)
thick[idx] = 0.0
thickness[idx] = 0.0
if mask is not None:
if isinstance(mask, (float, int)):
mask = [float(mask)]
for mask_value in mask:
thick[np.where(array == mask_value)] = np.nan
return thick
thickness[np.where(array == mask_value)] = np.nan
return thickness

def saturated_thick(self, array, mask=None):
"""
DEPRECATED method. saturated_thick will be removed in version 3.3.9
Get the saturated thickness for a structured, vertex, or unstructured
grid. If the optional array is passed then thickness is returned
relative to array values (saturated thickness). Returned values
ranges from zero to cell thickness if optional array is passed.
Parameters
----------
array : ndarray
array of elevations that will be used to adjust the cell thickness
mask: float, list, tuple, ndarray
array values to replace with a nan value.
Returns
-------
thick : calculated saturated thickness
"""
warnings.warn(
"saturated_thick has been replaced with saturated_thickness and will be removed in version 3.3.9,",
DeprecationWarning,
)
return self.saturated_thickness(array, mask)

@property
def units(self):
Expand Down
Loading

0 comments on commit 945be7c

Please sign in to comment.