move data update to package

mjreno · mjreno · commit 2747223f04f2 · 2025-09-03T08:13:23.000-04:00
diff --git a/.docs/Notebooks/mf6_netcdf01_tutorial.py b/.docs/Notebooks/mf6_netcdf01_tutorial.py
@@ -127,6 +127,12 @@
 # write simulation with structured NetCDF file
 sim.write_simulation(netcdf="structured")
 
+# ## Run MODFLOW 6 simulation with NetCDF input
+#
+# The simulation generated by this tutorial should be runnable by
+# Extended MODFLOW 6, available from the nightly-build repository
+# (linked above).
+
 # success, buff = sim.run_simulation(silent=True, report=True)
 # assert success, pformat(buff)
 
@@ -140,6 +146,12 @@
 # write simulation with with layered mesh NetCDF file
 sim.write_simulation(netcdf="layered")
 
+# ## Run MODFLOW 6 simulation with NetCDF input
+#
+# The simulation generated by this tutorial should be runnable by
+# Extended MODFLOW 6, available from the nightly-build repository
+# (linked above).
+
 # success, buff = sim.run_simulation(silent=True, report=True)
 # assert success, pformat(buff)
 
@@ -225,23 +237,6 @@
 dis = gwf.get_package("dis")
 ds = dis.update_dataset(ds)
 
-# ## Update array data
-#
-# We have created dataset array variables for the package but they do not yet
-# define the expected input data for MODFLOW 6. We will take advantage of the
-# existing simulation objects and update the dataset.
-#
-# Default dataset variable names are defined in the package `netcdf_info()`
-# dictionary. Here we will use the info dictionary to programmatically update
-# the dataset- for remaining packages we will hardcode the variable names
-# being updated for maximum clarity.
-
-nc_info = dis.netcdf_info()
-for v in nc_info:
-    name = nc_info[v]["attrs"]["modflow_input"].rsplit("/", 1)[1].lower()
-    d = getattr(dis, name)
-    ds[nc_info[v]["varname"]].values = d.get_data()
-
 # ## Access `NPF` package NetCDF attributes
 #
 # Access package scoped NetCDF details by storing the dictionary returned
@@ -270,12 +265,6 @@
 nc_info["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
 ds = npf.update_dataset(ds, netcdf_info=nc_info)
 
-# ## Update array data
-
-# update dataset from npf arrays
-ds["npf_icelltype"].values = npf.icelltype.get_data()
-ds["npf_k_updated"].values = npf.k.get_data()
-
 # ## Show dataset `NPF K` parameter with updates
 
 # print dataset npf k variable
@@ -286,37 +275,19 @@
 # ic
 ic = gwf.get_package("ic")
 ds = ic.update_dataset(ds)
-ds["ic_strt"].values = ic.strt.get_data()
 
 # ## Update the dataset with supported `STO` arrays
 
 # storage
 sto = gwf.get_package("sto")
 ds = sto.update_dataset(ds)
-ds["sto_iconvert"].values = sto.iconvert.get_data()
-ds["sto_sy"].values = sto.sy.get_data()
-ds["sto_ss"].values = sto.ss.get_data()
 
 # ## Update the dataset with supported `GHBG` arrays
 
 # update dataset with 'GHBG' arrays
 ghbg = gwf.get_package("ghbg_0")
 ds = ghbg.update_dataset(ds)
 
-# ## Update array data
-
-# update bhead netcdf array from flopy perioddata
-# timeseries step index is first of stress period
-for p in ghbg.bhead.get_data():
-    istp = sum(gwf.modeltime.nstp[0:p])
-    ds["ghbg_0_bhead"].values[istp] = ghbg.bhead.get_data()[p]
-
-# update cond netcdf array from flopy perioddata
-# timeseries step index is first of stress period
-for p in ghbg.cond.get_data():
-    istp = sum(gwf.modeltime.nstp[0:p])
-    ds["ghbg_0_cond"].values[istp] = ghbg.cond.get_data()[p]
-
 # ## Display generated dataset
 
 # show the dataset
@@ -337,3 +308,30 @@
 
 # success, buff = sim.run_simulation(silent=True, report=True)
 # assert success, pformat(buff)
+
+# ## Method 4: DIY with xarray
+#
+# The above method still uses FloPy objects to update the dataset arrays
+# to values consistent with the state of the objects. The `netcdf_info`
+# dictionary is intended to supported creation of the dataset without
+# an existing simulation defined. The base dataset can be defined with
+# `modelgrid` and `modeltime` objects, while the full package netcdf
+# dictionary can retrieved with a static call to a model or package mf6
+# type. The auxiliary input is optional but does show the variables that
+# would be required if auxiliary variables were defined in the package.
+
+# ## Demonstrate static call on MFPackage (structured dataset):
+
+netcdf_info = flopy.mf6.mfpackage.MFPackage.netcdf_package(
+    mtype="GWF",
+    ptype="NPF",
+    auxiliary=["CONCENTRATION"],
+)
+pprint(netcdf_info)
+
+# ## Demonstrate static call on MFPackage (layered dataset):
+
+netcdf_info = flopy.mf6.mfpackage.MFPackage.netcdf_package(
+    mtype="GWF", ptype="NPF", mesh="LAYERED", auxiliary=["CONCENTRATION"], nlay=2
+)
+pprint(netcdf_info)
diff --git a/autotest/regression/test_model_netcdf.py b/autotest/regression/test_model_netcdf.py
@@ -49,6 +49,8 @@ def check_netcdf(path, mobj, mesh=None):
         if p in packages:
             l = -1
             assert "modflow_input" in da.attrs
+            assert "longname" in da.attrs
+            assert da.attrs["longname"] != ""
             if mesh is None:
                 assert "layer" not in da.attrs
             else:
@@ -82,14 +84,6 @@ def check_netcdf(path, mobj, mesh=None):
                     assert np.allclose(ds.data_vars[varname].values, d.get_data())
 
 
-def update_dataset(dataset, pobj):
-    nc_info = pobj.netcdf_info()
-    for v in nc_info:
-        name = nc_info[v]["attrs"]["modflow_input"].rsplit("/", 1)[1].lower()
-        d = getattr(pobj, name)
-        dataset[nc_info[v]["varname"]].values = d.get_data()
-
-
 @pytest.mark.regression
 def test_uzf01_model_scope_nofile(function_tmpdir, example_data_path):
     sim_name = "uzf01"
@@ -350,16 +344,6 @@ def test_uzf01_pkg_scope(function_tmpdir, example_data_path):
     # add all packages and update data
     for p in gwf.packagelist:
         ds = p.update_dataset(ds)
-        nc_info = p.netcdf_info()
-        for v in nc_info:
-            name = nc_info[v]["attrs"]["modflow_input"].rsplit("/", 1)[1].lower()
-            d = getattr(p, name)
-            if d.repeating:
-                for per in d.get_data():
-                    istp = sum(gwf.modeltime.nstp[0:per])
-                    ds[nc_info[v]["varname"]].values[istp] = d.get_data()[per]
-            else:
-                ds[nc_info[v]["varname"]].values = d.get_data()
 
     # write dataset to netcdf
     ds.to_netcdf(ws / fname, format="NETCDF4", engine="netcdf4")
@@ -398,7 +382,6 @@ def test_uzf01_pkg_scope_modify(function_tmpdir, example_data_path):
     # update dataset with `DIS` arrays
     dis = gwf.get_package("dis")
     ds = dis.update_dataset(ds)
-    update_dataset(ds, dis)
 
     # get npf package netcdf info
     npf = gwf.get_package("npf")
@@ -410,36 +393,18 @@ def test_uzf01_pkg_scope_modify(function_tmpdir, example_data_path):
     nc_info["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
     ds = npf.update_dataset(ds, netcdf_info=nc_info)
 
-    # update dataset from npf arrays
-    ds["npf_icelltype"].values = npf.icelltype.get_data()
-    ds["npf_k_updated"].values = npf.k.get_data()
-
     # ic
     ic = gwf.get_package("ic")
     ds = ic.update_dataset(ds)
-    update_dataset(ds, ic)
 
     # storage
     sto = gwf.get_package("sto")
     ds = sto.update_dataset(ds)
-    update_dataset(ds, sto)
 
     # update dataset with 'GHBG' arrays
     ghbg = gwf.get_package("ghbg_0")
     ds = ghbg.update_dataset(ds)
 
-    # update bhead netcdf array from flopy perioddata
-    # timeseries step index is first of stress period
-    for p in ghbg.bhead.get_data():
-        istp = sum(gwf.modeltime.nstp[0:p])
-        ds["ghbg_0_bhead"].values[istp] = ghbg.bhead.get_data()[p]
-
-    # update cond netcdf array from flopy perioddata
-    # timeseries step index is first of stress period
-    for p in ghbg.cond.get_data():
-        istp = sum(gwf.modeltime.nstp[0:p])
-        ds["ghbg_0_cond"].values[istp] = ghbg.cond.get_data()[p]
-
     # write dataset to netcdf
     ds.to_netcdf(ws / fname, format="NETCDF4", engine="netcdf4")
 
diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py
@@ -1242,23 +1242,26 @@ def _get_file_entry(
                 layer_min = layer
                 layer_max = shape_ml.inc_shape_idx(layer)
 
-            if data_storage.netcdf:
-                file_entry_array.append(f"{indent}{self.structure.name}{indent}NETCDF\n")
-
-            elif layered_aux:
+            if layered_aux:
                 aux_var_names = (
                     self.data_dimensions.package_dim.get_aux_variables()[0]
                 )
                 for layer in range(0, len(aux_var_names)-1):
-                    file_entry_array.append(
-                        self._get_file_entry_layer(
-                            [layer],
-                            data_indent,
-                            data_storage.layer_storage[layer].data_storage_type,
-                            ext_file_action,
-                            layered_aux,
+                    if data_storage.netcdf:
+                        if data_storage.has_data():
+                            file_entry_array.append(f"{indent}{aux_var_names[layer+1]}{indent}NETCDF\n")
+                    else:
+                        file_entry_array.append(
+                            self._get_file_entry_layer(
+                                [layer],
+                                data_indent,
+                                data_storage.layer_storage[layer].data_storage_type,
+                                ext_file_action,
+                                layered_aux,
+                            )
                         )
-                    )
+            elif data_storage.netcdf:
+                file_entry_array.append(f"{indent}{self.structure.name}{indent}NETCDF\n")
 
             else:
                 for layer in shape_ml.indexes(layer_min, layer_max):
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
@@ -2371,8 +2371,15 @@ def _build_full_data(self, apply_multiplier=False):
         else:
             fill_value = None
         full_data = np.full(dimensions, fill_value, np_full_data_type)
-        is_aux = self.data_dimensions.structure.name == "aux"
-        if is_aux:
+        layer_aux = (
+            self.data_dimensions.structure.name == "aux"
+            and not self.layered
+        )
+        grid_aux = (
+            self.data_dimensions.structure.name == "aux"
+            and self.layered
+        )
+        if layer_aux:
             aux_data = []
         if not self.layered:
             layers_to_process = [0]
@@ -2398,7 +2405,7 @@ def _build_full_data(self, apply_multiplier=False):
                     or len(self.layer_storage[layer].internal_data) > 0
                     and self.layer_storage[layer].internal_data[0] is None
                 ):
-                    if is_aux:
+                    if layer_aux:
                         full_data = None
                     else:
                         return None
@@ -2409,9 +2416,14 @@ def _build_full_data(self, apply_multiplier=False):
                 ):
                     full_data = self.layer_storage[layer].internal_data * mult
                 else:
-                    full_data[layer] = (
-                        self.layer_storage[layer].internal_data * mult
-                    )
+                    if grid_aux:
+                        full_data = (
+                            self.layer_storage[layer].internal_data * mult
+                        )
+                    else:
+                        full_data[layer] = (
+                            self.layer_storage[layer].internal_data * mult
+                        )
             elif (
                 self.layer_storage[layer].data_storage_type
                 == DataStorageType.internal_constant
@@ -2469,11 +2481,11 @@ def _build_full_data(self, apply_multiplier=False):
                 ):
                     full_data = data_out
                 else:
-                    if is_aux and full_data.shape == data_out.shape:
+                    if layer_aux and full_data.shape == data_out.shape:
                         full_data = data_out
                     else:
                         full_data[layer] = data_out
-            if is_aux:
+            if layer_aux:
                 if full_data is not None:
                     all_none = False
                 aux_data.append(full_data)
@@ -2482,7 +2494,7 @@ def _build_full_data(self, apply_multiplier=False):
                     np.nan,
                     self.data_dimensions.structure.get_datum_type(True),
                 )
-        if is_aux:
+        if layer_aux:
             if all_none:
                 return None
             else:
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py
