Use ESRI SHP driver for reprojections to avoid off-by-one scenario wh…

…en mapping zonal stats to output vector

src/natcap/invest/urban_flood_risk_mitigation.py

@@ -173,14 +173,14 @@
                     "about": "Map of runoff volume.",
                     "bands": {1: {"type": "number", "units": u.meter**3}}
                 },
-                "reprojected_aoi.gpkg": {
+                "reprojected_aoi.shp": {
                     "about": (
                         "Copy of AOI vector reprojected to the same spatial "
                         "reference as the LULC."),
                     "geometries": spec_utils.POLYGONS,
                     "fields": {}
                 },
-                "structures_reprojected.gpkg": {
+                "structures_reprojected.shp": {
                     "about": (
                         "Copy of built infrastructure vector reprojected to "
                         "the same spatial reference as the LULC."),
@@ -408,14 +408,14 @@ def execute(args):
         task_name='calculate service built raster')
 
     reprojected_aoi_path = os.path.join(
-        intermediate_dir, 'reprojected_aoi.gpkg')
+        intermediate_dir, 'reprojected_aoi')
     reprojected_aoi_task = task_graph.add_task(
         func=pygeoprocessing.reproject_vector,
         args=(
             args['aoi_watersheds_path'],
             target_sr_wkt,
             reprojected_aoi_path),
-        kwargs={'driver_name': 'GPKG'},
+        kwargs={'driver_name': 'ESRI Shapefile'},
         target_path_list=[reprojected_aoi_path],
         task_name='reproject aoi/watersheds')
 
@@ -435,7 +435,7 @@ def execute(args):
             (runoff_retention_raster_path, 1),
             reprojected_aoi_path),
         store_result=True,
-        dependent_task_list=[runoff_retention_task],
+        dependent_task_list=[runoff_retention_task, reprojected_aoi_task],
         task_name='zonal_statistics over runoff_retention raster')
 
     runoff_retention_volume_stats_task = task_graph.add_task(
@@ -444,7 +444,7 @@ def execute(args):
             (runoff_retention_vol_raster_path, 1),
             reprojected_aoi_path),
         store_result=True,
-        dependent_task_list=[runoff_retention_vol_task],
+        dependent_task_list=[runoff_retention_vol_task, reprojected_aoi_task],
         task_name='zonal_statistics over runoff_retention_volume raster')
 
     damage_per_aoi_stats = None
@@ -457,13 +457,13 @@ def execute(args):
             args['built_infrastructure_vector_path'] not in ('', None)):
         # Reproject the built infrastructure vector to the target SRS.
         reprojected_structures_path = os.path.join(
-            intermediate_dir, 'structures_reprojected.gpkg')
+            intermediate_dir, 'structures_reprojected')
         reproject_built_infrastructure_task = task_graph.add_task(
             func=pygeoprocessing.reproject_vector,
             args=(args['built_infrastructure_vector_path'],
                   target_sr_wkt,
                   reprojected_structures_path),
-            kwargs={'driver_name': 'GPKG'},
+            kwargs={'driver_name': 'ESRI Shapefile'},
             target_path_list=[reprojected_structures_path],
             task_name='reproject built infrastructure to target SRS')
 
@@ -512,7 +512,7 @@ def _write_summary_vector(
         runoff_ret_vol_stats, flood_volume_stats, damage_per_aoi_stats=None):
     """Write a vector with summary statistics.
 
-    This vector will always contain two fields::
+    This vector will always contain three fields::
 
         * ``'flood_vol'``: The volume of flood (runoff), in m3, per watershed.
         * ``'rnf_rt_idx'``: Average of runoff retention values per watershed
@@ -571,34 +571,31 @@ def _write_summary_vector(
 
     target_watershed_layer.ResetReading()
     for target_feature in target_watershed_layer:
-        # Target vector is SHP, where FIDs start at 0, but stats were
-        # generated based on GPKG reprojection, where FIDs start at 1.
-        # Therefore, we need to reference stats at SHP FID + 1.
-        stat_key = target_feature.GetFID() + 1
+        feature_id = target_feature.GetFID()
 
-        pixel_count = runoff_ret_stats[stat_key]['count']
+        pixel_count = runoff_ret_stats[feature_id]['count']
         if pixel_count > 0:
             mean_value = (
-                runoff_ret_stats[stat_key]['sum'] / float(pixel_count))
+                runoff_ret_stats[feature_id]['sum'] / float(pixel_count))
             target_feature.SetField('rnf_rt_idx', float(mean_value))
 
         target_feature.SetField(
             'rnf_rt_m3', float(
-                runoff_ret_vol_stats[stat_key]['sum']))
+                runoff_ret_vol_stats[feature_id]['sum']))
 
         if damage_per_aoi_stats is not None:
-            pixel_count = runoff_ret_vol_stats[stat_key]['count']
+            pixel_count = runoff_ret_vol_stats[feature_id]['count']
             if pixel_count > 0:
-                damage_sum = damage_per_aoi_stats[stat_key]
+                damage_sum = damage_per_aoi_stats[feature_id]
                 target_feature.SetField('aff_bld', damage_sum)
 
                 # This is the service_built equation.
                 target_feature.SetField(
                     'serv_blt', (
-                        damage_sum * runoff_ret_vol_stats[stat_key]['sum']))
+                        damage_sum * runoff_ret_vol_stats[feature_id]['sum']))
 
         target_feature.SetField(
-            'flood_vol', float(flood_volume_stats[stat_key]['sum']))
+            'flood_vol', float(flood_volume_stats[feature_id]['sum']))
 
         target_watershed_layer.SetFeature(target_feature)