Watershed delineation will now always create multipolygons

HISTORY.rst

@@ -1,7 +1,11 @@
 Release History
 ===============
 
-.. Unreleased Changes
+Unreleased Changes
+------------------
+* Watershed delineation will now always create a layer of type MultiPolygon.
+  All features contained in this layer will subsequently be MultiPolygons.
+  https://github.com/natcap/pygeoprocessing/issues/435
 * ``raster_calculator`` ``calc_raster_stats=True`` will count the number
   of non-nodata pixels in the target raster and write
   'STATISTICS_VALID_PERCENT' metadata along with other stats.

src/pygeoprocessing/routing/watershed.pyx

@@ -729,13 +729,12 @@ def delineate_watersheds_d8(
     polygonized_watersheds_layer = watersheds_vector.CreateLayer(
         'polygonized_watersheds', watersheds_srs, ogr.wkbPolygon)
 
-    # Using wkbUnknown layer data type because it's possible for a single
-    # watershed to have several disjoint components, and inserting a
-    # Polygon geometry into a MultiPolygon layer (or vice versa) is
-    # technically not supported by the GPKG standard although GDAL
-    # allows it for the time being.
+    # Using wkbMultiPolygon for this layer because GDAL's polygonize function
+    # may create multiple polygons for a single watershed feature, and we want
+    # all geometries for a single watershed to be represented in a single
+    # feature.
     watersheds_layer = watersheds_vector.CreateLayer(
-        target_layer_name, watersheds_srs, ogr.wkbUnknown)
+        target_layer_name, watersheds_srs, ogr.wkbMultiPolygon)
     index_field = ogr.FieldDefn('ws_id', ogr.OFTInteger)
     index_field.SetWidth(24)
     polygonized_watersheds_layer.CreateField(index_field)
@@ -1000,7 +999,7 @@ def delineate_watersheds_d8(
                 new_geometry.AddGeometry(duplicate_geometry)
 
         watershed_feature = ogr.Feature(watersheds_layer.GetLayerDefn())
-        watershed_feature.SetGeometry(new_geometry)
+        watershed_feature.SetGeometry(ogr.ForceToMultiPolygon(new_geometry))
 
         source_feature = source_layer.GetFeature(ws_id)
         for field_name, field_value in source_feature.items().items():
@@ -1009,11 +1008,11 @@ def delineate_watersheds_d8(
     watersheds_layer.CommitTransaction()
 
     polygonized_watersheds_layer = None
+    watersheds_layer = None
     if remove_temp_files:
         watersheds_vector.DeleteLayer('polygonized_watersheds')
     LOGGER.info('Finished vector consolidation')
 
-    watersheds_layer = None
     watersheds_vector = None
     source_layer = None
     source_vector = None

tests/test_watershed_delineation.py

@@ -177,6 +177,8 @@ def test_watersheds_trivial(self):
                 watersheds_layer = watersheds_vector.GetLayer(
                     'watersheds_something')
                 self.assertEqual(watersheds_layer.GetFeatureCount(), 4)
+                self.assertEqual(watersheds_layer.GetGeomType(),
+                                 ogr.wkbMultiPolygon)
 
                 # All features should have the same watersheds, both in area
                 # and geometry.
@@ -193,11 +195,13 @@ def test_watersheds_trivial(self):
                 pygeoprocessing.shapely_geometry_to_vector(
                     [expected_watershed_geometry],
                     os.path.join(self.workspace_dir, 'foo.gpkg'), srs_wkt,
-                    'GPKG', ogr_geom_type=ogr.wkbGeometryCollection)
+                    'GPKG', ogr_geom_type=ogr.wkbPolygon)
 
                 id_to_fields = {}
                 for feature in watersheds_layer:
                     geometry = feature.GetGeometryRef()
+                    self.assertEqual(geometry.GetGeometryType(),
+                                     ogr.wkbMultiPolygon)
                     shapely_geom = shapely.wkb.loads(
                         bytes(geometry.ExportToWkb()))
                     self.assertEqual(