Merge branch 'main' into feature/plugins

.github/workflows/build-and-test.yml

@@ -313,7 +313,11 @@ jobs:
         uses: ./.github/actions/setup_env
         with:
           python-version: ${{ env.LATEST_SUPPORTED_PYTHON_VERSION }}
-          requirements-files: requirements.txt requirements-dev.txt requirements-docs.txt
+          requirements-files: |
+            requirements.txt
+            requirements-dev.txt
+            requirements-docs.txt
+            constraints_tests.txt
           requirements: ${{ env.CONDA_DEFAULT_DEPENDENCIES }} pandoc
 
       - name: Make install

.readthedocs.yml

@@ -18,5 +18,9 @@ build:
   tools:
     python: "mambaforge-4.10"
   jobs:
+    post_create_environment:
+      - pip install --upgrade-strategy=only-if-needed -r requirements.txt
+      - pip install --upgrade-strategy=only-if-needed -r requirements-dev.txt
+      - pip install --upgrade-strategy=only-if-needed -r requirements-docs.txt
     post_install:
       - make install

.readthedocs_environment.yml

@@ -13,7 +13,3 @@ dependencies:
 - python=3.11
 - gdal>=3.4.2
 - pip
-- pip:
-  - -r requirements.txt
-  - -r requirements-dev.txt
-  - -r requirements-docs.txt

HISTORY.rst

@@ -65,9 +65,17 @@ Unreleased Changes
     * Add support for latest GDAL versions; remove test-specific constraint on
       GDAL versions from invest requirements.
       https://github.com/natcap/invest/issues/916
+    * Updated to Cython 3 (https://github.com/natcap/invest/issues/556)
 * Annual Water Yield
     * Added the results_suffix to a few intermediate files where it was
       missing. https://github.com/natcap/invest/issues/1517
+* Coastal Blue Carbon
+    * Updated model validation to prevent the case where a user provides only
+      one snapshot year and no analysis year
+      (`#1534 <https://github.com/natcap/invest/issues/1534>`_).
+      Also enforces that the analysis year, if provided, is greater than the
+      latest snapshot year. An analysis year equal to the latest snapshot year
+      is no longer allowed.
 * Coastal Vulnerability
     * Fixed a bug in handling ``nan`` as the nodata value of the bathymetry
       raster. ``nan`` pixels will now be propertly ignored before calculating
@@ -87,6 +95,14 @@ Unreleased Changes
       to negative dimensions. https://github.com/natcap/invest/issues/1431
     * Optimized the creation of the summary vector by minimizing the number of
       times the target vector needs to be rasterized.
+* Seasonal Water Yield
+    * Fixed an issue with the precip directory units. Units for these input
+      rasters are now correctly stated as mm/month.
+      https://github.com/natcap/invest/issues/1571
+    * Fixed an issue where the monthly quickflow values were being summed over
+      a block area and not summed pixelwise. This caused the quickflow
+      output ``QF.tif`` to have malformed values.
+      https://github.com/natcap/invest/issues/1541
 * Wind Energy
     * Fixed a bug where some number inputs were not being properly cast to
       ``float`` or ``int`` types. If the inputs happened to be passed as

Makefile

@@ -2,15 +2,15 @@
 DATA_DIR := data
 GIT_SAMPLE_DATA_REPO        := https://bitbucket.org/natcap/invest-sample-data.git
 GIT_SAMPLE_DATA_REPO_PATH   := $(DATA_DIR)/invest-sample-data
-GIT_SAMPLE_DATA_REPO_REV    := 2e7cd618c661ec3f3b2a3bddfd2ce7d4704abc05
+GIT_SAMPLE_DATA_REPO_REV    := ab8c74a62a93fd0019de2bca064abc0a5a07afab
 
 GIT_TEST_DATA_REPO          := https://bitbucket.org/natcap/invest-test-data.git
 GIT_TEST_DATA_REPO_PATH     := $(DATA_DIR)/invest-test-data
 GIT_TEST_DATA_REPO_REV      := 324abde73e1d770ad75921466ecafd1ec6297752
 
 GIT_UG_REPO                 := https://github.com/natcap/invest.users-guide
 GIT_UG_REPO_PATH            := doc/users-guide
-GIT_UG_REPO_REV             := fa6b181d49136089dce56d4ff8f3dcaf12eb4ced
+GIT_UG_REPO_REV             := 0404bc5d4d43085cdc58f50f8fc29944b10cefb1
 
 ENV = "./env"
 ifeq ($(OS),Windows_NT)
@@ -251,6 +251,7 @@ $(INVEST_BINARIES_DIR): | $(DIST_DIR) $(BUILD_DIR)
 	-$(RMDIR) $(INVEST_BINARIES_DIR)
 	$(PYTHON) -m PyInstaller --workpath $(BUILD_DIR)/pyi-build --clean --distpath $(DIST_DIR) exe/invest.spec
 	$(CONDA) list > $(INVEST_BINARIES_DIR)/package_versions.txt
+	$(PYTHON) -m pip list >> $(INVEST_BINARIES_DIR)/package_versions.txt
 	$(INVEST_BINARIES_DIR)/invest list
 
 # Documentation.

pyproject.toml

@@ -45,7 +45,7 @@ invest = "natcap.invest.cli:main"
 # that we can provide a much easier build experience so long as GDAL is
 # available at runtime.
 requires = [
-    'setuptools>=61', 'wheel', 'setuptools_scm>=8.0', 'cython', 'babel',
+    'setuptools>=61', 'wheel', 'setuptools_scm>=8.0', 'cython>=3.0.0', 'babel',
     'oldest-supported-numpy'
 ]
 build-backend = "setuptools.build_meta"

requirements-dev.txt

@@ -10,7 +10,6 @@
 # scripts/convert-requirements-to-conda-yml.py as though it can only be found
 # on pip.
 
-Cython<3.0.0
 virtualenv>=12.0.1
 pytest
 pytest-subtests

setup.py

@@ -2,8 +2,8 @@
 import platform
 import subprocess
 
-import Cython.Build
 import numpy
+from Cython.Build import cythonize
 from setuptools import setup
 from setuptools.command.build_py import build_py as _build_py
 from setuptools.extension import Extension
@@ -46,25 +46,34 @@ def run(self):
 
 setup(
     install_requires=_REQUIREMENTS,
-    ext_modules=[
+    ext_modules=cythonize([
         Extension(
             name=f'natcap.invest.{package}.{module}',
             sources=[f'src/natcap/invest/{package}/{module}.pyx'],
-            include_dirs=[numpy.get_include()],
-            extra_compile_args=compiler_and_linker_args,
+            extra_compile_args=compiler_args + compiler_and_linker_args,
             extra_link_args=compiler_and_linker_args,
-            language='c++'
-        ) for package, module in [
-            ('delineateit', 'delineateit_core'),
-            ('recreation', 'out_of_core_quadtree'),
-            ('scenic_quality', 'viewshed'),
-            ('ndr', 'ndr_core'),
-            ('sdr', 'sdr_core'),
-            ('seasonal_water_yield', 'seasonal_water_yield_core')
+            language='c++',
+            define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
+        ) for package, module, compiler_args in [
+            ('delineateit', 'delineateit_core', []),
+            ('recreation', 'out_of_core_quadtree', []),
+            # clang-14 defaults to -ffp-contract=on, which causes the
+            # arithmetic of A*B+C to be implemented using a contraction, which
+            # causes an unexpected change in the precision in some viewshed
+            # tests on ARM64 (mac M1).  See these issues for more details:
+            #  * https://github.com/llvm/llvm-project/issues/91824
+            #  * https://github.com/natcap/invest/issues/1562
+            #  * https://github.com/natcap/invest/pull/1564/files
+            # Using this flag on gcc and on all versions of clang should work
+            # as expected, with consistent results.
+            ('scenic_quality', 'viewshed', ['-ffp-contract=off']),
+            ('ndr', 'ndr_core', []),
+            ('sdr', 'sdr_core', []),
+            ('seasonal_water_yield', 'seasonal_water_yield_core', [])
         ]
-    ],
+    ], compiler_directives={'language_level': '3'}),
+    include_dirs=[numpy.get_include()],
     cmdclass={
-        'build_ext': Cython.Build.build_ext,
         'build_py': build_py
     }
 )

src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py

@@ -113,8 +113,10 @@
 LOGGER = logging.getLogger(__name__)
 
 INVALID_ANALYSIS_YEAR_MSG = gettext(
-    "Analysis year {analysis_year} must be >= the latest snapshot year "
+    "Analysis year ({analysis_year}) must be greater than the latest snapshot year "
     "({latest_year})")
+MISSING_ANALYSIS_YEAR_MSG = gettext(
+    "Analysis year is required if only one snapshot year is provided.")
 INVALID_TRANSITION_VALUES_MSG = gettext(
     "The transition table expects values of {model_transitions} but found "
     "values of {transition_values}.")
@@ -2188,9 +2190,14 @@ def validate(args, limit_to=None):
             **MODEL_SPEC['args']['landcover_snapshot_csv']
         )['raster_path'].to_dict()
 
+        snapshot_years = set(snapshots.keys())
+        if len(snapshot_years) == 1 and "analysis_year" not in sufficient_keys:
+            validation_warnings.append(
+                (['analysis_year'], MISSING_ANALYSIS_YEAR_MSG))
+
         if ("analysis_year" not in invalid_keys
                 and "analysis_year" in sufficient_keys):
-            if max(set(snapshots.keys())) > int(args['analysis_year']):
+            if max(snapshot_years) >= int(args['analysis_year']):
                 validation_warnings.append((
                     ['analysis_year'],
                     INVALID_ANALYSIS_YEAR_MSG.format(

src/natcap/invest/delineateit/delineateit_core.pyx

@@ -1,4 +1,3 @@
-# cython: language_level=3
 import numpy
 import pygeoprocessing
 cimport numpy
@@ -105,5 +104,3 @@ cpdef cset[cpair[double, double]] calculate_pour_point_array(
     # return set of (x, y) coordinates referenced to the same coordinate system
     # as the original raster
     return pour_points
-
-

src/natcap/invest/ndr/ndr_core.pyx

@@ -1,5 +1,3 @@
-# cython: profile=False
-# cython: language_level=2
 import tempfile
 import logging
 import os
@@ -129,9 +127,9 @@ cdef class _ManagedRaster:
         self.block_xbits = numpy.log2(self.block_xsize)
         self.block_ybits = numpy.log2(self.block_ysize)
         self.block_nx = (
-            self.raster_x_size + (self.block_xsize) - 1) / self.block_xsize
+            self.raster_x_size + (self.block_xsize) - 1) // self.block_xsize
         self.block_ny = (
-            self.raster_y_size + (self.block_ysize) - 1) / self.block_ysize
+            self.raster_y_size + (self.block_ysize) - 1) // self.block_ysize
 
         self.lru_cache = new LRUCache[int, double*](MANAGED_RASTER_N_BLOCKS)
         self.raster_path = <bytes> raster_path
@@ -197,7 +195,7 @@ cdef class _ManagedRaster:
             if dirty_itr != self.dirty_blocks.end():
                 self.dirty_blocks.erase(dirty_itr)
                 block_xi = block_index % self.block_nx
-                block_yi = block_index / self.block_nx
+                block_yi = block_index // self.block_nx
 
                 # we need the offsets to subtract from global indexes for
                 # cached array
@@ -261,7 +259,7 @@ cdef class _ManagedRaster:
 
     cdef void _load_block(self, int block_index) except *:
         cdef int block_xi = block_index % self.block_nx
-        cdef int block_yi = block_index / self.block_nx
+        cdef int block_yi = block_index // self.block_nx
 
         # we need the offsets to subtract from global indexes for cached array
         cdef int xoff = block_xi << self.block_xbits
@@ -322,7 +320,7 @@ cdef class _ManagedRaster:
                     self.dirty_blocks.erase(dirty_itr)
 
                     block_xi = block_index % self.block_nx
-                    block_yi = block_index / self.block_nx
+                    block_yi = block_index // self.block_nx
 
                     xoff = block_xi << self.block_xbits
                     yoff = block_yi << self.block_ybits
@@ -488,7 +486,7 @@ def ndr_eff_calculation(
             # hasn't already been set for processing.
             flat_index = processing_stack.top()
             processing_stack.pop()
-            global_row = flat_index / n_cols
+            global_row = flat_index // n_cols
             global_col = flat_index % n_cols
 
             crit_len = <float>crit_len_raster.get(global_col, global_row)

src/natcap/invest/recreation/out_of_core_quadtree.pyx

@@ -1,5 +1,3 @@
-# cython: profile=True
-# cython: language_level=2
 """A hierarchical spatial index for fast culling of points in 2D space."""
 
 import os
@@ -22,7 +20,7 @@ from osgeo import osr
 cimport numpy
 
 MAX_BYTES_TO_BUFFER = 2**27  # buffer a little over 128 megabytes
-import buffered_numpy_disk_map
+from natcap.invest.recreation import buffered_numpy_disk_map
 _ARRAY_TUPLE_TYPE = (
     buffered_numpy_disk_map.BufferedNumpyDiskMap._ARRAY_TUPLE_TYPE)
 

src/natcap/invest/scenic_quality/viewshed.pyx

@@ -1,5 +1,3 @@
-# coding=UTF-8
-# cython: language_level=2
 """
 Implements the Wang et al (2000) viewshed based on reference planes.
 
@@ -307,9 +305,9 @@ cdef class _ManagedRaster:
         self.block_xbits = numpy.log2(self.block_xsize)
         self.block_ybits = numpy.log2(self.block_ysize)
         self.block_nx = (
-            self.raster_x_size + (self.block_xsize) - 1) / self.block_xsize
+            self.raster_x_size + (self.block_xsize) - 1) // self.block_xsize
         self.block_ny = (
-            self.raster_y_size + (self.block_ysize) - 1) / self.block_ysize
+            self.raster_y_size + (self.block_ysize) - 1) // self.block_ysize
 
         self.lru_cache = new LRUCache[int, double*](MANAGED_RASTER_N_BLOCKS)
         self.raster_path = <bytes> raster_path
@@ -375,7 +373,7 @@ cdef class _ManagedRaster:
             if dirty_itr != self.dirty_blocks.end():
                 self.dirty_blocks.erase(dirty_itr)
                 block_xi = block_index % self.block_nx
-                block_yi = block_index / self.block_nx
+                block_yi = block_index // self.block_nx
 
                 # we need the offsets to subtract from global indexes for
                 # cached array
@@ -439,7 +437,7 @@ cdef class _ManagedRaster:
 
     cdef void _load_block(self, int block_index) except *:
         cdef int block_xi = block_index % self.block_nx
-        cdef int block_yi = block_index / self.block_nx
+        cdef int block_yi = block_index // self.block_nx
 
         # we need the offsets to subtract from global indexes for cached array
         cdef int xoff = block_xi << self.block_xbits
@@ -500,7 +498,7 @@ cdef class _ManagedRaster:
                     self.dirty_blocks.erase(dirty_itr)
 
                     block_xi = block_index % self.block_nx
-                    block_yi = block_index / self.block_nx
+                    block_yi = block_index // self.block_nx
 
                     xoff = block_xi << self.block_xbits
                     yoff = block_yi << self.block_ybits
@@ -815,23 +813,7 @@ def viewshed(dem_raster_path_band,
             if target_distance > max_visible_radius:
                 break
 
-            # This is a weird platform-specific workaround addressing
-            # https://github.com/natcap/invest/issues/1562
-            # On M1 macs, the all-in-one-line addition of _product and r_v
-            # would create small but noticeable numerical error.  Breaking the
-            # calculation onto two lines eliminates the numerical error.  This
-            # behavior is reproducible in C, outside of Cython on an M1 mac.
-            # So, this calculation would introduce error:
-            #   z = (((previous_height-r_v)/slope_distance) * target_distance) + r_v
-            # while the formlation below does not.
-            # For the script used for testing, see
-            #  https://gist.github.com/phargogh/c4264b37e7f0beed31661eacce53d14a
-            #
-            # Some of this may be related to the fact that x86 chips have
-            # extended precision for FPU-based calculations while M1 ARM chips
-            # do not.  Still, that doesn't explain why the error is introduced.
-            _product = (((previous_height-r_v)/slope_distance) * target_distance)
-            z = _product + r_v
+            z = (((previous_height-r_v)/slope_distance) * target_distance) + r_v
 
             # add on refractivity/curvature-of-earth calculations.
             adjustment = 0.0  # increase in required height due to curvature

src/natcap/invest/sdr/sdr_core.pyx

@@ -1,5 +1,3 @@
-# cython: profile=False
-# cython: language_level=3
 import logging
 import os
 
@@ -196,7 +194,7 @@ cdef class _ManagedRaster:
             if dirty_itr != self.dirty_blocks.end():
                 self.dirty_blocks.erase(dirty_itr)
                 block_xi = block_index % self.block_nx
-                block_yi = block_index / self.block_nx
+                block_yi = block_index // self.block_nx
 
                 # we need the offsets to subtract from global indexes for
                 # cached array

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

@@ -84,7 +84,14 @@
             "contents": {
                 # monthly precipitation maps, each file ending in a number 1-12
                 "[MONTH]": {
-                    **spec_utils.PRECIP,
+                    "type": "raster",
+                    "bands": {
+                        1: {
+                            "type": "number",
+                            "units": u.millimeter/u.month,
+                        },
+                    },
+                    "name": gettext("precipitation"),
                     "about": gettext(
                         "Twelve files, one for each month. File names must "
                         "end with the month number (1-12). For example, "
@@ -916,7 +923,7 @@ def execute(args):
             ],
             dependent_task_list=[
                 align_task, flow_dir_task, stream_threshold_task,
-                fill_pit_task, qf_task] + quick_flow_task_list,
+                fill_pit_task] + quick_flow_task_list,
             task_name='calculate local recharge')
 
     # calculate Qb as the sum of local_recharge_avail over the AOI, Eq [9]
@@ -988,7 +995,7 @@ def execute(args):
 
 
 # raster_map equation: sum the monthly qfis
-def qfi_sum_op(*qf_values): return numpy.sum(qf_values)
+def qfi_sum_op(*qf_values): return numpy.sum(qf_values, axis=0)
 
 
 def _calculate_l_avail(l_path, gamma, target_l_avail_path):