fix(model_splitter.py): add trap for obs packages in bc packages (#2493)

* fix(model_splitter.py): add trap for obs packages in standard bc packages

* fix indexing error in observation remapping routine

* update testing syntax for ruff

* change testing syntax for ruff

* update test_package_observations

Author: jlarsen-usgs

autotest/test_model_splitter.py

@@ -1489,3 +1489,112 @@ def test_timeseries(function_tmpdir):
 
     if not success:
         raise AssertionError("Timeseries split simulation did not properly run")
+
+
+def test_package_observations():
+    XL = 10000
+    YL = 10500
+    dx = dy = 500
+    nlay = 3
+    nrow = int(YL / dy)
+    ncol = int(XL / dx)
+
+    delc = np.full((nrow,), dy)
+    delr = np.full((ncol,), dx)
+
+    top = np.full((nrow, ncol), 330)
+
+    botm = np.zeros((nlay, nrow, ncol))
+    botm[0] = 220
+    botm[1] = 200
+    botm[2] = 0
+
+    sim = flopy.mf6.MFSimulation()
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="ps1a_mod")
+
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+        idomain=np.ones(botm.shape, dtype=int),
+    )
+
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=330)
+
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf, k=[50, 0.01, 200], k33=[10, 0.01, 20], icelltype=[1, 0, 0]
+    )
+
+    # build the CHD stress period data
+    records = []
+    for row in range(nrow):
+        rec = ((0, row, 0), 330)
+        records.append(rec)
+
+    for row in range(nrow):
+        rec = ((0, row, ncol - 1), 320)
+        records.append(rec)
+
+    spd = {0: records}
+
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=spd)
+
+    observations = {
+        "chd_obs_out.csv": [
+            ("C1", "chd", (0, 4, 0)),  # 80 (end up as same)
+            ("C2", "chd", (0, 16, 0)),  # 320 (end up as m2 120)
+            ("R1", "chd", (0, 3, 19)),  # 79 (end up as same)
+            ("R2", "chd", (0, 20, 19)),  # 419 (end up as m2 219)
+        ]
+    }
+
+    # do observations on the chd package!!!!
+    obs = flopy.mf6.ModflowUtlobs(
+        chd,
+        continuous=observations,
+    )
+
+    # now do model splitter???
+    array = np.ones((nrow, ncol), dtype=int)
+    split_loc = ncol // 2
+    array[split_loc:, :] = 2
+
+    mfs = flopy.mf6.utils.Mf6Splitter(sim)
+    new_sim = mfs.split_model(array)
+
+    gwf1 = new_sim.get_model("ps1a_mod_1")
+    gwf2 = new_sim.get_model("ps1a_mod_2")
+
+    obsdata1 = gwf1.obs.continuous.data["chd_obs_out_1.csv"]
+    obsdata2 = gwf2.obs.continuous.data["chd_obs_out_2.csv"]
+
+    vdict1 = {"C1": (0, 4, 0), "R1": (0, 3, 19)}
+    vdict2 = {"C2": (0, 6, 0), "R2": (0, 10, 19)}
+
+    for row in obsdata1:
+        obsname = row.obsname
+        cellid = row.id
+        if obsname not in vdict1:
+            raise AssertionError("Observations not correctly mapped to new model")
+
+        vcid = vdict1[obsname]
+        if vcid != cellid:
+            raise AssertionError("Observation cellid not correctly mapped to new model")
+
+    for row in obsdata2:
+        obsname = row.obsname
+        cellid = row.id
+        if obsname not in vdict2:
+            raise AssertionError("Observations not correctly mapped to new model")
+
+        vcid = vdict2[obsname]
+        if vcid != cellid:
+            raise AssertionError("Observation cellid not correctly mapped to new model")

flopy/mf6/utils/model_splitter.py

@@ -2735,11 +2735,11 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None):
                             dtype=object,
                         )
                         for mkey, model in self._model_dict.items():
-                            idx = np.asarray(new_model1 == mkey).nonzero()
+                            idx = np.asarray(new_model1 == mkey).nonzero()[0]
                             idx = [
                                 ix
-                                for ix, i in enumerate(recarray.id[idx])
-                                if not isinstance(i, str)
+                                for ix in idx
+                                if not isinstance(recarray.id[ix], str)
                             ]
                             tmp_cellid = self._new_node_to_cellid(
                                 model, new_node1, layers1, idx
@@ -3382,9 +3382,28 @@ def _remap_package(self, package, ismvr=False):
                             "interpolation_methodrecord": tspkg.interpolation_methodrecord.array["interpolation_method"][0]
                         }
                         mapped_data[mkey]["timeseries"] = tsdict
+
                 else:
                     pass
 
+            if hasattr(package, "obs"):
+                obs_map = {"cellid": self._node_map}
+                for mkey, mdict in mapped_data.items():
+                    if "stress_period_data" in mdict:
+                        for _, ra in mdict["stress_period_data"].items():
+                            if isinstance(ra, dict):
+                                continue
+                            obs_map = self._set_boundname_remaps(
+                                ra, obs_map, list(obs_map.keys()), mkey
+                            )
+                for obspak in package.obs._packages:
+                    mapped_data = self._remap_obs(
+                        obspak,
+                        mapped_data,
+                        obs_map["cellid"],
+                        pkg_type=package.package_type,
+                    )
+
         observations = mapped_data.pop("observations", None)
 
         if "options" in package.blocks: