alignment with mf6

Author: mjreno

autotest/regression/test_model_netcdf.py

@@ -32,13 +32,9 @@ def compare_netcdf(base, gen, projection=False, update=None):
 
     # coordinates
     for coordname, da in xrb.coords.items():
-        assert coordname in xrg.coords
-        assert np.allclose(xrb.coords[coordname].data, xrg.coords[coordname].data)
-        for a in da.attrs:
-            if a == "grid_mapping" and not projection:
-                continue
-            assert a in xrg.coords[coordname].attrs
-            assert da.attrs[a] == xrg.coords[coordname].attrs[a]
+        compare_netcdf_var(
+            coordname, xrb.coords, xrg.data_vars, xrg.coords, projection, update
+        )
 
     # variables
     for varname, da in xrb.data_vars.items():
@@ -147,7 +143,7 @@ def test_load_gwfsto01(function_tmpdir, example_data_path):
 
         # load example
         sim = flopy.mf6.MFSimulation.load(sim_ws=base_path)
-        gwf = sim.get_model("gwf_sto01")
+        # gwf = sim.get_model("gwf_sto01")
 
         # set simulation path and write simulation
         sim.set_sim_path(test_path)
@@ -200,8 +196,7 @@ def test_update_gwfsto01(function_tmpdir, example_data_path):
 
     nlay, nrow, ncol = 3, 10, 10
 
-    dis_delr = np.array([1010, 1010, 1010, 1010, 1010, 1010, 1010, 1010, 1010, 1010])
-    dis_delc = [2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090]
+    dis_top = np.full((nrow, ncol), 50)
 
     # ic
     strt1 = np.full((nrow, ncol), 0.15)
@@ -210,8 +205,7 @@ def test_update_gwfsto01(function_tmpdir, example_data_path):
     ic_strt = np.array([strt1, strt2, strt3])
 
     update = {
-        "dis_delr": dis_delr,
-        "dis_delc": dis_delc,
+        "dis_top": dis_top.flatten()[0],
         "ic_strt": ic_strt,
         "ic_strt_l1": ic_strt[0].flatten(),
         "ic_strt_l2": ic_strt[1].flatten(),
@@ -233,9 +227,8 @@ def test_update_gwfsto01(function_tmpdir, example_data_path):
         # get model instance
         gwf = sim.get_model("gwf_sto01")
 
-        # update dis delr and delc
-        gwf.dis.delr = dis_delr
-        gwf.dis.delc = dis_delc
+        # update dis top
+        gwf.dis.top = dis_top
 
         # update ic strt
         gwf.ic.strt.set_data(ic_strt)
@@ -1090,7 +1083,9 @@ def test_dis_transform(function_tmpdir, example_data_path):
         ws = transform_ws / t
         sim.set_sim_path(ws)
         sim.write_simulation(netcdf=t)
-        compare_netcdf_data(cmp_pth / f"transform.{t}.nc", ws / "transform.in.nc")
+        compare_netcdf(
+            cmp_pth / f"transform.{t}.nc", ws / "transform.in.nc", projection=True
+        )
 
 
 @requires_exe("triangle")
@@ -1184,8 +1179,7 @@ def test_disv_transform(function_tmpdir, example_data_path):
 
     sim.write_simulation(netcdf=nc_type)
 
-    compare_netcdf_data(cmp_pth / f"tri.{nc_type}.nc", mf6_ws / "tri.in.nc")
-    # compare_netcdf_data(cmp_pth / f"tri.{nc_type}.nc", vertex_ws / "tri.nc")
+    compare_netcdf(cmp_pth / f"tri.{nc_type}.nc", mf6_ws / "tri.in.nc", projection=True)
 
 
 @pytest.mark.regression

examples/data/mf6/netcdf/test_dis_transform/transform.structured.nc

@@ -1839,7 +1839,6 @@ def store_external(
                     )
 
                     if self.netcdf:
-                        # TODO: multiplier NETCDF-DEV
                         file_access.set_netcdf_array(
                             self.layer_storage[layer].nc_dataset,
                             self.data_dimensions.structure.get_package(),
@@ -3018,7 +3017,10 @@ def _layer_prep(self, layer):
         if layer is None:
             # layer is none means the data provided is for all layers or this
             # is not layered data
-            layer = (0,)
+            if self.netcdf:
+                layer = (-1,)
+            else:
+                layer = (0,)
             self.layer_storage.list_shape = (1,)
             self.layer_storage.multi_dim_list = [
                 self.layer_storage.first_item()

examples/data/mf6/netcdf/test_gwf_sto01_structured/gwf_sto01/gwf_sto01.in.nc

@@ -930,6 +930,7 @@ def __init__(self):
         self.parameter_name = None
         self.one_per_pkg = False
         self.jagged_array = None
+        self.netcdf = False
 
     def set_value(self, line, common):
         arr_line = line.strip().split()
@@ -1107,6 +1108,8 @@ def set_value(self, line, common):
                 self.one_per_pkg = bool(arr_line[1])
             elif arr_line[0] == "jagged_array":
                 self.jagged_array = arr_line[1]
+            elif arr_line[0] == "netcdf":
+                self.netcdf = arr_line[1]
 
     def get_type_string(self):
         return f"[{self.type_string}]"
@@ -1425,14 +1428,7 @@ def __init__(self, data_item, model_data, package_type, dfn_list):
             or "nodes" in data_item.shape
             or len(data_item.layer_dims) > 1
         )
-        self.netcdf = (
-            ("ncol" in data_item.shape
-            or "nrow" in data_item.shape
-            or "ncpl" in data_item.shape
-            or "nlay" in data_item.shape
-            or "nodes" in data_item.shape)
-            and data_item.block_name == 'griddata'
-        )
+        self.netcdf = data_item.netcdf
         self.num_data_items = len(data_item.data_items)
         self.record_within_record = False
         self.file_data = False

flopy/mf6/data/mfdatastorage.py

@@ -1348,7 +1348,11 @@ def write(
         """
 
         # write netcdf file
-        if netcdf or self._nc_dataset is not None:
+        if (netcdf or self._nc_dataset is not None) and (
+            self.model_type == "gwf6"
+            or self.model_type == "gwt6"
+            or self.model_type == "gwe6"
+        ):
             kwargs = {}
             if self._nc_dataset is None:
                 from ..utils.model_netcdf import create_dataset
@@ -1361,7 +1365,7 @@ def write(
 
                 # create netcdf dataset
                 self._nc_dataset = create_dataset(
-                    self.model_type, self.name, netcdf, nc_fname, self._modelgrid
+                    self.model_type, self.name, netcdf, nc_fname, self.modelgrid
                 )
 
                 # reset data storage and populate netcdf file
@@ -1375,6 +1379,7 @@ def write(
                         kwargs["chunk_y"] = pp.chunk_y.get_data()
                         kwargs["chunk_z"] = pp.chunk_z.get_data()
                         kwargs["wkt"] = pp.wkt.get_data()
+
                     pp._set_netcdf_storage(
                         self._nc_dataset, create=True
                     )

flopy/mf6/data/mfstructure.py

@@ -19,6 +19,7 @@
 
 FILLNA_INT32 = np.int32(-2147483647)
 FILLNA_DBL = 9.96920996838687e36
+lenunits = {0: "m", 1: "ft", 2: "m", 3: "m"}
 
 
 class ModelNetCDFDataset:
@@ -45,7 +46,6 @@ class ModelNetCDFDataset:
     """
 
     def __init__(self):
-        # TODO: grid offsets
         self._modelname = None
         self._modeltype = None
         self._dataset = None
@@ -157,6 +157,9 @@ def set_array(
     ):
         """
         Set data in an existing array. Override this function in a derived class.
+        Do not use to update variables that establish the vertices and cell
+        structure of the grid, e.g. delr/delc or vertex/cell parameters, only
+        model data associated with the grid.
 
         Args:
             package (str): A package name provided as an argument to create_array().
@@ -286,64 +289,59 @@ def _set_encoding(self, **kwargs):
         return encoding
 
     def _set_projection(self):
-        if self._dis:
-            crs = None
-            wkt = None
-            if self._dis.crs:
-                try:
-                    crs = CRS.from_user_input(self._dis.crs)
-                    wkt = crs.to_wkt()
-                except Exception as e:
-                    # crs = None
-                    # wkt = None
-                    warnings.warn(
-                        f"Cannot generate CRS object from user input: {e}",
-                        UserWarning,
-                    )
-
-                if wkt is not None:
-                    # add projection variable
-                    self._dataset = self._dataset.assign(
-                        {"projection": ([], np.int32(1))}
-                    )
-                    if self._nc_type == "structured":
-                        self._dataset["projection"].attrs["crs_wkt"] = wkt
-                        self._dataset["x"].attrs["grid_mapping"] = "projection"
-                        self._dataset["y"].attrs["grid_mapping"] = "projection"
-                    elif self._nc_type == "mesh2d":
-                        self._dataset["projection"].attrs["wkt"] = wkt
-                        self._dataset["mesh_node_x"].attrs["grid_mapping"] = (
-                            "projection"
-                        )
-                        self._dataset["mesh_node_y"].attrs["grid_mapping"] = (
-                            "projection"
-                        )
-                        self._dataset["mesh_face_x"].attrs["grid_mapping"] = (
-                            "projection"
-                        )
-                        self._dataset["mesh_face_y"].attrs["grid_mapping"] = (
+        if not self._dis:
+            return
+
+        crs = None
+        wkt = None
+        projection = False
+        if self._dis.crs:
+            try:
+                crs = CRS.from_user_input(self._dis.crs)
+                wkt = crs.to_wkt()
+                projection = True
+            except Exception as e:
+                warnings.warn(
+                    f"Cannot generate CRS object from user input: {e}",
+                    UserWarning,
+                )
+
+        # update coords based on crs
+        coords = self._set_coords(crs)
+
+        # Don't define projection and grid mapping if using
+        # geographic coordinates in the structured type
+        if self._nc_type == "structured" and coords == "lon lat":
+            projection = False
+
+        if projection:
+            # add projection variable
+            self._dataset = self._dataset.assign({"projection": ([], np.int32(1))})
+            if self._nc_type == "structured":
+                self._dataset["projection"].attrs["crs_wkt"] = wkt
+                self._dataset["x"].attrs["grid_mapping"] = "projection"
+                self._dataset["y"].attrs["grid_mapping"] = "projection"
+            elif self._nc_type == "mesh2d":
+                self._dataset["projection"].attrs["wkt"] = wkt
+                self._dataset["mesh_node_x"].attrs["grid_mapping"] = "projection"
+                self._dataset["mesh_node_y"].attrs["grid_mapping"] = "projection"
+                self._dataset["mesh_face_x"].attrs["grid_mapping"] = "projection"
+                self._dataset["mesh_face_y"].attrs["grid_mapping"] = "projection"
+
+        # set grid_mapping and coordinates attributes
+        for p in self._tags:
+            for l in self._tags[p]:
+                dims = self._dataset[self._tags[p][l]].dims
+                if (self._nc_type == "structured" and len(dims) > 1) or (
+                    self._nc_type == "mesh2d"
+                    and ("nmesh_face" in dims or "nmesh_node" in dims)
+                ):
+                    if projection:
+                        self._dataset[self._tags[p][l]].attrs["grid_mapping"] = (
                             "projection"
                         )
-
-            # update coords based on crs
-            coords = self._set_coords(crs)
-
-            # set grid_mapping and coordinates attributes
-            for p in self._tags:
-                for l in self._tags[p]:
-                    dims = self._dataset[self._tags[p][l]].dims
-                    if (self._nc_type == "structured" and len(dims) > 1) or (
-                        self._nc_type == "mesh2d"
-                        and ("nmesh_face" in dims or "nmesh_node" in dims)
-                    ):
-                        if crs is not None:
-                            self._dataset[self._tags[p][l]].attrs["grid_mapping"] = (
-                                "projection"
-                            )
-                        if coords is not None:
-                            self._dataset[self._tags[p][l]].attrs["coordinates"] = (
-                                coords
-                            )
+                    if coords is not None:
+                        self._dataset[self._tags[p][l]].attrs["coordinates"] = coords
 
     def _set_modflow_attrs(self):
         if self._modeltype.endswith("6"):
@@ -545,10 +543,16 @@ def array(self, package: str, param: str, layer=None):
             return self._dataset[self._tags[path][-1]].data
 
     def _set_grid(self, dis):
-        lenunits = {0: "m", 1: "feet", 2: "m", 3: "cm"}
+        if dis.angrot != 0.0:
+            xoff = 0.0
+            yoff = 0.0
+        else:
+            xoff = dis.xoffset
+            yoff = dis.yoffset
 
+        # set coordinate var bounds
         x_bnds = []
-        xv = dis.xoffset + dis.xyedges[0]
+        xv = xoff + dis.xyedges[0]
         for idx, val in enumerate(xv):
             if idx + 1 < len(xv):
                 bnd = []
@@ -557,18 +561,17 @@ def _set_grid(self, dis):
                 x_bnds.append(bnd)
 
         y_bnds = []
-        yv = dis.yoffset + dis.xyedges[1]
+        yv = yoff + dis.xyedges[1]
         for idx, val in enumerate(yv):
             if idx + 1 < len(yv):
                 bnd = []
                 bnd.append(yv[idx + 1])
                 bnd.append(yv[idx])
                 y_bnds.append(bnd)
 
-        # x = dis.xcellcenters[0]
-        # y = dis.ycellcenters[:, 0]
-        x = dis.xoffset + dis.xycenters[0]
-        y = dis.yoffset + dis.xycenters[1]
+        # set coordinate vars
+        x = xoff + dis.xycenters[0]
+        y = yoff + dis.xycenters[1]
         z = [float(x) for x in range(1, dis.nlay + 1)]
 
         # create coordinate vars
@@ -715,8 +718,6 @@ def array(self, package: str, param: str, layer=None):
         return None
 
     def _set_grid(self, dis):
-        lenunits = {0: "m", 1: "feet", 2: "m", 3: "cm"}
-
         # mesh container variable
         self._dataset = self._dataset.assign({"mesh": ([], np.int32(1))})
         self._dataset["mesh"].attrs["cf_role"] = "mesh_topology"
@@ -770,7 +771,6 @@ def _set_grid(self, dis):
                     y_bnds.append(bnd)
 
         var_d = {
-            # TODO modflow6 and flopy results differ for mesh_face_x gwf_sto01
             "mesh_face_x": (["nmesh_face"], dis.xcellcenters.flatten()),
             "mesh_face_xbnds": (["nmesh_face", "max_nmesh_face_nodes"], x_bnds),
             "mesh_face_y": (["nmesh_face"], dis.ycellcenters.flatten()),
@@ -861,9 +861,6 @@ def array(self, package: str, param: str, layer=None):
         return None
 
     def _set_grid(self, disv):
-        # default metric "m" when undefined
-        lenunits = {0: "m", 1: "feet", 2: "m", 3: "cm"}
-
         # mesh container variable
         self._dataset = self._dataset.assign({"mesh": ([], np.int32(1))})
         self._dataset["mesh"].attrs["cf_role"] = "mesh_topology"