Feat exclude feature ids

docs/source/config-file.rst

@@ -143,8 +143,25 @@ Some additional notes on YAML
     * ``[12,1.2]`` parses to ``[12,1.2]``
     * ``(12,1.2)`` parses to ``"(12,1.2)"``
     * ``{12,1.2}`` parses to ``{12: None, 1.2: None}``
+* dictionaries can be represented with indentation, but spaces are needed after the colon(s):
 
+  .. code-block:: yaml
 
+      items:
+        0:1
+        1:2
+
+  parses to ``'0:1 1:2'``
+
+  .. code-block:: yaml
+
+      items:
+        0: 1
+        1: 2
+
+  parses to ``{0:1, 1:2}``
+
+Using a YAML-aware text editor such as VS Code can help with these issues, for example by changing the highlighting color to indicate a string in the first dictionary example above and an interpreted python data structure in the second dictionary example.
 
 .. _JSON: https://www.json.org/json-en.html
 .. _pyyaml: https://pyyaml.org/

mfsetup/grid.py

@@ -634,7 +634,7 @@ def write_bbox_shapefile(modelgrid, outshp):
 
 
 def rasterize(feature, grid, id_column=None,
-              include_ids=None, names_column=None,
+              include_ids=None, exclude_ids=None, names_column=None,
               crs=None, **kwargs):
     """Rasterize a feature onto the model grid, using
     the rasterio.features.rasterize method. Features are intersected
@@ -649,6 +649,8 @@ def rasterize(feature, grid, id_column=None,
         from this column will be assigned to the output raster.
     include_ids : sequence
         Subset of IDs in id_column to include
+    exclude_ids : sequence
+        Subset of IDs in id_column to exclude
     names_column : str, optional
         By default, the IDs in id_column, or sequential integers
         are returned. This option allows another column of strings
@@ -738,7 +740,8 @@ def rasterize(feature, grid, id_column=None,
     # subset to include_ids
     if id_column is not None and include_ids is not None:
         df = df.loc[df[id_column].isin(include_ids)].copy()
-
+    if id_column is not None and exclude_ids is not None:
+        df = df.loc[~df[id_column].isin(exclude_ids)].copy()
     # create list of GeoJSON features, with unique value for each feature
     if id_column is None:
         numbers = list(range(1, len(df)+1))
@@ -986,14 +989,14 @@ def setup_structured_grid(xoff=None, yoff=None, xul=None, yul=None,
     if parent_model is not None and snap_to_parent:
         to_grid_units_parent = convert_length_units(get_model_length_units(parent_model), grid_units)
         # parent model grid spacing in meters
-        parent_delr_grid = np.round(parent_model.dis.delr.array[0] * to_grid_units_parent, 4)
-        if not parent_delr_grid % delr_grid % parent_delr_grid == 0:
-            raise ValueError('inset delr spacing of {} must be factor of parent spacing of {}'.format(delr_grid,
-                                                                                                      parent_delr_grid))
-        parent_delc_grid = np.round(parent_model.dis.delc.array[0] * to_grid_units_parent, 4)
-        if not parent_delc_grid % delc_grid % parent_delc_grid == 0:
-            raise ValueError('inset delc spacing of {} must be factor of parent spacing of {}'.format(delc_grid,
-                                                                                                      parent_delc_grid))
+        #parent_delr_grid = np.round(parent_model.dis.delr.array[0] * to_grid_units_parent, 4)
+        #if not parent_delr_grid % delr_grid % parent_delr_grid == 0:
+        #    raise ValueError('inset delr spacing of {} must be factor of parent spacing of {}'.format(delr_grid,
+        #                                                                                              parent_delr_grid))
+        #parent_delc_grid = np.round(parent_model.dis.delc.array[0] * to_grid_units_parent, 4)
+        #if not parent_delc_grid % delc_grid % parent_delc_grid == 0:
+        #    raise ValueError('inset delc spacing of {} must be factor of parent spacing of {}'.format(delc_grid,
+        #                                                                                              parent_delc_grid))
 
     # option 1: make grid from xoff, yoff and specified dimensions
     if xoff is not None and yoff is not None:
@@ -1057,6 +1060,10 @@ def setup_structured_grid(xoff=None, yoff=None, xul=None, yul=None,
                 df.to_crs(crs, inplace=True)
             # use all features by default
             features = df.geometry.tolist()
+        elif features is None and features_shapefile is not None:
+            raise ValueError(
+                "setup_grid: need one of xoff/yoff, xul/yul, features_shapefile or "
+                "features inputs")
         # alternatively, accept features as an argument
         # convert multiple features to a MultiPolygon
         if isinstance(features, list):

mfsetup/mf6model.py

@@ -405,7 +405,7 @@ def create_lgr_models(self):
             last_refined_layer = max(np.where(is_refined > 0)[0])
             consecutive = all(np.diff(is_refined)[:last_refined_layer] == 0)
             if (is_refined[0] != 1) | (not consecutive):
-                raise ValueError("Configuration input: layer_refinement must"
+                raise ValueError("Configuration input: layer_refinement must "
                                  "include consecutive sequence of layers, "
                                  "starting with the top layer.")
             # check the specified DIS package input is consistent

mfsetup/mfmodel.py

@@ -1417,6 +1417,12 @@ def setup_grid(self):
             self._modelgrid.cfg = self.cfg['grid']
         else:
             kwargs = get_input_arguments(cfg, setup_structured_grid)
+            if not set(kwargs.keys()).intersection({
+                'features_shapefile', 'features', 'xoff', 'yoff', 'xul', 'yul'}):
+                raise ValueError(
+                    "No features_shapefile or xoff, yoff supplied "
+                    "to setup_grid: block. Check configuration file input, "
+                    "including for accidental indentation of the setup_grid: block.")
             self._modelgrid = setup_structured_grid(**kwargs)
             self.cfg['grid'] = self._modelgrid.cfg
             # update DIS package configuration

mfsetup/tests/data/shellmound.yml

@@ -317,10 +317,10 @@ drn:
     shapefile:
       filename: '../../../mfsetup/tests/data/shellmound/shps/waterbodies.shp'
       id_column: 'COMID'
-      # Note: to include all features in a shapefile,
-      # simply omit the include_ids: key
-      # id_column: can also be omitted
-      include_ids: [18047154, 18046236]
+      # Include all features in the above shapefile,
+      # except those associated with these COMIDs
+      exclude_ids: [18046230, 18046226, 18046238, 17953939, 18046140, 18046162]
+      boundname_column: 'COMID'
     elev:
       filename: 'shellmound/rasters/meras_100m_dem.tif'
       elevation_units: 'feet'

mfsetup/tests/test_mf6_shellmound.py

@@ -650,6 +650,11 @@ def test_basic_stress_package_setup(shellmound_model_with_dis, pckg_abbrv,
         # heads in CHD package should match the CSV input,
         # after conversion to meters
         assert np.allclose(in_df['head'].values, spd_heads[1:]/.3048)
+    if pckg_abbrv == 'drn':
+        exclude_ids = list(map(str, m.cfg[pckg_abbrv]['source_data']['shapefile']['exclude_ids']))
+        comids = [s.replace('feature-','') for s in np.unique(m.drn.stress_period_data.data[0]['boundname'])]
+        assert not set(exclude_ids).intersection(comids)
+        assert comids == ['18046236', '18047154']
     # more advanced transient input case with csv
     # and conductance values specified via a raster
     if pckg_abbrv == 'ghb':