From a4192b04b9e69c2d1fd2a18fd3f13c28c6beaa50 Mon Sep 17 00:00:00 2001 From: mjreno Date: Thu, 12 Feb 2026 17:41:49 -0500 Subject: [PATCH 1/8] fix grb v2 reading verbose logging --- autotest/test_binarygrid_util.py | 16 ++++++++++++++++ examples/data/mfgrd_test/gwf_uzf01.dis.grb | Bin 0 -> 5174 bytes flopy/mf6/utils/binarygrid_util.py | 5 ++++- 3 files changed, 20 insertions(+), 1 deletion(-) create mode 100644 examples/data/mfgrd_test/gwf_uzf01.dis.grb diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index 5d56e03203..1e5c67e5b3 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -31,6 +31,22 @@ def test_mfgrddis_MfGrdFile(mfgrd_test_path): assert isinstance(modelgrid, StructuredGrid) +def test_mfgrddis_MfGrdFile_v2(mfgrd_test_path): + grb = MfGrdFile(mfgrd_test_path / "gwf_uzf01.dis.grb", verbose=True) + nodes = grb.nodes + ia = grb.ia + shape = ia.shape[0] + assert shape == nodes + 1, f"ia size ({shape}) not equal to {nodes + 1}" + + nnz = ia[-1] + ja = grb.ja + shape = ja.shape[0] + assert shape == nnz, f"ja size ({shape}) not equal to {nnz}" + + modelgrid = grb.modelgrid + assert isinstance(modelgrid, StructuredGrid) + + def test_mfgrddis_modelgrid(mfgrd_test_path): fn = mfgrd_test_path / "nwtp3.dis.grb" modelgrid = StructuredGrid.from_binary_grid_file(fn, verbose=True) diff --git a/examples/data/mfgrd_test/gwf_uzf01.dis.grb b/examples/data/mfgrd_test/gwf_uzf01.dis.grb new file mode 100644 index 0000000000000000000000000000000000000000..fd2649528c9a8021b771854e0df7a705921f3936 GIT binary patch literal 5174 zcmeIze`p(Z9LDkMbXspZH?3~YO>b-Mwr;CQyY<()HED87)22zAq`#)q>4qQ(f*=Tj zoFE8-AP9mW2!bF8!axuNAs`5XAP9mW2!bF8!axuNL7!)R7iqO7lh(f)`r_&5`6m75 za(8K=TFgiJVl8U*5zfceTCq}&`dWTPxxQGB(u1u&KZ?Rqy!@I}s>K(D4U)uLj# z9v9+jRL&Qtqg2!xrCa_JO?L{*rR>tC?w}(tS1Swuut(E3cO6gAXH_X}`W}s&!c=c+)uY)TTW397`XBblmJ8KNeKQXu3iENP zx@w1XlwRkU|7-0%a+|$Jy)v`OTSVb_r9QoOla2ctx%Tmf#WlBwPDg5kPe{1-7EDJQ zdjsl7!eYKMoz*?K>WkM0)9N>O^Q$c^`n$cpG!t*$E@7@(i*gg$YBpD2`S;w<%bAcFG9ViZtdkisARF_nlNmA~8*{9a88RRn zb=JuY8IX+{>tu!u$VQcQGD8MrW0rL?Lk47HhIKMS1}lAKc|9PLGBcS>Gm}Y)nM|gb z$z+n5Op458GQmtH1!gii$4n+MGnwR>$t1^2CgaRxl4T~74D(8kF`p&RkfY>ja)cZv zhsZ&4fILN>B>TxeGEJt)Uh+8ELmnZ!$u9B`8Ik+Ry<`Wui)<$Ya^?P5zPD}&g5V~( zMQ)YbK%12Qb5;+t?>VwsRBDa)+XWkJr%1-U4fx8>5xO)KZy5a6H+#P|tqj1*)cgNuFINY6pyI#0^5bjcNmxjAOxa)_zlW=zm z?grp)5blQHZW!)H;O;csjl$g-xH}7XV{n&&yDZ#|!(9&U@^BZ!-8r}`z}*Dg72$3Y z?xx_b1b5SLSBAR^+|9t Date: Thu, 12 Feb 2026 19:35:41 -0500 Subject: [PATCH 2/8] Xadd version and crs properties --- autotest/test_binarygrid_util.py | 12 +++++++- .../{gwf_uzf01.dis.grb => flow_v2.dis.grb} | Bin flopy/mf6/utils/binarygrid_util.py | 27 +++++++++++++++++- 3 files changed, 37 insertions(+), 2 deletions(-) rename examples/data/mfgrd_test/{gwf_uzf01.dis.grb => flow_v2.dis.grb} (100%) diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index 1e5c67e5b3..0dcfbb192e 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -17,6 +17,8 @@ def mfgrd_test_path(example_data_path): def test_mfgrddis_MfGrdFile(mfgrd_test_path): grb = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=True) + assert grb.version == 1 + nodes = grb.nodes ia = grb.ia shape = ia.shape[0] @@ -30,9 +32,13 @@ def test_mfgrddis_MfGrdFile(mfgrd_test_path): modelgrid = grb.modelgrid assert isinstance(modelgrid, StructuredGrid) + assert grb.crs is None + def test_mfgrddis_MfGrdFile_v2(mfgrd_test_path): - grb = MfGrdFile(mfgrd_test_path / "gwf_uzf01.dis.grb", verbose=True) + grb = MfGrdFile(mfgrd_test_path / "flow_v2.dis.grb", verbose=True) + assert grb.version == 2 + nodes = grb.nodes ia = grb.ia shape = ia.shape[0] @@ -46,6 +52,8 @@ def test_mfgrddis_MfGrdFile_v2(mfgrd_test_path): modelgrid = grb.modelgrid assert isinstance(modelgrid, StructuredGrid) + assert grb.crs == "EPSG:26916" + def test_mfgrddis_modelgrid(mfgrd_test_path): fn = mfgrd_test_path / "nwtp3.dis.grb" @@ -78,6 +86,7 @@ def test_mfgrddis_modelgrid(mfgrd_test_path): def test_mfgrddisv_MfGrdFile(mfgrd_test_path): fn = mfgrd_test_path / "flow.disv.grb" grb = MfGrdFile(fn, verbose=True) + assert grb.version == 1 nodes = grb.nodes ia = grb.ia @@ -130,6 +139,7 @@ def test_mfgrddisv_modelgrid(mfgrd_test_path): def test_mfgrddisu_MfGrdFile(mfgrd_test_path): fn = mfgrd_test_path / "keating.disu.grb" grb = MfGrdFile(fn, verbose=True) + assert grb.version == 1 nodes = grb.nodes ia = grb.ia diff --git a/examples/data/mfgrd_test/gwf_uzf01.dis.grb b/examples/data/mfgrd_test/flow_v2.dis.grb similarity index 100% rename from examples/data/mfgrd_test/gwf_uzf01.dis.grb rename to examples/data/mfgrd_test/flow_v2.dis.grb diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 0ed8c235da..ea8e8cf963 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -81,7 +81,7 @@ def __init__(self, filename, precision="double", verbose=False): # version line = self.read_text(self._initial_len).strip() t = line.split() - self._version = t[1] + self._version = int(t[1]) # ntxt line = self.read_text(self._initial_len).strip() @@ -347,6 +347,17 @@ def _get_cellcenters(self): return xycellcenters # properties + @property + def version(self): + """ + MODFLOW 6 grid file version. + + Returns + ------- + version : str + """ + return self._version + @property def grid_type(self): """ @@ -934,3 +945,17 @@ def export(self, filename, precision=None, version=1, verbose=False): if verbose: print(f"Successfully wrote {filename}") + + @property + def crs(self): + """ + CRS user input string (version 2 GRB file only). + + Returns + ------- + crs : str or None + """ + crs = None + if "CRS" in self._datadict: + crs = self._datadict["CRS"] + return crs From c13a2f1789aa9f07235dd8a499b8be7a1f72e3ff Mon Sep 17 00:00:00 2001 From: mjreno Date: Thu, 12 Feb 2026 20:06:51 -0500 Subject: [PATCH 3/8] version is int not str --- flopy/mf6/utils/binarygrid_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index ea8e8cf963..11b0ec1e18 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -354,7 +354,7 @@ def version(self): Returns ------- - version : str + version : int """ return self._version From e0264736eed7aa4858ec0ce5533d5a4c31bc8920 Mon Sep 17 00:00:00 2001 From: mjreno Date: Fri, 26 Jun 2026 13:40:31 -0400 Subject: [PATCH 4/8] support grb v2 export --- autotest/test_binarygrid_util.py | 180 +++++++++++++++++++++++++ flopy/mf6/utils/binarygrid_util.py | 206 ++++++++++++++++------------- 2 files changed, 293 insertions(+), 93 deletions(-) diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index 0dcfbb192e..ce85609a08 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -422,3 +422,183 @@ def test_write_grb_disu_precision_conversion(tmp_path, mfgrd_test_path): # Double precision should match exactly (same precision as original) np.testing.assert_allclose(grb_double.top, grb.top, rtol=1e-12) np.testing.assert_allclose(grb_double.bot, grb.bot, rtol=1e-12) + + +def test_write_grb_v2_roundtrip(tmp_path, mfgrd_test_path): + """Round-trip a version 2 GRB: version and CRS are preserved.""" + grb_orig = MfGrdFile(mfgrd_test_path / "flow_v2.dis.grb", verbose=False) + assert grb_orig.version == 2 + assert grb_orig.crs == "EPSG:26916" + + output_file = tmp_path / "flow_v2_copy.dis.grb" + grb_orig.export(output_file, verbose=False) + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + assert grb_new.nodes == grb_orig.nodes + np.testing.assert_allclose(grb_new.top, grb_orig.top) + np.testing.assert_allclose(grb_new.bot, grb_orig.bot) + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + np.testing.assert_array_equal(grb_new.ja, grb_orig.ja) + + +def test_write_grb_v1_upgrade_to_v2(tmp_path, mfgrd_test_path): + """Exporting a v1 GRB with crs= automatically produces a v2 file.""" + grb_orig = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + assert grb_orig.version == 1 + assert grb_orig.crs is None + + output_file = tmp_path / "nwtp3_v2.dis.grb" + grb_orig.export(output_file, crs="EPSG:26916", verbose=False) + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + assert grb_new.nodes == grb_orig.nodes + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + + +def test_write_grb_v1_force_override(tmp_path, mfgrd_test_path): + """Explicit version=1 suppresses CRS even when crs= is provided.""" + grb = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + output_file = tmp_path / "forced_v1.dis.grb" + grb.export(output_file, version=1, crs="EPSG:26916") + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 1 + assert grb_new.crs is None + assert grb_new.nodes == grb.nodes + + +def test_write_grb_v2_requires_crs(tmp_path, mfgrd_test_path): + """Requesting version=2 without a CRS raises ValueError.""" + grb = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + with pytest.raises(ValueError, match="version=2 requires a CRS"): + grb.export(tmp_path / "out.grb", version=2) + + +def test_write_grb_v2_crs_too_long(tmp_path, mfgrd_test_path): + """A CRS string exceeding 5000 characters raises ValueError.""" + grb = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + with pytest.raises(ValueError, match="exceeds the MODFLOW 6 maximum"): + grb.export(tmp_path / "out.grb", crs="X" * 5001) + + +def test_crs_to_string(): + """_crs_to_string normalizes CRS inputs to strings.""" + from flopy.mf6.utils.binarygrid_util import _crs_to_string + + # str pass-through + assert _crs_to_string("EPSG:26916") == "EPSG:26916" + assert _crs_to_string("some WKT...") == "some WKT..." + + # int → EPSG authority string + assert _crs_to_string(26916) == "EPSG:26916" + + # pyproj.CRS with a known EPSG code → prefers short EPSG form + pyproj = pytest.importorskip("pyproj") + crs_obj = pyproj.CRS.from_epsg(26916) + assert _crs_to_string(crs_obj) == "EPSG:26916" + + # pyproj.CRS built from WKT that has an EPSG code still normalizes + crs_from_wkt = pyproj.CRS.from_wkt(crs_obj.to_wkt()) + assert _crs_to_string(crs_from_wkt) == "EPSG:26916" + + +def test_write_grb_v2_crs_int(tmp_path, mfgrd_test_path): + """Passing an integer EPSG code writes 'EPSG:N' to the GRB.""" + grb_orig = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + output_file = tmp_path / "nwtp3_v2_int.dis.grb" + grb_orig.export(output_file, crs=26916) + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + + +def test_write_grb_v2_crs_pyproj(tmp_path, mfgrd_test_path): + """Passing a pyproj.CRS object normalizes to EPSG string in the GRB.""" + pyproj = pytest.importorskip("pyproj") + grb_orig = MfGrdFile(mfgrd_test_path / "nwtp3.dis.grb", verbose=False) + output_file = tmp_path / "nwtp3_v2_pyproj.dis.grb" + grb_orig.export(output_file, crs=pyproj.CRS.from_epsg(26916)) + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + + +def test_write_grb_disv_v1_upgrade_to_v2(tmp_path, mfgrd_test_path): + """Upgrading a v1 DISV GRB to v2 preserves geometry and writes CRS.""" + grb_orig = MfGrdFile(mfgrd_test_path / "flow.disv.grb", verbose=False) + assert grb_orig.version == 1 + assert grb_orig.crs is None + + output_file = tmp_path / "flow_v2.disv.grb" + grb_orig.export(output_file, crs="EPSG:26916") + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + assert grb_new.grid_type == "DISV" + assert grb_new.nodes == grb_orig.nodes + assert grb_new.nja == grb_orig.nja + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + np.testing.assert_array_equal(grb_new.ja, grb_orig.ja) + np.testing.assert_allclose(grb_new.top, grb_orig.top) + np.testing.assert_allclose(grb_new.bot, grb_orig.bot) + np.testing.assert_allclose( + grb_new._datadict["VERTICES"], grb_orig._datadict["VERTICES"] + ) + np.testing.assert_allclose( + grb_new._datadict["CELLX"], grb_orig._datadict["CELLX"] + ) + np.testing.assert_allclose( + grb_new._datadict["CELLY"], grb_orig._datadict["CELLY"] + ) + + +def test_write_grb_disv_v2_roundtrip(tmp_path, mfgrd_test_path): + """Round-trip a v2 DISV GRB: version and CRS are preserved.""" + # First create a v2 DISV file by upgrading the v1 source + v2_file = tmp_path / "flow_v2.disv.grb" + MfGrdFile(mfgrd_test_path / "flow.disv.grb").export( + v2_file, crs="EPSG:26916" + ) + + grb_orig = MfGrdFile(v2_file, verbose=False) + assert grb_orig.version == 2 + + output_file = tmp_path / "flow_v2_copy.disv.grb" + grb_orig.export(output_file, verbose=False) + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + assert grb_new.nodes == grb_orig.nodes + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + np.testing.assert_allclose( + grb_new._datadict["VERTICES"], grb_orig._datadict["VERTICES"] + ) + + +def test_write_grb_disu_v1_upgrade_to_v2(tmp_path, mfgrd_test_path): + """Upgrading a v1 DISU GRB to v2 preserves connectivity and writes CRS.""" + grb_orig = MfGrdFile(mfgrd_test_path / "flow.disu.grb", verbose=False) + assert grb_orig.version == 1 + assert grb_orig.crs is None + + output_file = tmp_path / "flow_v2.disu.grb" + grb_orig.export(output_file, crs="EPSG:26916") + + grb_new = MfGrdFile(output_file, verbose=False) + assert grb_new.version == 2 + assert grb_new.crs == "EPSG:26916" + assert grb_new.grid_type == "DISU" + assert grb_new.nodes == grb_orig.nodes + assert grb_new.nja == grb_orig.nja + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + np.testing.assert_array_equal(grb_new.ja, grb_orig.ja) + np.testing.assert_allclose(grb_new.top, grb_orig.top) + np.testing.assert_allclose(grb_new.bot, grb_orig.bot) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 11b0ec1e18..5a734b1b82 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -14,6 +14,43 @@ warnings.simplefilter("always", DeprecationWarning) +def _crs_to_string(crs): + """Convert a CRS value to a string suitable for writing to a GRB file. + + Strings are returned unchanged (pass-through for DIS package / GRB + round-trip cases). Integers are treated as EPSG codes. Objects with + pyproj-compatible methods (``to_epsg``, ``to_authority``, ``to_wkt``) are + converted in that priority order so the shortest interoperable form is + preferred. + + Parameters + ---------- + crs : str, int, or pyproj.CRS + CRS value to convert. + + Returns + ------- + str + CRS as a string ready for writing to the GRB CRS field. + """ + if isinstance(crs, str): + return crs + if isinstance(crs, int): + return f"EPSG:{crs}" + # Duck-type pyproj.CRS or any compatible object — no hard dependency. + if hasattr(crs, "to_epsg"): + epsg = crs.to_epsg() + if epsg is not None: + return f"EPSG:{epsg}" + if hasattr(crs, "to_authority"): + authority = crs.to_authority() + if authority is not None: + return ":".join(authority) + if hasattr(crs, "to_wkt"): + return crs.to_wkt() + return str(crs) + + class MfGrdFile(FlopyBinaryData): """ The MfGrdFile class. @@ -761,7 +798,7 @@ def cell2d(self): vertices, cell2d = None, None return vertices, cell2d - def export(self, filename, precision=None, version=1, verbose=False): + def export(self, filename, precision=None, version=None, crs=None, verbose=False): """ Export the binary grid file to a new file. @@ -773,7 +810,14 @@ def export(self, filename, precision=None, version=1, verbose=False): 'single' or 'double'. If None, uses the precision from the original file (default None) version : int, optional - Grid file version (default 1) + Grid file version to write. If None (default), version 2 is + written when a CRS string is available, otherwise version 1. + Pass ``version=1`` to force version 1 even when a CRS is set. + crs : str, optional + CRS user input string to write (e.g. ``"EPSG:26916"`` or OGC + WKT). If None, uses the CRS from the source file if present. + Providing this argument with a version 1 source file upgrades + the output to version 2. verbose : bool, optional Print progress messages (default False) @@ -782,94 +826,60 @@ def export(self, filename, precision=None, version=1, verbose=False): >>> from flopy.mf6.utils import MfGrdFile >>> grb = MfGrdFile('model.dis.grb') >>> grb.export('model_copy.dis.grb') - >>> # Convert to single precision - >>> grb.export('model_single.dis.grb', precision='single') + >>> # Round-trip a v2 file preserving CRS + >>> grb2 = MfGrdFile('model_v2.dis.grb') + >>> grb2.export('model_v2_copy.dis.grb') + >>> # Upgrade a v1 file to v2 with an explicit CRS + >>> grb.export('model_v2.dis.grb', crs='EPSG:26916') """ if precision is None: precision = self.precision - # Build data dictionary from instance - data_dict = {} - for key in self._recordkeys: - if key in ("IA", "JA"): - # Use original 1-based arrays - data_dict[key] = self._datadict[key] - elif key == "TOP": - data_dict[key] = self.top - elif key == "BOTM": - data_dict[key] = self.bot - elif key in self._datadict: - data_dict[key] = self._datadict[key] - - # Define variable metadata based on grid type + raw_crs = crs if crs is not None else self.crs + effective_crs = _crs_to_string(raw_crs) if raw_crs is not None else None + if version is None: + version = 2 if effective_crs else 1 + float_type = "SINGLE" if precision.lower() == "single" else "DOUBLE" - if self.grid_type == "DIS": - var_list = [ - ("NCELLS", "INTEGER", 0, []), - ("NLAY", "INTEGER", 0, []), - ("NROW", "INTEGER", 0, []), - ("NCOL", "INTEGER", 0, []), - ("NJA", "INTEGER", 0, []), - ("XORIGIN", float_type, 0, []), - ("YORIGIN", float_type, 0, []), - ("ANGROT", float_type, 0, []), - ("DELR", float_type, 1, [self.ncol]), - ("DELC", float_type, 1, [self.nrow]), - ("TOP", float_type, 1, [self.nodes]), - ("BOTM", float_type, 1, [self.nodes]), - ("IA", "INTEGER", 1, [self.nodes + 1]), - ("JA", "INTEGER", 1, [self.nja]), - ("IDOMAIN", "INTEGER", 1, [self.nodes]), - ("ICELLTYPE", "INTEGER", 1, [self.nodes]), - ] - elif self.grid_type == "DISV": - # Get dimensions for DISV arrays - nvert = self._datadict["NVERT"] - njavert = self._datadict["NJAVERT"] - var_list = [ - ("NCELLS", "INTEGER", 0, []), - ("NLAY", "INTEGER", 0, []), - ("NCPL", "INTEGER", 0, []), - ("NVERT", "INTEGER", 0, []), - ("NJAVERT", "INTEGER", 0, []), - ("NJA", "INTEGER", 0, []), - ("XORIGIN", float_type, 0, []), - ("YORIGIN", float_type, 0, []), - ("ANGROT", float_type, 0, []), - ("TOP", float_type, 1, [self.nodes]), - ("BOTM", float_type, 1, [self.nodes]), - ("VERTICES", float_type, 2, [nvert, 2]), - ("CELLX", float_type, 1, [self.nodes]), - ("CELLY", float_type, 1, [self.nodes]), - ("IAVERT", "INTEGER", 1, [self.nodes + 1]), - ("JAVERT", "INTEGER", 1, [njavert]), - ("IA", "INTEGER", 1, [self.nodes + 1]), - ("JA", "INTEGER", 1, [self.nja]), - ("IDOMAIN", "INTEGER", 1, [self.nodes]), - ("ICELLTYPE", "INTEGER", 1, [self.nodes]), - ] - elif self.grid_type == "DISU": - var_list = [ - ("NODES", "INTEGER", 0, []), - ("NJA", "INTEGER", 0, []), - ("XORIGIN", float_type, 0, []), - ("YORIGIN", float_type, 0, []), - ("ANGROT", float_type, 0, []), - ("TOP", float_type, 1, [self.nodes]), - ("BOT", float_type, 1, [self.nodes]), - ("IA", "INTEGER", 1, [self.nodes + 1]), - ("JA", "INTEGER", 1, [self.nja]), - ("ICELLTYPE", "INTEGER", 1, [self.nodes]), - ] - # IDOMAIN is optional for DISU - if "IDOMAIN" in self._datadict: - var_list.insert(-1, ("IDOMAIN", "INTEGER", 1, [self.nodes])) - else: - raise NotImplementedError( - f"Grid type {self.grid_type} not yet implemented. " - "Supported grid types: DIS, DISV, DISU" - ) + # Build var_list and data_dict from the parsed record structure so that + # array dimensions come directly from the source file rather than being + # recomputed. This avoids shape mismatches and naturally handles any + # grid type. + var_list = [] + data_dict = {} + for key in self._recordkeys: + if key == "CRS": + continue # handled by the version check below + dt, nd, shp = self._recorddict[key] + if dt == np.int32: + dtype_str = "INTEGER" + elif dt in (np.float32, np.float64): + dtype_str = float_type + elif dt == str: + dtype_str = "CHARACTER" + else: + dtype_str = float_type + # shp is stored in reversed (numpy) order; the write loop reverses + # it back to Fortran (definition-line) order via dims[::-1]. + dims = list(shp) if nd > 0 else [] + var_list.append((key, dtype_str, nd, dims)) + data_dict[key] = self._datadict[key] + + if version >= 2: + if not effective_crs: + raise ValueError( + "version=2 requires a CRS string. Provide crs= or use a " + "version 2 source file." + ) + _MF6_CRS_MAXLEN = 5000 + if len(effective_crs) > _MF6_CRS_MAXLEN: + raise ValueError( + f"CRS string length {len(effective_crs)} exceeds the " + f"MODFLOW 6 maximum of {_MF6_CRS_MAXLEN} characters." + ) + var_list.append(("CRS", "CHARACTER", 1, [len(effective_crs)])) + data_dict["CRS"] = effective_crs ntxt = len(var_list) lentxt = 100 @@ -887,23 +897,31 @@ def export(self, filename, precision=None, version=1, verbose=False): with open(filename, "wb") as f: writer.file = f - # Write text header lines (50 chars each, newline terminated) + # Write text header lines and definition lines in MODFLOW 6 format: + # content left-aligned, space-padded to (nchar-1), then newline at + # position (nchar-1). This matches what MF6 writes and is required + # for Fortran readers that parse these as fixed-length text records. + def write_line(text, nchar): + writer.file.write( + text.encode("ascii").ljust(nchar - 1) + b"\n" + ) + header_len = 50 - writer.write_text(f"GRID {self.grid_type}\n", header_len) - writer.write_text(f"VERSION {version}\n", header_len) - writer.write_text(f"NTXT {ntxt}\n", header_len) - writer.write_text(f"LENTXT {lentxt}\n", header_len) + write_line(f"GRID {self.grid_type}", header_len) + write_line(f"VERSION {version}", header_len) + write_line(f"NTXT {ntxt}", header_len) + write_line(f"LENTXT {lentxt}", header_len) # Write variable definition lines (100 chars each) for name, dtype_str, ndim, dims in var_list: if ndim == 0: - line = f"{name} {dtype_str} NDIM {ndim}\n" + line = f"{name} {dtype_str} NDIM {ndim}" else: dims_str = " ".join( str(d) for d in dims[::-1] ) # Reverse for Fortran order - line = f"{name} {dtype_str} NDIM {ndim} {dims_str}\n" - writer.write_text(line, lentxt) + line = f"{name} {dtype_str} NDIM {ndim} {dims_str}" + write_line(line, lentxt) # Write binary data for each variable for name, dtype_str, ndim, dims in var_list: @@ -924,7 +942,9 @@ def export(self, filename, precision=None, version=1, verbose=False): print(f" Writing {name}") # Write scalar or array data - if ndim == 0: + if dtype_str == "CHARACTER": + writer.write_text(value, len(value)) + elif ndim == 0: # Scalar value if dtype_str == "INTEGER": writer.write_integer(int(value)) From af6ea252d770c1880cf7f7b6b65f55198134f16b Mon Sep 17 00:00:00 2001 From: mjreno Date: Fri, 26 Jun 2026 13:47:06 -0400 Subject: [PATCH 5/8] format --- autotest/test_binarygrid_util.py | 12 +++--------- flopy/mf6/utils/binarygrid_util.py | 13 +++++-------- 2 files changed, 8 insertions(+), 17 deletions(-) diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index ce85609a08..ed5ef09a49 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -551,21 +551,15 @@ def test_write_grb_disv_v1_upgrade_to_v2(tmp_path, mfgrd_test_path): np.testing.assert_allclose( grb_new._datadict["VERTICES"], grb_orig._datadict["VERTICES"] ) - np.testing.assert_allclose( - grb_new._datadict["CELLX"], grb_orig._datadict["CELLX"] - ) - np.testing.assert_allclose( - grb_new._datadict["CELLY"], grb_orig._datadict["CELLY"] - ) + np.testing.assert_allclose(grb_new._datadict["CELLX"], grb_orig._datadict["CELLX"]) + np.testing.assert_allclose(grb_new._datadict["CELLY"], grb_orig._datadict["CELLY"]) def test_write_grb_disv_v2_roundtrip(tmp_path, mfgrd_test_path): """Round-trip a v2 DISV GRB: version and CRS are preserved.""" # First create a v2 DISV file by upgrading the v1 source v2_file = tmp_path / "flow_v2.disv.grb" - MfGrdFile(mfgrd_test_path / "flow.disv.grb").export( - v2_file, crs="EPSG:26916" - ) + MfGrdFile(mfgrd_test_path / "flow.disv.grb").export(v2_file, crs="EPSG:26916") grb_orig = MfGrdFile(v2_file, verbose=False) assert grb_orig.version == 2 diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 5a734b1b82..8f10025757 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -189,10 +189,7 @@ def __init__(self, filename, precision="double", verbose=False): self._datadict[key] = v if self.verbose: - if ( - nd == 0 or - dt == str - ): + if nd == 0 or dt == str: print(f" {key} = {v}") else: print(f" {key}: min = {v.min()} max = {v.max()}") @@ -902,9 +899,7 @@ def export(self, filename, precision=None, version=None, crs=None, verbose=False # position (nchar-1). This matches what MF6 writes and is required # for Fortran readers that parse these as fixed-length text records. def write_line(text, nchar): - writer.file.write( - text.encode("ascii").ljust(nchar - 1) + b"\n" - ) + writer.file.write(text.encode("ascii").ljust(nchar - 1) + b"\n") header_len = 50 write_line(f"GRID {self.grid_type}", header_len) @@ -926,7 +921,9 @@ def write_line(text, nchar): # Write binary data for each variable for name, dtype_str, ndim, dims in var_list: if name not in data_dict: - raise ValueError(f"Required variable '{name}' not found in grid file") + raise ValueError( + f"Required variable '{name}' not found in grid file" + ) value = data_dict[name] From f1680e8f47326347d6ba9e1b1b7bdce032719bc6 Mon Sep 17 00:00:00 2001 From: mjreno Date: Tue, 30 Jun 2026 08:45:48 -0400 Subject: [PATCH 6/8] fix byte order --- flopy/mf6/utils/binarygrid_util.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 8f10025757..f7235b85a5 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -957,8 +957,7 @@ def write_line(text, nchar): elif dtype_str == "SINGLE": arr = arr.astype(np.float32) - # Write array in column-major (Fortran) order - writer.write_record(arr.flatten(order="F"), dtype=arr.dtype) + writer.write_record(arr.flatten(order="C"), dtype=arr.dtype) if verbose: print(f"Successfully wrote {filename}") From 9641378ae32c6014e567104e1f4a7a2a927fe6f2 Mon Sep 17 00:00:00 2001 From: mjreno Date: Tue, 30 Jun 2026 09:21:53 -0400 Subject: [PATCH 7/8] add disv compare test --- autotest/test_binarygrid_util.py | 33 ++++++++++++++++++ .../data/mfgrd_test/disv_layered.disv.grb | Bin 0 -> 3988 bytes 2 files changed, 33 insertions(+) create mode 100644 examples/data/mfgrd_test/disv_layered.disv.grb diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index ed5ef09a49..6c88fd4c86 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -190,11 +190,13 @@ def test_mfgrddisu_modelgrid(mfgrd_test_path): def test_write_grb_instance_method(tmp_path, mfgrd_test_path): original_file = mfgrd_test_path / "nwtp3.dis.grb" grb_orig = MfGrdFile(original_file, verbose=False) + _ = grb_orig.modelgrid # trigger _set_modelgrid so BOTM is reshaped before export output_file = tmp_path / "test_instance.dis.grb" grb_orig.export(output_file, verbose=False) grb_new = MfGrdFile(output_file, verbose=False) + _ = grb_new.modelgrid assert grb_new.grid_type == grb_orig.grid_type assert grb_new.nodes == grb_orig.nodes @@ -217,6 +219,37 @@ def test_write_grb_instance_method(tmp_path, mfgrd_test_path): np.testing.assert_array_equal(grb_new.idomain, grb_orig.idomain) +def test_write_grb_instance_method_disv(tmp_path, mfgrd_test_path): + # disv_layered.disv.grb has distinct BOTM values per layer (-10/-20/-30), + # which distinguishes C from F flatten order and guards against order regression. + original_file = mfgrd_test_path / "disv_layered.disv.grb" + grb_orig = MfGrdFile(original_file, verbose=False) + _ = grb_orig.modelgrid # trigger _set_modelgrid so BOTM is reshaped before export + + output_file = tmp_path / "test_instance.disv.grb" + grb_orig.export(output_file, verbose=False) + + grb_new = MfGrdFile(output_file, verbose=False) + _ = grb_new.modelgrid + + assert grb_new.grid_type == grb_orig.grid_type + assert grb_new.nodes == grb_orig.nodes + assert grb_new.nlay == grb_orig.nlay + assert grb_new.ncpl == grb_orig.ncpl + assert grb_new.nja == grb_orig.nja + + np.testing.assert_allclose(grb_new.xorigin, grb_orig.xorigin) + np.testing.assert_allclose(grb_new.yorigin, grb_orig.yorigin) + np.testing.assert_allclose(grb_new.angrot, grb_orig.angrot) + + np.testing.assert_allclose(grb_new.top, grb_orig.top) + np.testing.assert_allclose(grb_new.bot, grb_orig.bot) + + np.testing.assert_array_equal(grb_new.ia, grb_orig.ia) + np.testing.assert_array_equal(grb_new.ja, grb_orig.ja) + np.testing.assert_array_equal(grb_new.idomain, grb_orig.idomain) + + def test_write_grb_instance_method_precision_conversion(tmp_path, mfgrd_test_path): original_file = mfgrd_test_path / "nwtp3.dis.grb" grb = MfGrdFile(original_file, verbose=False) diff --git a/examples/data/mfgrd_test/disv_layered.disv.grb b/examples/data/mfgrd_test/disv_layered.disv.grb new file mode 100644 index 0000000000000000000000000000000000000000..069f9ab4fe5abe130d2903937edc3511cfd78912 GIT binary patch literal 3988 zcmchY%Wl&^6o$hchD*7W`+||$MOjeAra+MpC~d6BRcc3xTO}J72#FON76=ImHoOC` z!&~q={9o`wr|ISeP8i$`%!-or-=_L7I>Fr!_0@x zs`=i#ZZxh5npKNjvK@8110N?@)QS2&X~(<1=~sOCsygkAOqX<9qe+J6Gs@-)OGmFe z$?y4EOTrs98*+;HOlaAnD*mNocMOSx7Cdf zzAHLwt~U$MvE3dWbfZ(?XeFI~nk`a>cUjsi5Cp!~@0J~OTWPjiBB}9+k6#iuTo!Lf zg91o!0zvJ&HRc#Cm?NHY%mW`(y|<4(6$cSF?J5a;eHRJeRb0AC0v~Q3 zv{1`PI*Qko_Nms(YOJH3?zVJa75)e%roO*H?G+zyHtI!-V{wlgjmA0Qif~hSCj1ai z=UM&ppXa1@BVXsF_HDk-$?UJczgyDvEw^hWzgO}t&uk-)Tb|XMU$$0f^Hy(lWovmh zZ}sL|z1i}-)mz-^&0o&@f0Nsfxt(|3>dMyEEn8c^V9)49TokSfeBXD4N5T_fRd_0_ z3*?;@$h~mL=WJZc=ZWK=7cL9To`XQhAQa(={ zpE@)!k5A9^0$7XqE#aPE8do%46FAqQuqnI~B4J19346jwcrSbqJ_(`SYwO(B7=HNJ1Gw(wW9I1x zuH}3TpY0{@o_zSoB@a&ai!r>mCcVLZUtnxn$)kq{0=1@>+=l|z^b#kRcY!s{ Date: Tue, 30 Jun 2026 10:17:29 -0400 Subject: [PATCH 8/8] cleanup --- flopy/mf6/utils/binarygrid_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index f7235b85a5..360a4bdaab 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -957,7 +957,7 @@ def write_line(text, nchar): elif dtype_str == "SINGLE": arr = arr.astype(np.float32) - writer.write_record(arr.flatten(order="C"), dtype=arr.dtype) + writer.write_record(np.ravel(arr, order="C"), dtype=arr.dtype) if verbose: print(f"Successfully wrote {filename}")