Repository: natcap/invest
Pull request #724 (2021-11-15 20:44:55)

Urban Nature Access: boilerplate and population resampling
This PR is the first of several relating to #722 .  This PR specifically starts out by:

* Adding an `ARGS_SPEC` (mostly there, not 100% complete yet)
* Implements and tests a function to resample a population raster without gaining or losing significant numbers of people (this is the only serious computation in the model right now)
* Adds the start to a test for the model.

I haven't yet made any notes in the README and there isn't a UG chapter yet, so both of those checklist items are pending a later PR.

I've already created a new target branch (`feature/urban-nature-access`, based on `release/3.10`) that this PR is set to merge into.  If we end up waiting to merge this into `release/3.11` after the 3.10 release, that'll work too.




Patch:
diff --git a/pyproject.toml b/pyproject.toml
index 2acfe34c7..642226188 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,3 +6,23 @@
 # that we can provide a much easier build experience so long as GDAL is
 # available at runtime.
 requires = ['setuptools', 'wheel', 'setuptools_scm', 'numpy', 'cython']
+
+[tool.isort]
+# isort can be installed from conda-forge
+# flake8-isort can also be used for flake8 linting integration
+#
+# To test isort configuration, run:
+#     isort . --show-config
+# A warning will present if config can't be loaded.
+#
+# Each section of imports will be sorted alphabetically
+force_alphabetical_sort_within_sections = true
+
+# "from osgeo import gdal" and "from osgeo import osr" are on separate lines
+force_single_line = true
+
+# This defines the sections used, and clarifies that we don't want to have
+# lines between the FIRSTPARTY (absolute natcap.invest imports) and the
+# LOCALFOLDER (relative imports, e.g. "from . import validation) section.
+sections = 'FUTURE,STDLIB,THIRDPARTY,FIRSTPARTY,LOCALFOLDER'
+no_lines_before = 'LOCALFOLDER'
diff --git a/src/natcap/invest/__init__.py b/src/natcap/invest/__init__.py
index 9cf31f773..5d01a1c01 100644
--- a/src/natcap/invest/__init__.py
+++ b/src/natcap/invest/__init__.py
@@ -197,6 +197,12 @@ MODEL_METADATA = {
         gui='urban_cooling_model.UrbanCoolingModel',
         userguide='urban_cooling_model.html',
         aliases=('ucm',)),
+    'urban_nature_access': _MODELMETA(
+        model_title='Urban Nature Access',
+        pyname='natcap.invest.urban_nature_access',
+        gui='',  # TODO
+        userguide='',  # TODO,
+        aliases=('una',)),
 }
 
 
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
new file mode 100644
index 000000000..64cc4248e
--- /dev/null
+++ b/src/natcap/invest/urban_nature_access.py
@@ -0,0 +1,321 @@
+import logging
+import os
+import shutil
+import tempfile
+
+import numpy
+import pygeoprocessing
+import taskgraph
+from osgeo import gdal
+from osgeo import osr
+
+from natcap.invest import MODEL_METADATA
+from . import spec_utils
+from . import utils
+from . import validation
+from .spec_utils import u
+
+LOGGER = logging.getLogger(__name__)
+UINT32_NODATA = int(numpy.iinfo(numpy.uint32).max)
+FLOAT32_NODATA = float(numpy.finfo(numpy.float32).max)
+ARGS_SPEC = {
+    'model_name': MODEL_METADATA['urban_nature_access'].model_title,
+    'pyname': MODEL_METADATA['urban_nature_access'].pyname,
+    'userguide_html': MODEL_METADATA['urban_nature_access'].userguide,
+    'args_with_spatial_overlap': {
+        'spatial_keys': [
+            'lulc_raster_path', 'population_raster_path',
+            'admin_unit_vector_path'],
+        'different_projections_ok': True,
+    },
+    'args': {
+        'workspace_dir': spec_utils.WORKSPACE,
+        'results_suffix': spec_utils.SUFFIX,
+        'n_workers': spec_utils.N_WORKERS,
+        'lulc_raster_path': {
+            **spec_utils.LULC,
+            'projected': True,
+            'projection_units': u.meter,
+            'about': (
+                "A map of LULC codes. "
+                "All values in this raster must have corresponding entries "
+                "in the LULC attribute table."),
+        },
+        'lulc_attribute_table': {
+            'name': 'LULC attribute table',
+            'type': 'csv',
+            'columns': {
+                'lucode': {
+                    'type': 'integer',
+                    'about': (
+                        "LULC code.  Every value in the LULC map must have a "
+                        "corresponding entry in this column."),
+                },
+                'greenspace': {
+                    'type': 'number',
+                    'units': u.none,
+                    'about': (
+                        "1 if this landcover code represents greenspace, 0 "
+                        "if not."
+                    ),
+                }
+            },
+            'about': (
+                "A table identifying which LULC codes represent greenspace."
+            ),
+        },
+        'population_raster_path': {
+            'type': 'raster',
+            'name': 'population raster',
+            'bands': {
+                1: {'type': 'number', 'units': u.none}
+            },
+            'projected': True,
+            'projection_units': u.meter,
+            'about': (
+                "A raster representing the number of people who live in each "
+                "pixel."
+            ),
+        },
+        'admin_unit_vector_path': {
+            'type': 'vector',
+            'name': 'administrative boundaries',
+            'geometries': spec_utils.POLYGONS,
+            'fields': {},  # TODO, complete required fields (if any)
+            'about': "",  # TODO, will know more about this when I implement.
+        },
+        'greenspace_demand': {
+            'type': 'number',
+            '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.
+        },
+        'search_radius': {
+            'type': 'number',
+            'name': 'search radius',
+            'units': u.m,
+            'expression': "value > 0",
+            'about': "",  # TODO, will know more about this when I implement.
+        }
+    }
+}
+
+
+_OUTPUT_BASE_FILES = {}
+_INTERMEDIATE_BASE_FILES = {
+    'resampled_population': 'population_resampled.tif',
+}
+
+
+def execute(args):
+    """Urban Nature Access.
+
+    Args:
+        args['workspace_dir'] (string): (required) Output directory for
+            intermediate, temporary and final files.
+        args['results_suffix'] (string): (optional) String to append to any
+            output file.
+        args['n_workers'] (int): (optional) The number of worker processes to
+            use for executing the tasks of this model.  If omitted, computation
+            will take place in the current process.
+        args['lulc_raster_path'] (string): (required) A string path to a
+            GDAL-compatible land-use/land-cover raster containing integer
+            landcover codes.  Must be linearly projected in meters.
+        args['lulc_attribute_table'] (string): (required) A string path to a
+            CSV with the following columns:
+
+            * ``lucode``: the integer landcover code represented.
+            * ``greenspace``: ``0`` or ``1`` indicating whether this landcover
+              code is (``1``) or is not (``0``) a greenspace pixel.
+
+        args['population_raster_path'] (string): (required) A string path to a
+            GDAL-compatible raster where pixels represent the population of
+            that pixel.  Must be linearly projected in meters.
+        args['admin_unit_vector_path'] (string): (required) A string path to a
+            GDAL-compatible vector containing polygon administrative
+            boundaries.
+        args['greenspace_demand'] (number): (required) A positive, nonzero
+            number indicating the required greenspace, in m² per capita.
+        args['search_radius'] (number): (required) A positive, nonzero number
+            indicating the maximum distance that people travel for recreation.
+
+    Returns:
+        ``None``
+    """
+    LOGGER.info('Starting Urban Nature Access Model')
+
+    output_dir = os.path.join(args['workspace_dir'], 'output')
+    intermediate_dir = os.path.join(args['workspace_dir'], 'intermediate')
+    utils.make_directories([output_dir, intermediate_dir])
+
+    file_suffix = utils.make_suffix_string(args, 'results_suffix')
+    file_registry = utils.build_file_registry(
+        [(_OUTPUT_BASE_FILES, output_dir),
+         (_INTERMEDIATE_BASE_FILES, intermediate_dir)],
+        file_suffix)
+
+    work_token_dir = os.path.join(intermediate_dir, '_taskgraph_working_dir')
+    try:
+        n_workers = int(args['n_workers'])
+    except (KeyError, ValueError, TypeError):
+        # KeyError when n_workers is not present in args
+        # ValueError when n_workers is an empty string.
+        # TypeError when n_workers is None.
+        n_workers = -1  # Synchronous execution
+    graph = taskgraph.TaskGraph(work_token_dir, n_workers)
+
+    # Align the population raster to the LULC.
+    lulc_raster_info = pygeoprocessing.get_raster_info(
+        args['lulc_raster_path'])
+
+    population_alignment_task = graph.add_task(
+        _resample_population_raster,
+        kwargs={
+            'source_population_raster_path': args['population_raster_path'],
+            'target_population_raster_path': file_registry[
+                'resampled_population'],
+            'lulc_pixel_size': lulc_raster_info['pixel_size'],
+            'lulc_bb': lulc_raster_info['bounding_box'],
+            'lulc_projection_wkt': lulc_raster_info['projection_wkt'],
+            'working_dir': intermediate_dir,
+        },
+        target_path_list=[file_registry['resampled_population']],
+        task_name='Resample population to LULC resolution')
+
+    graph.close()
+    graph.join()
+    LOGGER.info('Finished Urban Nature Access Model')
+
+
+def _resample_population_raster(
+        source_population_raster_path, target_population_raster_path,
+        lulc_pixel_size, lulc_bb, lulc_projection_wkt, working_dir):
+    """Resample a population raster without losing or gaining people.
+
+    Population rasters are an interesting special case where the data are
+    neither continuous nor categorical, and the total population count
+    typically matters.  Common resampling methods for continuous
+    (interpolation) and categorical (nearest-neighbor) datasets leave room for
+    the total population of a resampled raster to significantly change.  This
+    function resamples a population raster with the following steps:
+
+        1. Convert a population count raster to population density per pixel
+        2. Warp the population density raster to the target spatial reference
+           and pixel size using bilinear interpolation.
+        3. Convert the warped density raster back to population counts.
+
+    Args:
+        source_population_raster_path (string): The source population raster.
+            Pixel values represent the number of people occupying the pixel.
+            Must be linearly projected in meters.
+        target_population_raster_path (string): The path to where the target,
+            warped population raster will live on disk.
+        lulc_pixel_size (tuple): A tuple of the pixel size for the target
+            raster.  Passed directly to ``pygeoprocessing.warp_raster``.
+        lulc_bb (tuple): A tuple of the bounding box for the target raster.
+            Passed directly to ``pygeoprocessing.warp_raster``.
+        lulc_projection_wkt (string): The Well-Known Text of the target
+            spatial reference fro the target raster.  Passed directly to
+            ``pygeoprocessing.warp_raster``.  Assumed to be a linear projection
+            in meters.
+        working_dir (string): The path to a directory on disk.  A new directory
+            is created within this directory for the storage of temporary files
+            and then deleted upon successful completion of the function.
+
+    Returns:
+        ``None``
+    """
+    if not os.path.isdir(working_dir):
+        os.makedirs(working_dir)
+    tmp_working_dir = tempfile.mkdtemp(dir=working_dir)
+    population_raster_info = pygeoprocessing.get_raster_info(
+        source_population_raster_path)
+    pixel_area = numpy.multiply(*population_raster_info['pixel_size'])
+    population_nodata = population_raster_info['nodata'][0]
+
+    population_srs = osr.SpatialReference()
+    population_srs.ImportFromWkt(population_raster_info['projection_wkt'])
+
+    # Convert population pixel area to square km
+    population_pixel_area = (
+        pixel_area * population_srs.GetLinearUnits()) / 1e6
+
+    def _convert_population_to_density(population):
+        """Convert population counts to population per square km.
+
+        Args:
+            population (numpy.array): A numpy array where pixel values
+                represent the number of people who reside in a pixel.
+
+        Returns:
+            """
+        out_array = numpy.full(
+            population.shape, FLOAT32_NODATA, dtype=numpy.float32)
+
+        valid_mask = slice(None)
+        if population_nodata is not None:
+            valid_mask = ~numpy.isclose(population, population_nodata)
+
+        out_array[valid_mask] = population[valid_mask] / population_pixel_area
+        return out_array
+
+    # Step 1: convert the population raster to population density per sq. km
+    density_raster_path = os.path.join(tmp_working_dir, 'pop_density.tif')
+    pygeoprocessing.raster_calculator(
+        [(source_population_raster_path, 1)],
+        _convert_population_to_density,
+        density_raster_path, gdal.GDT_Float32, FLOAT32_NODATA)
+
+    # Step 2: align to the LULC
+    warped_density_path = os.path.join(tmp_working_dir, 'warped_density.tif')
+    pygeoprocessing.warp_raster(
+        density_raster_path,
+        target_pixel_size=lulc_pixel_size,
+        target_raster_path=warped_density_path,
+        resample_method='bilinear',
+        target_bb=lulc_bb,
+        target_projection_wkt=lulc_projection_wkt)
+
+    # Step 3: convert the warped population raster back from density to the
+    # population per pixel
+    target_srs = osr.SpatialReference()
+    target_srs.ImportFromWkt(lulc_projection_wkt)
+    # Calculate target pixel area in km to match above
+    target_pixel_area = (
+        numpy.multiply(*lulc_pixel_size) * target_srs.GetLinearUnits()) / 1e6
+
+    def _convert_density_to_population(density):
+        """Convert a population density raster back to population counts.
+
+        Args:
+            density (numpy.array): An array of the population density per
+                square kilometer.
+
+        Returns:
+            A ``numpy.array`` of the population counts given the target pixel
+            size of the output raster."""
+        # We're using a float32 array here because doing these unit
+        # conversions is likely to end up with partial people spread out
+        # between multiple pixels.  So it's preserving an unrealistic degree of
+        # precision, but that's probably OK because pixels are imprecise
+        # measures anyways.
+        out_array = numpy.full(
+            density.shape, FLOAT32_NODATA, dtype=numpy.float32)
+
+        # We already know that the nodata value is FLOAT32_NODATA
+        valid_mask = ~numpy.isclose(density, FLOAT32_NODATA)
+        out_array[valid_mask] = density[valid_mask] * target_pixel_area
+        return out_array
+
+    pygeoprocessing.raster_calculator(
+        [(warped_density_path, 1)],
+        _convert_density_to_population,
+        target_population_raster_path, gdal.GDT_Float32, FLOAT32_NODATA)
+
+    shutil.rmtree(tmp_working_dir, ignore_errors=True)
+
+
+def validate(args, limit_to=None):
+    return validation.validate(
+        args, ARGS_SPEC['args'], ARGS_SPEC['args_with_spatial_overlap'])


Test patch:
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
new file mode 100644
index 000000000..99dede28c
--- /dev/null
+++ b/tests/test_urban_nature_access.py
@@ -0,0 +1,134 @@
+# coding=UTF-8
+"""Tests for the Urban Nature Access Model."""
+import os
+import random
+import shutil
+import tempfile
+import unittest
+
+import numpy
+import pygeoprocessing
+from osgeo import gdal
+from osgeo import osr
+
+_DEFAULT_ORIGIN = (444720, 3751320)
+_DEFAULT_PIXEL_SIZE = (30, -30)
+_DEFAULT_EPSG = 3116
+
+
+class UNATests(unittest.TestCase):
+    """Tests for the Urban Nature Access Model."""
+
+    def setUp(self):
+        """Override setUp function to create temp workspace directory."""
+        # this lets us delete the workspace after its done no matter the
+        # the test result
+        self.workspace_dir = tempfile.mkdtemp(suffix='\U0001f60e')  # smiley
+
+    def tearDown(self):
+        """Override tearDown function to remove temporary directory."""
+        shutil.rmtree(self.workspace_dir)
+
+    def test_resample_population_raster(self):
+        """UNA: Test population raster resampling."""
+        from natcap.invest import urban_nature_access
+
+        random.seed(-1)  # for our random number generation
+
+        source_population_raster_path = os.path.join(
+            self.workspace_dir, 'population.tif')
+        population_pixel_size = (90, -90)
+        population_array_shape = (10, 10)
+
+        array_of_100s = numpy.full(
+            population_array_shape, 100, dtype=numpy.uint32)
+        array_of_random_ints = numpy.array(
+            random.choices(range(0, 100), k=100),
+            dtype=numpy.uint32).reshape(population_array_shape)
+
+        for population_array in (
+                array_of_100s, array_of_random_ints):
+            population_srs = osr.SpatialReference()
+            population_srs.ImportFromEPSG(_DEFAULT_EPSG)
+            population_wkt = population_srs.ExportToWkt()
+            pygeoprocessing.numpy_array_to_raster(
+                base_array=population_array,
+                target_nodata=-1,
+                pixel_size=population_pixel_size,
+                origin=_DEFAULT_ORIGIN,
+                projection_wkt=population_wkt,
+                target_path=source_population_raster_path)
+
+            for target_pixel_size in (
+                    (30, -30),  # 1/3 the pixel size
+                    (4, -4),  # way smaller
+                    (100, -100)):  # bigger
+                target_population_raster_path = os.path.join(
+                    self.workspace_dir, 'resampled_population.tif')
+                urban_nature_access._resample_population_raster(
+                    source_population_raster_path,
+                    target_population_raster_path,
+                    lulc_pixel_size=target_pixel_size,
+                    lulc_bb=pygeoprocessing.get_raster_info(
+                        source_population_raster_path)['bounding_box'],
+                    lulc_projection_wkt=population_wkt,
+                    working_dir=os.path.join(self.workspace_dir, 'working'))
+
+                resampled_population_array = (
+                    pygeoprocessing.raster_to_numpy_array(
+                        target_population_raster_path))
+
+                # There should be no significant loss or gain of population due
+                # to warping, but the fact that this is aggregating across the
+                # whole raster (lots of pixels) means we need to lower the
+                # relative tolerance.
+                numpy.testing.assert_allclose(
+                    population_array.sum(), resampled_population_array.sum(),
+                    rtol=1e-3)
+
+    def test_model(self):
+        """UNA: Run through the model."""
+        from natcap.invest import urban_nature_access
+
+        args = {
+            'workspace_dir': os.path.join(self.workspace_dir, 'workspace'),
+            'results_suffix': 'suffix',
+            'population_raster_path': os.path.join(
+                self.workspace_dir, 'population.tif'),
+            'lulc_raster_path': os.path.join(self.workspace_dir, 'lulc.tif'),
+        }
+
+        random.seed(-1)  # for our random number generation
+        population_pixel_size = (90, -90)
+        population_array_shape = (10, 10)
+        population_array = numpy.array(
+            random.choices(range(0, 100), k=100),
+            dtype=numpy.int32).reshape(population_array_shape)
+        population_srs = osr.SpatialReference()
+        population_srs.ImportFromEPSG(_DEFAULT_EPSG)
+        population_wkt = population_srs.ExportToWkt()
+        pygeoprocessing.numpy_array_to_raster(
+            base_array=population_array,
+            target_nodata=-1,
+            pixel_size=population_pixel_size,
+            origin=_DEFAULT_ORIGIN,
+            projection_wkt=population_wkt,
+            target_path=args['population_raster_path'])
+
+        lulc_pixel_size = _DEFAULT_PIXEL_SIZE
+        lulc_array_shape = (30, 30)
+        lulc_array = numpy.array(
+            random.choices(range(0, 10), k=900),
+            dtype=numpy.int32).reshape(lulc_array_shape)
+        pygeoprocessing.numpy_array_to_raster(
+            base_array=lulc_array,
+            target_nodata=-1,
+            pixel_size=lulc_pixel_size,
+            origin=_DEFAULT_ORIGIN,
+            projection_wkt=population_wkt,
+            target_path=args['lulc_raster_path'])
+
+        urban_nature_access.execute(args)
+
+        # TODO: Assertions will be added here later when I have something to
+        # test.