transpose dis arrays

Author: mjreno

docs/examples/quickstart.py

@@ -33,7 +33,7 @@
 )
 
 # model interface
-gwfif = Flopy3Model(gwf)
+fp3gwf = Flopy3Model(gwf)
 
 # check CHD
 assert chd.data["head"][0, 0] == 1.0
@@ -55,15 +55,21 @@
 # set data
 bobj = flopy.utils.CellBudgetFile(bpth, precision="double")
 spdis = bobj.get_data(text="DATA-SPDIS")[0]
-# spdis = imod.mf6.open_cbc(bpth, grbpth, merge_to_dataset=True)
+spdisimod = imod.mf6.open_cbc(bpth, grbpth, merge_to_dataset=True)
+print(spdis)
+print(type(spdisimod))
+# print(spdisimod)
+for d in spdisimod.variables:
+    print(d)
+    print(spdisimod.variables[d].data)
 heads = imod.mf6.open_hds(hpth, grbpth)
 
 # discharge vectors
-qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwfif)
+qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, fp3gwf)
 
 # plot
 fig, ax = plt.subplots()
-pmv = flopy.plot.PlotMapView(model=gwfif, ax=ax)
+pmv = flopy.plot.PlotMapView(model=fp3gwf, ax=ax)
 pmv.plot_array(heads[0][0])
 pmv.plot_grid(colors="white")
 pmv.contour_array(heads[0][0], levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)

flopy4/mf6/gwf/dis.py

@@ -58,19 +58,19 @@ class Dis(Package):
     )
     top: NDArray[np.floating] = array(
         default=1.0,
-        dims=("ncol", "nrow"),
+        dims=("nrow", "ncol"),
         metadata={"block": "griddata"},
         converter=Converter(convert_array, takes_self=True, takes_field=True),
     )
     botm: NDArray[np.floating] = array(
         default=0.0,
-        dims=("ncol", "nrow", "nlay"),
+        dims=("nlay", "nrow", "ncol"),
         metadata={"block": "griddata"},
         converter=Converter(convert_array, takes_self=True, takes_field=True),
     )
     idomain: NDArray[np.integer] = array(
         default=1,
-        dims=("ncol", "nrow", "nlay"),
+        dims=("nlay", "nrow", "ncol"),
         metadata={"block": "griddata"},
         converter=Converter(convert_array, takes_self=True, takes_field=True),
     )

flopy4/mf6/interface/flopy3.py

@@ -1,15 +1,14 @@
 import re
-from typing import Optional
+from typing import Optional, Union
 
-import numpy as np
 from flopy.datbase import DataInterface, DataListInterface, DataType
 from flopy.discretization import StructuredGrid
 from flopy.discretization.grid import Grid
 from flopy.discretization.modeltime import ModelTime
 from flopy.mbase import ModelInterface
 from flopy.pakbase import PackageInterface
 from flopy.plot.plotutil import PlotUtilities
-from xattree import _get_xatspec
+from xattree import _Array, _Attr, _get_xatspec
 
 from flopy4.mf6.model import Model
 from flopy4.mf6.package import Package
@@ -51,9 +50,9 @@ def __init__(
                 self._grid = StructuredGrid(
                     delc=model.dis.delc.data,
                     delr=model.dis.delr.data,
-                    top=np.transpose(model.dis.top.data),
-                    botm=np.transpose(model.dis.botm.data),
-                    idomain=np.transpose(model.dis.idomain.data),
+                    top=model.dis.top.data,
+                    botm=model.dis.botm.data,
+                    idomain=model.dis.idomain.data,
                     lenuni=lenuni,
                     crs=None,
                     prjfile=None,
@@ -104,8 +103,8 @@ def solver_tols(self):
         """
         Solver inner hclose and rclose values.
         """
-        assert self._ims
-        return self._ims.inner_hclose, self._ims.inner_rclose
+        if self._ims is not None:
+            return self._ims.inner_hclose, self._ims.inner_rclose
 
         return None
 
@@ -154,8 +153,8 @@ def get_package_list(self, ftype=None):
         return [p.name for p in self._plist]
 
     def plot(self, packages: Optional[list] = None, **kwargs):
-        if packages is None:
-            packages = self.get_package_list()
+        # if packages is None:
+        #    packages = self.get_package_list()
         return PlotUtilities._plot_model_helper(
             self, SelPackList=packages, **kwargs
         )
@@ -190,25 +189,25 @@ def __init__(
                 and self._spec.flat[a].type is not None
             ):
                 d_fp3 = Flopy3Data(
+                    data=self._data.attrs[a],
+                    spec=self._spec.flat[a],
                     name=a,
                     modelname=self.parent,
                     modelgrid=self._grid,
                     modeltime=modeltime,
-                    data=self._data.attrs[a],
-                    spec=self._spec.flat[a],
                 )
                 self.__dict__[f"{a}"] = d_fp3
                 self._dlist.append(d_fp3)
 
         for v in self._data.data_vars:
             if self._data.data_vars[v] is not None:
                 d_fp3 = Flopy3Data(
+                    data=self._data.data_vars[v],
+                    spec=self._spec.flat[v],
                     name=v,
                     modelname=self.parent,
                     modelgrid=self._grid,
                     modeltime=modeltime,
-                    data=self._data.data_vars[v],
-                    spec=self._spec.flat[v],
                 )
                 self.__dict__[f"{v}"] = d_fp3
                 self._dlist.append(d_fp3)
@@ -257,22 +256,18 @@ def check(self, f=None, verbose=True, level=1, checktype=None):
         return None
 
     def plot(self, **kwargs):
-        # kwargs = {}
-        # kwargs["filename_base"] = "modelif"
-
         return PlotUtilities._plot_package_helper(self, **kwargs)
 
 
 class Flopy3Data(DataInterface):
     def __init__(
-        # TODO: types of data and spec are unions
         self,
+        data,
+        spec: Union[_Attr, _Array],
         name: Optional[str] = None,
         modelname: Optional[str] = None,
         modelgrid: Optional[Grid] = None,
         modeltime: Optional[ModelTime] = None,
-        data=None,
-        spec=None,
     ):
         assert data is not None
         assert spec is not None
@@ -326,14 +321,12 @@ def data_type(self):
 
     @property
     def dtype(self):
-        # return self._spec.type
         return self._spec.type.__name__
 
     @property
     def array(self):
         if self._spec.type.__name__ == "ndarray":
             if "nnodes" in self._data.dims:
-                # TODO: transpose?
                 if "nper" in self._data.dims:
                     shape = (
                         self._time.nper,
@@ -348,7 +341,7 @@ def array(self):
 
                 return self._data.data.reshape(shape)
             else:
-                return np.transpose(self._data.data)
+                return self._data.data
         return None
 
     @property