Repository: natcap/invest (python)

Issue:
Urban Nature Access: support a proportion of urban nature instead of binary 1/0
@lmandle suggested that we modify the Urban Nature Access to support a proportion of urban nature in a pixel rather than a strict binary 1/0 classification.

It would be helpful to spec this out a bit more and be sure that we update all the needed places:
- [ ] model source
- [ ] MODEL_SPEC
- [ ] sample data
- [ ] user's guide

Pull request:
Redefine Urban Nature pixels to have a proportion of Urban Nature
This PR implements urban nature as a proportion of the pixel area rather than assuming that the pixel is entirely urban nature.  UG updates are in progress and will be complete before the release.  No updates to sample data are necessary because we'll be using peer-reviewed values in a forthcoming PR (see https://github.com/natcap/invest/issues/1354).


Fixes #1180

## 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 e6df8937d..b48c6f833 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -100,11 +100,7 @@ Unreleased Changes
       consequence criteria were skipped for a single habitat. The model now
       correctly handles this case. https://github.com/natcap/invest/issues/1250
     * Tables in the .xls format are no longer supported. This format was
-      deprecated by ``pandas``.
-      (`#1271 <https://github.com/natcap/invest/issues/1271>`_)
-    * Fixed a bug where vector inputs could be rasterized onto a grid that is
-      not exactly aligned with other raster inputs.
-      (`#1312 <https://github.com/natcap/invest/issues/1312>`_)
+      deprecated by ``pandas``. (`#1271 <https://github.com/natcap/invest/issues/1271>`_)
 * NDR
     * Fixed a bug where results were calculated incorrectly if the runoff proxy
       raster (or the DEM or LULC) had no nodata value
@@ -186,6 +182,15 @@ Unreleased Changes
         * When defining search radii for population groups, one new output
           raster is created for each population group.  These files are named
           ``accessible_urban_nature_to_[POP_GROUP].tif``.
+
+    * Urban nature classes can now be defined to occupy a proportion of a
+      pixel, such as a park that is semi-developed. This proportion is
+      provided through user input as a proportion (0-1) in the
+      ``urban_nature`` column of the LULC Attribute Table.  A value of ``0``
+      indicates that there is no urban nature in this class, ``0.333``
+      indicates that a third of the area of this LULC class is urban nature,
+      and ``1`` would indicate that the entire LULC class's area is urban
+      nature.  https://github.com/natcap/invest/issues/1180
 * Visitation: Recreation and Tourism
     * Fixed a bug where overlapping predictor polygons would be double-counted
       in ``polygon_area_coverage`` and ``polygon_percent_coverage`` calculations.
diff --git a/src/natcap/invest/hra.py b/src/natcap/invest/hra.py
index 481c04f65..10b54b6eb 100644
--- a/src/natcap/invest/hra.py
+++ b/src/natcap/invest/hra.py
@@ -1585,7 +1585,7 @@ def _align(raster_path_map, vector_path_map, target_pixel_size,
             layer = None
             vector = None
 
-            pygeoprocessing.create_raster_from_bounding_box(
+            _create_raster_from_bounding_box(
                 target_raster_path=target_raster_path,
                 target_bounding_box=target_bounding_box,
                 target_pixel_size=target_pixel_size,
@@ -1600,6 +1600,74 @@ def _align(raster_path_map, vector_path_map, target_pixel_size,
                 burn_values=burn_values, option_list=rasterize_option_list)
 
 
+def _create_raster_from_bounding_box(
+        target_raster_path, target_bounding_box, target_pixel_size,
+        target_pixel_type, target_srs_wkt, target_nodata=None,
+        fill_value=None):
+    """Create a raster from a given bounding box.
+
+    Args:
+        target_raster_path (string): The path to where the new raster should be
+            created on disk.
+        target_bounding_box (tuple): a 4-element iterable of (minx, miny,
+            maxx, maxy) in projected units matching the SRS of
+            ``target_srs_wkt``.
+        target_pixel_size (tuple): A 2-element tuple of the (x, y) pixel size
+            of the target raster.  Elements are in units of the target SRS.
+        target_pixel_type (int): The GDAL GDT_* type of the target raster.
+        target_srs_wkt (string): The SRS of the target raster, in Well-Known
+            Text format.
+        target_nodata (float): If provided, the nodata value of the target
+            raster.
+        fill_value=None (number): If provided, the value that the target raster
+            should be filled with.
+
+    Returns:
+        ``None``
+    """
+    bbox_minx, bbox_miny, bbox_maxx, bbox_maxy = target_bounding_box
+
+    driver = gdal.GetDriverByName('GTiff')
+    n_bands = 1
+    n_cols = int(numpy.ceil(
+        abs((bbox_maxx - bbox_minx) / target_pixel_size[0])))
+    n_rows = int(numpy.ceil(
+        abs((bbox_maxy - bbox_miny) / target_pixel_size[1])))
+
+    raster = driver.Create(
+        target_raster_path, n_cols, n_rows, n_bands, target_pixel_type,
+        options=['TILED=YES', 'BIGTIFF=YES', 'COMPRESS=DEFLATE',
+                 'BLOCKXSIZE=256', 'BLOCKYSIZE=256'])
+    raster.SetProjection(target_srs_wkt)
+
+    # Set the transform based on the upper left corner and given pixel
+    # dimensions.  Bounding box is in format [minx, miny, maxx, maxy]
+    if target_pixel_size[0] < 0:
+        x_source = bbox_maxx
+    else:
+        x_source = bbox_minx
+    if target_pixel_size[1] < 0:
+        y_source = bbox_maxy
+    else:
+        y_source = bbox_miny
+    raster_transform = [
+        x_source, target_pixel_size[0], 0.0,
+        y_source, 0.0, target_pixel_size[1]]
+    raster.SetGeoTransform(raster_transform)
+
+    # Fill the band if requested.
+    band = raster.GetRasterBand(1)
+    if fill_value is not None:
+        band.Fill(fill_value)
+
+    # Set the nodata value.
+    if target_nodata is not None:
+        band.SetNoDataValue(float(target_nodata))
+
+    band = None
+    raster = None
+
+
 def _simplify(source_vector_path, tolerance, target_vector_path,
               preserve_columns=None):
     """Simplify a geometry to a given tolerance.
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index 4da0cd6a8..bff4ddafc 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -1,5 +1,4 @@
 import collections
-import functools
 import logging
 import math
 import os
@@ -80,11 +79,13 @@ MODEL_SPEC = {
             'columns': {
                 'lucode': spec_utils.LULC_TABLE_COLUMN,
                 'urban_nature': {
-                    'type': 'number',
-                    'units': u.none,
+                    'type': 'ratio',
                     'about': (
-                        "Binary code indicating whether the LULC type is "
-                        "(1) or is not (0) an urban nature type."
+                        "The proportion (0-1) indicating how much of the land "
+                        "in this LULC type is urban nature. "
+                        "0 indicates no area this LULC type is urban nature, "
+                        "1 indicates that this LULC type is entirely urban "
+                        "nature."
                     ),
                 },
                 'search_radius_m': {
@@ -648,9 +649,10 @@ def execute(args):
             CSV with the following columns:
 
             * ``lucode``: (required) the integer landcover code represented.
-            * ``urban_nature``: (required) ``0`` or ``1`` indicating whether
-              this landcover code is (``1``) or is not (``0``) an urban nature
-              pixel.
+            * ``urban_nature``: (required) a proportion (0-1) representing
+              how much of this landcover type is urban nature.  ``0``
+              indicates none of this type's area is urban nature, ``1``
+              indicates all of this type's area is urban nature.
             * ``search_radius_m``: (conditionally required) the search radius
               for this urban nature LULC class in meters. Required for all
               urban nature LULC codes if ``args['search_radius_mode'] ==
@@ -939,7 +941,7 @@ def execute(args):
     if args['search_radius_mode'] == RADIUS_OPT_UNIFORM:
         search_radii = set([float(args['search_radius'])])
     elif args['search_radius_mode'] == RADIUS_OPT_URBAN_NATURE:
-        urban_nature_attrs = attr_table[attr_table['urban_nature'] == 1]
+        urban_nature_attrs = attr_table[attr_table['urban_nature'] > 0]
         try:
             search_radii = set(urban_nature_attrs['search_radius_m'].unique())
         except KeyError as missing_key:
@@ -1802,13 +1804,16 @@ def _reclassify_urban_nature_area(
     """Reclassify LULC pixels into the urban nature area they represent.
 
     After execution, urban nature pixels will have values representing the
-    pixel's area, while pixels that are not urban nature will have a pixel
-    value of 0.  Nodata values will propagate to the output raster.
+    pixel's area of urban nature (pixel area * proportion of urban nature),
+    while pixels that are not urban nature will have a pixel value of 0.
+    Nodata values will propagate to the output raster.
 
     Args:
         lulc_raster_path (string): The path to a land-use/land-cover raster.
         lulc_attribute_table (string): The path to a CSV table representing
             LULC attributes.  Must have "lucode" and "urban_nature" columns.
+            The "urban_nature" column represents a proportion 0-1 of how much
+            of the pixel's area represents urban nature.
         target_raster_path (string): Where the reclassified urban nature raster
             should be written.
         only_these_urban_nature_codes=None (iterable or None): If ``None``, all
@@ -1830,13 +1835,15 @@ def _reclassify_urban_nature_area(
         valid_urban_nature_codes = set(only_these_urban_nature_codes)
     else:
         valid_urban_nature_codes = set(
-            lulc_attribute_df[lulc_attribute_df['urban_nature'] == 1].index)
+            lulc_attribute_df[lulc_attribute_df['urban_nature'] > 0].index)
 
     urban_nature_area_map = {}
-    for lucode in lulc_attribute_df.index:
+    for row in lulc_attribute_df[['urban_nature']].itertuples():
+        lucode = row.Index
+        urban_nature_proportion = row.urban_nature
         urban_nature_area = 0
         if lucode in valid_urban_nature_codes:
-            urban_nature_area = squared_pixel_area
+            urban_nature_area = squared_pixel_area * urban_nature_proportion
         urban_nature_area_map[lucode] = urban_nature_area
 
     lulc_raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)


Test patch:
diff --git a/tests/test_hra.py b/tests/test_hra.py
index 0999810ea..43a8c4e4a 100644
--- a/tests/test_hra.py
+++ b/tests/test_hra.py
@@ -514,6 +514,40 @@ class HRAUnitTests(unittest.TestCase):
             (source_array != nodata).astype(numpy.uint8)
         )
 
+    def test_create_raster_from_bounding_box(self):
+        """HRA: test creation of a raster from a bbox."""
+        from natcap.invest import hra
+
+        # [minx, miny, maxx, maxy]
+        bounding_box = [
+            ORIGIN[0],
+            ORIGIN[1] - 100,  # force rounding up of pixel dimensions
+            ORIGIN[0] + 90,  # no rounding up needed
+            ORIGIN[1],
+        ]
+        pixel_size = (30, -30)
+        target_raster_path = os.path.join(self.workspace_dir, 'raster.tif')
+        hra._create_raster_from_bounding_box(
+            target_raster_path, bounding_box, pixel_size, gdal.GDT_Byte,
+            SRS_WKT, target_nodata=2, fill_value=2)
+
+        try:
+            raster = gdal.OpenEx(target_raster_path)
+            band = raster.GetRasterBand(1)
+            self.assertEqual(
+                raster.GetGeoTransform(),
+                (ORIGIN[0], pixel_size[0], 0.0, ORIGIN[1], 0.0, pixel_size[1])
+            )
+            self.assertEqual(raster.RasterXSize, 3)
+            self.assertEqual(raster.RasterYSize, 4)
+            self.assertEqual(band.GetNoDataValue(), 2)
+            numpy.testing.assert_array_equal(
+                band.ReadAsArray(),
+                numpy.full((4, 3), 2, dtype=numpy.uint8))
+        finally:
+            band = None
+            raster = None
+
     def test_align(self):
         """HRA: test alignment function."""
         from natcap.invest import hra
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index 49cd88d47..eac33a3ea 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -85,7 +85,8 @@ def _build_model_args(workspace):
             6,0,100
             7,1,100
             8,0,100
-            9,1,100"""))
+            9,1,100
+            """))
 
     admin_geom = [
         shapely.geometry.box(
@@ -956,9 +957,76 @@ class UNATests(unittest.TestCase):
         numpy.testing.assert_allclose(
             numpy.sum(weighted_sum_array[~nodata_pixels]), 1122.5)
 
+    def test_urban_nature_proportion(self):
+        """UNA: Run the model with urban nature proportion."""
+        from natcap.invest import urban_nature_access
+
+        args = _build_model_args(self.workspace_dir)
+        args['search_radius_mode'] = urban_nature_access.RADIUS_OPT_UNIFORM
+        args['search_radius'] = 1000
+        with open(args['lulc_attribute_table'], 'a') as attr_table:
+            attr_table.write("10,0.5,100\n")
+
+        # make sure our inputs validate
+        validation_results = urban_nature_access.validate(args)
+        self.assertEqual(validation_results, [])
+
+        urban_nature_access.execute(args)
+
+    def test_reclassify_urban_nature(self):
+        """UNA: Test for urban nature area reclassification."""
+        from natcap.invest import urban_nature_access
+        args = _build_model_args(self.workspace_dir)
+
+        # Rewrite the lulc attribute table to use proportions of urban nature.
+        with open(args['lulc_attribute_table'], 'w') as attr_table:
+            attr_table.write(textwrap.dedent(
+                """\
+                lucode,urban_nature,search_radius_m
+                0,0,100
+                1,0.1,100
+                2,0,100
+                3,0.3,100
+                4,0,100
+                5,0.5,100
+                6,0,100
+                7,0.7,100
+                8,0,100
+                9,0.9,100
+                """))
+
+        urban_nature_area_path = os.path.join(
+            self.workspace_dir, 'urban_nature_area.tif')
+
+        for limit_to_lucodes in (None, set([1, 3])):
+            urban_nature_access._reclassify_urban_nature_area(
+                args['lulc_raster_path'], args['lulc_attribute_table'],
+                urban_nature_area_path,
+                only_these_urban_nature_codes=limit_to_lucodes)
+
+            # The source lulc is randomized, so need to programmatically build
+            # up the expected array.
+            source_lulc_array = pygeoprocessing.raster_to_numpy_array(
+                args['lulc_raster_path'])
+            pixel_area = abs(_DEFAULT_PIXEL_SIZE[0] * _DEFAULT_PIXEL_SIZE[1])
+            expected_array = numpy.zeros(source_lulc_array.shape,
+                                         dtype=numpy.float32)
+            for i in range(1, 10, 2):
+                if limit_to_lucodes is not None:
+                    if i not in limit_to_lucodes:
+                        continue
+                factor = float(f"0.{i}")
+                expected_array[source_lulc_array == i] = factor * pixel_area
+
+            reclassified_array = pygeoprocessing.raster_to_numpy_array(
+                urban_nature_area_path)
+            numpy.testing.assert_array_almost_equal(
+                reclassified_array, expected_array)
+
     def test_validate(self):
         """UNA: Basic test for validation."""
         from natcap.invest import urban_nature_access
         args = _build_model_args(self.workspace_dir)
-        args['search_radius_mode'] = urban_nature_access.RADIUS_OPT_URBAN_NATURE
+        args['search_radius_mode'] = (
+            urban_nature_access.RADIUS_OPT_URBAN_NATURE)
         self.assertEqual(urban_nature_access.validate(args), [])