Repository: natcap/invest (python)

Issue:
Urban Nature Access: handle the case where an admin boundary doesn't overlap any people
A [user on the forums](https://community.naturalcapitalproject.org/t/something-went-wrong-when-adding-task-aggregate-supply-demand-to-admin-units-single-rasters-17/4405) encountered a `ZeroDivisionError` that appears to be caused by an admin boundary polygon that doesn't overlap any people in the population raster.  The workaround suggested to the user was to locate that polygon and remove it from the vector until a fix can be implemented.

Pull request:
UNA: Handle the case where an admin geometry doesn't overlap any people
This PR addresses a case where an admin geometry doesn't overlap any people in the population raster, as might happen with a rural park.  Per the science team, `NaN` (which OGR saves as NULL, it appears) is used as a value in this case.

Fixes #1503

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
- [x] Updated the user's guide (if needed)
~~- [ ] Tested the Workbench UI (if relevant)~~


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index c5c79f61d..7a0341fc5 100644
--- a/HISTORY.rst
+++ b/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)
 -------------------
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index 05928c990..f856da5dd 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -347,9 +347,12 @@ MODEL_SPEC = {
                             "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 _supply_demand_vector_for_single_raster_modes(
     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]


Test patch:
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index 59e278b41..4de4c58b9 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -915,6 +915,48 @@ class UNATests(unittest.TestCase):
             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