Skip to content

Commit 6bb9c1a

Browse files
authored
fix(HeadUFile): disable get_alldata without layer (#2686)
Close #1502. The issue is that get_alldata() assumes regular structure that USG files don't have. We could return e.g. a list of arrays or some other ragged representation but I figure better to acknowledge the limitation clearly than to return something technically valid but confusing. So, block usage of the method without providing a layer.
1 parent 5e64da8 commit 6bb9c1a

2 files changed

Lines changed: 76 additions & 0 deletions

File tree

autotest/test_headufile.py

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,3 +149,36 @@ def get_expected():
149149
expected = get_expected()
150150
for layer, nn in lni:
151151
assert expected[layer][nn] == head[layer][nn]
152+
153+
154+
@requires_exe("mfusg", "gridgen")
155+
@requires_pkg("shapely", "shapefile")
156+
def test_get_alldata_not_supported(mfusg_model):
157+
# test for https://github.com/modflowpy/flopy/issues/1502
158+
# get_alldata should raise NotImplementedError when mflay=None
159+
160+
model, head_file = mfusg_model
161+
162+
# Should raise NotImplementedError when mflay=None (all layers)
163+
with pytest.raises(NotImplementedError) as excinfo:
164+
head_file.get_alldata()
165+
166+
assert "mflay=None" in str(excinfo.value)
167+
assert "not supported" in str(excinfo.value).lower()
168+
169+
# But should work fine when mflay is specified (single layer)
170+
import numpy as np
171+
172+
alldata = head_file.get_alldata(mflay=0)
173+
174+
# Verify it returns a numpy array
175+
assert isinstance(alldata, np.ndarray)
176+
177+
# Verify shape is (ntimes, ncells_in_layer)
178+
times = head_file.get_times()
179+
assert alldata.shape[0] == len(times)
180+
181+
# Verify the data matches what get_data returns for each timestep
182+
for i, time in enumerate(times):
183+
data = head_file.get_data(totim=time, mflay=0)
184+
np.testing.assert_array_almost_equal(alldata[i], data)

flopy/utils/binaryfile/__init__.py

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -863,6 +863,49 @@ def get_ts(self, idx):
863863

864864
return np.array(result)
865865

866+
def get_alldata(self, mflay=None, nodata=-9999):
867+
"""
868+
Get all data from the USG head file.
869+
870+
Parameters
871+
----------
872+
mflay : integer
873+
MODFLOW zero-based layer number to return. For USG files, this
874+
parameter is required (cannot be None). (Default is None.)
875+
876+
nodata : float
877+
The nodata value in the data array. All array values that have the
878+
nodata value will be assigned np.nan.
879+
880+
Returns
881+
-------
882+
data : numpy array
883+
Array has size (ntimes, ncells_in_layer) when mflay is specified.
884+
885+
Raises
886+
------
887+
NotImplementedError
888+
Raised when mflay=None. USG head files contain ragged arrays with
889+
variable-sized data per layer, which cannot be converted to a
890+
uniform numpy array when retrieving all layers.
891+
892+
Notes
893+
-----
894+
For MODFLOW-USG files, the mflay parameter must be specified to retrieve
895+
data for a single layer across all timesteps. To get all layers for a
896+
specific timestep, use get_data() instead.
897+
898+
"""
899+
if mflay is None:
900+
raise NotImplementedError(
901+
"get_alldata() with mflay=None is not supported for"
902+
"MODFLOW-USG head files. These contain variably-size "
903+
"data per layer which cannot be stacked into a numpy "
904+
"array. Specify mflay to get data for a single layer "
905+
"or use get_data() for specific timesteps."
906+
)
907+
return super().get_alldata(mflay=mflay, nodata=nodata)
908+
866909

867910
class BudgetIndexError(Exception):
868911
pass

0 commit comments

Comments
 (0)