Parse MODFLOW6 input (partially)

imod/mf6/__init__.py

@@ -4,6 +4,7 @@
 
 from imod.mf6.chd import ConstantHead
 from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disu import LowLevelUnstructuredDiscretization
 from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.drn import Drainage
 from imod.mf6.evt import Evapotranspiration
@@ -18,7 +19,7 @@
 from imod.mf6.model import GroundwaterFlowModel
 from imod.mf6.npf import NodePropertyFlow
 from imod.mf6.oc import OutputControl
-from imod.mf6.out import open_cbc, open_hds, read_cbc_headers, read_grb
+from imod.mf6.out import open_cbc, open_hds, open_hds_like, read_cbc_headers, read_grb
 from imod.mf6.rch import Recharge
 from imod.mf6.riv import River
 from imod.mf6.simulation import Modflow6Simulation

imod/mf6/dis.py

@@ -1,8 +1,13 @@
+from pathlib import Path
+
 import numpy as np
+import xarray as xr
 
 import imod
 from imod.mf6.pkgbase import Package, VariableMetaData
 
+from .read_input import read_dis_blockfile
+
 
 class StructuredDiscretization(Package):
     """
@@ -60,6 +65,40 @@ def _delrc(self, dx):
         else:
             raise ValueError(f"Unhandled type of {dx}")
 
+    @classmethod
+    def open(cls, path: str, simroot: Path) -> "StructuredDiscretization":
+        content = read_dis_blockfile(simroot / path, simroot)
+        nlay = content["nlay"]
+
+        griddata = content["griddata"]
+        xmin = float(content.get("xorigin", 0.0))
+        ymin = float(content.get("yorigin", 0.0))
+        dx = griddata.pop("delr")
+        dy = griddata.pop("delc")
+        x = (xmin + np.cumsum(dx)) - 0.5 * dx
+        y = ((ymin + np.cumsum(dy)) - 0.5 * dy)[::-1]
+
+        # Top is a 2D array unlike the others
+        top = xr.DataArray(
+            data=griddata.pop("top"),
+            coords={"y": y, "x": x},
+            dims=("y", "x"),
+        )
+
+        coords = {"layer": np.arange(1, nlay + 1), "y": y, "x": x}
+        dims = ("layer", "y", "x")
+        inverse_keywords = {v: k for k, v in cls._keyword_map.items()}
+        variables = {"top": top}
+        for key, value in griddata.items():
+            invkey = inverse_keywords.get(key, key)
+            variables[invkey] = xr.DataArray(
+                value,
+                coords,
+                dims,
+            )
+
+        return cls(**variables)
+
     def render(self, directory, pkgname, globaltimes, binary):
         disdirectory = directory / pkgname
         d = {}