diff --git a/mfsetup/grid.py b/mfsetup/grid.py
index f28f9fcb..487e5101 100644
--- a/mfsetup/grid.py
+++ b/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))
diff --git a/mfsetup/tests/data/shellmound.yml b/mfsetup/tests/data/shellmound.yml
index 3917a438..6ebae823 100644
--- a/mfsetup/tests/data/shellmound.yml
+++ b/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'
diff --git a/mfsetup/tests/test_mf6_shellmound.py b/mfsetup/tests/test_mf6_shellmound.py
index 6391d32d..a3095cea 100644
--- a/mfsetup/tests/test_mf6_shellmound.py
+++ b/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':