Merge pull request #1510 from phargogh/bugfix/1503-admin-boundary-ove…

…rlaps-noone

UNA: Handle the case where an admin geometry doesn't overlap any people

HISTORY.rst

@@ -42,6 +42,9 @@ Unreleased Changes
     * Fixed a ``NameError`` that occurred when running the model using
       search radii defined per population group with an exponential search
       kernel. https://github.com/natcap/invest/issues/1502
+    * Fixed an issue where Urban Nature Access would crash if an administrative
+      boundary geometry did not overlap any people in the population raster.
+      https://github.com/natcap/invest/issues/1503
 
 3.14.1 (2023-12-18)
 -------------------

src/natcap/invest/urban_nature_access.py

@@ -347,9 +347,12 @@
                             "type": "number",
                             "units": u.m**2/u.person,
                             "about": (
-                                "The averge urban nature supply/demand "
+                                "The average urban nature supply/demand "
                                 "balance available per person within this "
-                                "administrative unit.")
+                                "administrative unit. If no people reside "
+                                "within this administrative unit, this field "
+                                "will have no value (NaN, NULL or None, "
+                                "depending on your GIS software).")
                         },
                         "Pund_adm": {
                             "type": "number",
@@ -2050,12 +2053,19 @@ def _get_zonal_stats(raster_path):
     stats_by_feature = {}
     for fid in urban_nature_budget_stats.keys():
         stats = {
-            'SUP_DEMadm_cap': (
-                urban_nature_budget_stats[fid]['sum'] /
-                population_stats[fid]['sum']),
             'Pund_adm': undersupplied_stats[fid]['sum'],
             'Povr_adm': oversupplied_stats[fid]['sum'],
         }
+
+        # Handle the case where an administrative unit might overlap no people
+        if population_stats[fid]['sum'] == 0:
+            per_capita_supply = float('nan')
+        else:
+            per_capita_supply = (
+                urban_nature_budget_stats[fid]['sum'] /
+                population_stats[fid]['sum'])
+        stats['SUP_DEMadm_cap'] = per_capita_supply
+
         for pop_group_field in pop_group_fields:
             group = group_names[pop_group_field]
             group_proportion = pop_proportions_by_fid[fid][group]

tests/test_urban_nature_access.py

@@ -915,6 +915,48 @@ def test_write_vector(self):
             layer = None
             vector = None
 
+    def test_write_vector_for_single_raster_modes(self):
+        """UNA: create a summary vector for single-raster summary stats."""
+
+        from natcap.invest import urban_nature_access
+
+        args = _build_model_args(self.workspace_dir)
+
+        # Overwrite all population pixels with 0
+        try:
+            raster = gdal.Open(args['population_raster_path'], gdal.GA_Update)
+            band = raster.GetRasterBand(1)
+            array = band.ReadAsArray()
+            array.fill(0.0)
+            band.WriteArray(array)
+        finally:
+            raster = band = None
+
+        args['search_radius_mode'] = urban_nature_access.RADIUS_OPT_UNIFORM
+        args['search_radius'] = 1000
+
+        urban_nature_access.execute(args)
+
+        summary_vector = gdal.OpenEx(
+            os.path.join(args['workspace_dir'], 'output',
+                         'admin_boundaries_suffix.gpkg'))
+        summary_layer = summary_vector.GetLayer()
+        self.assertEqual(summary_layer.GetFeatureCount(), 1)
+        summary_feature = summary_layer.GetFeature(1)
+
+        expected_field_values = {
+            'adm_unit_id': 0,
+            'Pund_adm': 0,
+            'Povr_adm': 0,
+            'SUP_DEMadm_cap': None,  # OGR converts NaN to None.
+        }
+        self.assertEqual(
+            set(defn.GetName() for defn in summary_layer.schema),
+            set(expected_field_values.keys()))
+        for fieldname, expected_value in expected_field_values.items():
+            self.assertAlmostEqual(
+                expected_value, summary_feature.GetField(fieldname))
+
     def test_urban_nature_proportion(self):
         """UNA: Run the model with urban nature proportion."""
         from natcap.invest import urban_nature_access