Repository: natcap/invest
Pull request #743 (2021-12-01 17:04:52)

Urban Nature Access - Greenspace Budget and Summarizing by Admin Units 
This PR implements the Greenspace Budgets part of UNA and also summarizes the results by the admin units vector.

I have an outstanding question to the model devs about whether we're likely to see admin unit geometries that are complex enough to need simplification, but since the sample data definitely doesn't need simplification I've left it out so far.


Patch:
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index e8a6f6618..457c85c1c 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -5,9 +5,11 @@ import shutil
 import tempfile
 
 import numpy
+import numpy.testing
 import pygeoprocessing
 import taskgraph
 from osgeo import gdal
+from osgeo import ogr
 from osgeo import osr
 
 from natcap.invest import MODEL_METADATA
@@ -91,7 +93,11 @@ ARGS_SPEC = {
             'name': 'greenspace demand per capita',
             'units': u.m**2,  # defined as m² per capita
             'expression': "value > 0",
-            'about': "",  # TODO, will know more about this when I implement.
+            'about': (
+                "The amount of greenspace that each resident should have "
+                "access to. This is often defined by local urban planning "
+                "documents."
+            )
         },
         'search_radius': {
             'type': 'number',
@@ -118,6 +124,7 @@ ARGS_SPEC = {
 
 _OUTPUT_BASE_FILES = {
     'greenspace_supply': 'greenspace_supply.tif',
+    'admin_units': 'admin_units.gpkg',
 }
 _INTERMEDIATE_BASE_FILES = {
     'aligned_population': 'aligned_population.tif',
@@ -126,6 +133,8 @@ _INTERMEDIATE_BASE_FILES = {
     'kernel': 'kernel.tif',
     'greenspace_population_ratio': 'greenspace_population_ratio.tif',
     'convolved_population': 'convolved_population.tif',
+    'greenspace_budget': 'greenspace_budget.tif',
+    'greenspace_supply_demand_budget': 'greenspace_supply_demand_budget.tif',
 }
 
 
@@ -329,6 +338,10 @@ def execute(args):
             'kernel_path_band': (file_registry['kernel'], 1),
             'target_path': file_registry['greenspace_supply'],
             'working_dir': intermediate_dir,
+            # Insurance against future pygeoprocessing API changes.  The target
+            # nodata right now is the minimum possible numpy float32 value,
+            # which is also what we use here as FLOAT32_NODATA.
+            'target_nodata': FLOAT32_NODATA,
         },
         task_name='2SFCA - greenspace supply',
         target_path_list=[file_registry['greenspace_supply']],
@@ -337,11 +350,179 @@ def execute(args):
             greenspace_population_ratio_task,
         ])
 
+    # This is "SUP_DEMi_cap" from the user's guide
+    per_capita_greenspace_budget_task = graph.add_task(
+        pygeoprocessing.raster_calculator,
+        kwargs={
+            'base_raster_path_band_const_list': [
+                (file_registry['greenspace_supply'], 1),
+                (float(args['greenspace_demand']), 'raw')
+            ],
+            'local_op': _greenspace_budget_op,
+            'target_raster_path': file_registry['greenspace_budget'],
+            'datatype_target': gdal.GDT_Float32,
+            'nodata_target': FLOAT32_NODATA
+        },
+        task_name='Calculate per-capita greenspace budget',
+        target_path_list=[file_registry['greenspace_budget']],
+        dependent_task_list=[
+            convolved_greenspace_population_ratio_task
+        ])
+
+    # This is "SUP_DEMi" from the user's guide
+    greenspace_supply_demand_task = graph.add_task(
+        pygeoprocessing.raster_calculator,
+        kwargs={
+            'base_raster_path_band_const_list': [
+                (file_registry['greenspace_budget'], 1),
+                (file_registry['aligned_population'], 1)
+            ],
+            'local_op': _greenspace_supply_demand_op,
+            'target_raster_path': (
+                file_registry['greenspace_supply_demand_budget']),
+            'datatype_target': gdal.GDT_Float32,
+            'nodata_target': FLOAT32_NODATA
+        },
+        task_name='Calculate per-capita greenspace supply-demand',
+        target_path_list=[file_registry['greenspace_supply_demand_budget']],
+        dependent_task_list=[
+             per_capita_greenspace_budget_task,
+             population_alignment_task,
+        ])
+
+    aggregate_admin_units_task = graph.add_task(
+        _admin_level_supply_demand,
+        kwargs={
+            'greenspace_budget_path': file_registry[
+                'greenspace_supply_demand_budget'],
+            'population_path': file_registry['aligned_population'],
+            'admin_unit_vector_path': args['admin_unit_vector_path'],
+            'target_admin_unit_vector_path': file_registry['admin_units'],
+        },
+        task_name='Aggregate supply-demand to the admin units',
+        target_path_list=[file_registry['admin_units']],
+        dependent_task_list=[
+            greenspace_supply_demand_task,
+            population_alignment_task,
+        ])
+
     graph.close()
     graph.join()
     LOGGER.info('Finished Urban Nature Access Model')
 
 
+def _admin_level_supply_demand(
+        greenspace_budget_path, population_path, admin_unit_vector_path,
+        target_admin_unit_vector_path):
+    """Calculate average greenspace supply per admin unit.
+
+    Note:
+        The greenspace budget raster and population raster must align
+        perfectly in terms of pixel sizes and raster dimensions, and must have
+        the same projection.
+
+    Args:
+        greenspace_budget_path (string): The path to the greenspace budget
+            raster on disk.
+        population_path (string): The path to the population raster on disk.
+        admin_unit_vector_path (string): The path to the administrative units
+            vector.
+        target_admin_unit_vector_path (string): The path to where the target
+            administrative units vector will be created on disk.
+
+    Returns:
+        ``None``
+    """
+    raster_info = pygeoprocessing.get_raster_info(greenspace_budget_path)
+
+    # Reprojecting the vector here within this function allows us to produce
+    # the vector where it is used.  Otherwise, we'd be reprojecting in a
+    # separate task and then copying it to another filepath.
+    pygeoprocessing.reproject_vector(
+        admin_unit_vector_path, raster_info['projection_wkt'],
+        target_admin_unit_vector_path, driver_name='GPKG')
+
+    target_vector = gdal.OpenEx(target_admin_unit_vector_path, gdal.GA_Update)
+    target_layer = target_vector.GetLayer()
+
+    supply_sum_fieldname = 'average_greenspace_budget'
+    supply_sum_field = ogr.FieldDefn(supply_sum_fieldname, ogr.OFTReal)
+    supply_sum_field.SetWidth(24)
+    supply_sum_field.SetPrecision(11)
+    target_layer.CreateField(supply_sum_field)
+
+    greenspace_stats = pygeoprocessing.zonal_statistics(
+        (greenspace_budget_path, 1), target_admin_unit_vector_path)
+    population_stats = pygeoprocessing.zonal_statistics(
+        (population_path, 1), target_admin_unit_vector_path)
+
+    target_layer.StartTransaction()
+    for feature in target_layer:
+        feature_id = feature.GetFID()
+
+        avg_greenspace_supply_demand = (
+            greenspace_stats[feature_id]['sum'] /
+            population_stats[feature_id]['sum'])
+
+        feature.SetField(supply_sum_fieldname, avg_greenspace_supply_demand)
+        target_layer.SetFeature(feature)
+    target_layer.CommitTransaction()
+
+    target_layer = None
+    target_vector = None
+
+
+def _greenspace_budget_op(greenspace_supply, greenspace_demand):
+    """Calculate the per-capita greenspace budget.
+
+    This is the amount of greenspace that each pixel has above (positive
+    values) or below (negative values) the user-defined ``greenspace_demand``
+    value.
+
+    Args:
+        greenspace_supply (numpy.array): The supply of greenspace available to
+            each person in the population.  This is ``Ai`` in the User's Guide.
+            This matrix must have ``FLOAT32_NODATA`` as its nodata value.
+        greenspace_demand (float): The per-person greenspace requirement, in
+
+    Returns:
+        A ``numpy.array`` of the calculated greenspace budget.
+    """
+    budget = numpy.full(
+        greenspace_supply.shape, FLOAT32_NODATA, dtype=numpy.float32)
+    valid_pixels = ~numpy.isclose(greenspace_supply, FLOAT32_NODATA)
+    budget[valid_pixels] = greenspace_supply[valid_pixels] - greenspace_demand
+    return budget
+
+
+def _greenspace_supply_demand_op(greenspace_budget, population):
+    """Calculate the supply/demand of greenspace per person.
+
+    Args:
+        greenspace_budget (numpy.array): The area of greenspace budgeted to
+            each person, relative to a minimum required per-person area of
+            greenspace.  This matrix must have ``FLOAT32_NODATA`` as its nodata
+            value.  This matrix must be the same size and shape as
+            ``population``.
+        population (numpy.array): Pixel values represent the population count
+            of the pixel.  This matrix must be the same size and shape as
+            ``greenspace_budget``, and must have ``FLOAT32_NODATA`` as its
+            nodata value.
+
+    Returns:
+        A ``numpy.array`` of the area (in square meters) of greenspace supplied
+        to each individual in each pixel.
+    """
+    supply_demand = numpy.full(
+        greenspace_budget.shape, FLOAT32_NODATA, dtype=numpy.float32)
+    valid_pixels = (
+        ~numpy.isclose(greenspace_budget, FLOAT32_NODATA) &
+        ~numpy.isclose(population, FLOAT32_NODATA))
+    supply_demand[valid_pixels] = (
+        greenspace_budget[valid_pixels] * population[valid_pixels])
+    return supply_demand
+
+
 def _greenspace_population_ratio(
         greenspace_area, convolved_population, greenspace_nodata,
         population_nodata):


Test patch:
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index dc57a101e..fdbb0a238 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -10,6 +10,7 @@ import unittest
 
 import numpy
 import pygeoprocessing
+import shapely.geometry
 from osgeo import gdal
 from osgeo import osr
 
@@ -163,6 +164,36 @@ class UNATests(unittest.TestCase):
         self.assertEqual(0.75, kernel_array.max())
         self.assertEqual(0, kernel_array.min())
 
+    def test_greenspace_budgets(self):
+        """UNA: Test the per-capita greenspace budgets functions."""
+        from natcap.invest import urban_nature_access
+
+        nodata = urban_nature_access.FLOAT32_NODATA
+        greenspace_supply = numpy.array([
+            [nodata, 100.5],
+            [75, 100]], dtype=numpy.float32)
+        greenspace_demand = 50
+
+        population = numpy.array([
+            [50, 100],
+            [40.75, nodata]], dtype=numpy.float32)
+
+        greenspace_budget = urban_nature_access._greenspace_budget_op(
+                greenspace_supply, greenspace_demand)
+        expected_greenspace_budget = numpy.array([
+            [nodata, 50.5],
+            [25, 50]], dtype=numpy.float32)
+        numpy.testing.assert_allclose(
+            greenspace_budget, expected_greenspace_budget)
+
+        supply_demand = urban_nature_access._greenspace_supply_demand_op(
+            greenspace_budget, population)
+        expected_supply_demand = numpy.array([
+            [nodata, 100 * 50.5],
+            [25 * 40.75, nodata]], dtype=numpy.float32)
+        numpy.testing.assert_allclose(
+            supply_demand, expected_supply_demand)
+
     def test_model(self):
         """UNA: Run through the model."""
         from natcap.invest import urban_nature_access
@@ -177,6 +208,9 @@ class UNATests(unittest.TestCase):
                 self.workspace_dir, 'lulc_attributes.csv'),
             'decay_function': 'gaussian',
             'search_radius': 100.0,  # meters
+            'greenspace_demand': 100,  # square meters
+            'admin_unit_vector_path': os.path.join(
+                self.workspace_dir, 'admin_units.geojson'),
         }
 
         random.seed(-1)  # for our random number generation
@@ -223,11 +257,21 @@ class UNATests(unittest.TestCase):
                 8,0
                 9,1"""))
 
+        admin_geom = [
+            shapely.geometry.box(
+                *pygeoprocessing.get_raster_info(
+                    args['lulc_raster_path'])['bounding_box'])]
+        pygeoprocessing.shapely_geometry_to_vector(
+            admin_geom, args['admin_unit_vector_path'],
+            population_wkt, 'GeoJSON')
+
         urban_nature_access.execute(args)
 
         # Since we're doing a semi-manual alignment step, assert that the
         # aligned LULC and population rasters have the same pixel sizes, origin
         # and raster dimensions.
+        # TODO: Remove these assertions once we're using align_and_resize and
+        # it works as expected.
         aligned_lulc_raster_info = pygeoprocessing.get_raster_info(
             os.path.join(args['workspace_dir'], 'intermediate',
                          f"aligned_lulc_{args['results_suffix']}.tif"))
@@ -246,3 +290,17 @@ class UNATests(unittest.TestCase):
         numpy.testing.assert_allclose(
             aligned_lulc_raster_info['bounding_box'],
             aligned_population_raster_info['bounding_box'])
+
+        # Check that we're getting the appropriate summary values in the admin
+        # units vector.
+        admin_vector_path = os.path.join(
+            args['workspace_dir'], 'output',
+            f"admin_units_{args['results_suffix']}.gpkg")
+        admin_vector = gdal.OpenEx(admin_vector_path)
+        admin_layer = admin_vector.GetLayer()
+        self.assertEqual(admin_layer.GetFeatureCount(), 1)
+
+        admin_feature = admin_layer.GetFeature(1)
+        numpy.testing.assert_allclose(
+            admin_feature.GetField('average_greenspace_budget'),
+            -17.9078)  # from eyeballing the results; random seed = 1