Skip to content

Commit

Permalink
Merge branch 'bugfix/1509-una-validation-missing-uniform-search-radiu…
Browse files Browse the repository at this point in the history 
…s' of github.com:phargogh/invest into bugfix/1509-una-validation-missing-uniform-search-radius
  • Loading branch information
@phargogh
phargogh committed Feb 14, 2024
2 parents 21ae820 + 85d8028 commit b5a98a4
Show file tree
Hide file tree
Showing 5 changed files with 135 additions and 86 deletions.
16 changes: 16 additions & 0 deletions HISTORY.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,22 @@ Unreleased Changes
* Fixed an issue where validation was failing to catch missing values in
the uniform search radius args key when using uniform search radii.
https://github.com/natcap/invest/issues/1509
* Fixed an issue where the output administrative units vector's
``Pund_adm`` and ``Povr_adm`` fields representing undersupplied and
oversupplied populations, respectively, had values of 0 when running the
model with search radii defined per population group. The output
administrative units vector now has the correct values for these fields,
consistent with the user's guide chapter.
https://github.com/natcap/invest/issues/1512
* Fixed an issue where certain nodata values were not being handled
correctly, leading to pixel values of +/- infinity in the urban nature
balance output raster. https://github.com/natcap/invest/issues/1519
* SDR
* Fixed an issue encountered in the sediment deposition function where
rasters with more than 2^32 pixels would raise a cryptic error relating
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.

3.14.1 (2023-12-18)
-------------------
Expand Down
35 changes: 23 additions & 12 deletions src/natcap/invest/sdr/sdr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@

from .. import gettext
from .. import spec_utils
from .. import urban_nature_access
from .. import utils
from .. import validation
from ..model_metadata import MODEL_METADATA
Expand Down Expand Up @@ -609,7 +610,9 @@ def execute(args):
target_pixel_size = (min_pixel_size, -min_pixel_size)

target_sr_wkt = dem_raster_info['projection_wkt']
vector_mask_options = {'mask_vector_path': args['watersheds_path']}
vector_mask_options = {
'mask_vector_path': args['watersheds_path'],
}
align_task = task_graph.add_task(
func=pygeoprocessing.align_and_resize_raster_stack,
args=(
Expand Down Expand Up @@ -1491,18 +1494,26 @@ def _generate_report(
target_layer = target_vector.GetLayer()
target_layer.SyncToDisk()

# It's worth it to check if the geometries don't significantly overlap.
# On large rasters, this can save a TON of time rasterizing even a
# relatively simple vector.
geometries_might_overlap = urban_nature_access._geometries_overlap(
watershed_results_sdr_path)
fields_and_rasters = [
('usle_tot', usle_path), ('sed_export', sed_export_path),
('sed_dep', sed_deposition_path), ('avoid_exp', avoided_export_path),
('avoid_eros', avoided_erosion_path)]

# Using the list option for raster path bands so that we can reduce
# rasterizations, which are costly on large datasets.
zonal_stats_results = pygeoprocessing.zonal_statistics(
[(raster_path, 1) for (_, raster_path) in fields_and_rasters],
watershed_results_sdr_path,
polygons_might_overlap=geometries_might_overlap)

field_summaries = {
'usle_tot': pygeoprocessing.zonal_statistics(
(usle_path, 1), watershed_results_sdr_path),
'sed_export': pygeoprocessing.zonal_statistics(
(sed_export_path, 1), watershed_results_sdr_path),
'sed_dep': pygeoprocessing.zonal_statistics(
(sed_deposition_path, 1), watershed_results_sdr_path),
'avoid_exp': pygeoprocessing.zonal_statistics(
(avoided_export_path, 1), watershed_results_sdr_path),
'avoid_eros': pygeoprocessing.zonal_statistics(
(avoided_erosion_path, 1), watershed_results_sdr_path),
}
field: stats for ((field, _), stats) in
zip(fields_and_rasters, zonal_stats_results)}

for field_name in field_summaries:
field_def = ogr.FieldDefn(field_name, ogr.OFTReal)
Expand Down
6 changes: 3 additions & 3 deletions src/natcap/invest/sdr/sdr_core.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -436,17 +436,17 @@ def calculate_sediment_deposition(
flow_dir_info = pygeoprocessing.get_raster_info(mfd_flow_direction_path)
n_cols, n_rows = flow_dir_info['raster_size']
cdef int mfd_nodata = 0
cdef stack[int] processing_stack
cdef stack[long] processing_stack
cdef float sdr_nodata = pygeoprocessing.get_raster_info(
sdr_path)['nodata'][0]
cdef float e_prime_nodata = pygeoprocessing.get_raster_info(
e_prime_path)['nodata'][0]
cdef long col_index, row_index, win_xsize, win_ysize, xoff, yoff
cdef long global_col, global_row, j, k
cdef unsigned long flat_index
cdef long flat_index
cdef long seed_col = 0
cdef long seed_row = 0
cdef long neighbor_row, neighbor_col
cdef long neighbor_row, neighbor_col, ds_neighbor_row, ds_neighbor_col
cdef int flow_val, neighbor_flow_val, ds_neighbor_flow_val
cdef int flow_weight, neighbor_flow_weight
cdef float flow_sum, neighbor_flow_sum
Expand Down
139 changes: 78 additions & 61 deletions src/natcap/invest/urban_nature_access.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,15 +358,23 @@
"about": (
"The total population within the "
"administrative unit that is undersupplied "
"with urban nature.")
"with urban nature. If aggregating by "
"population groups, this will be the sum "
"of undersupplied populations across all "
"population groups within this administrative "
"unit.")
},
"Povr_adm": {
"type": "number",
"units": u.people,
"about": (
"The total population within the "
"administrative unit that is oversupplied "
"with urban nature.")
"with urban nature. If aggregating by "
"population groups, this will be the sum "
"of oversupplied populations across all "
"population groups within this administrative "
"unit.")
},
"SUP_DEMadm_cap_[POP_GROUP]": {
"type": "number",
Expand Down Expand Up @@ -1318,20 +1326,15 @@ def decay_func(dist_array):
output_dir,
f'urban_nature_balance_percapita_{pop_group}{suffix}.tif')
per_cap_urban_nature_balance_pop_group_task = graph.add_task(
pygeoprocessing.raster_calculator,
_calculate_urban_nature_balance_percapita,
kwargs={
'base_raster_path_band_const_list': [
(urban_nature_supply_percapita_to_group_path, 1),
(float(args['urban_nature_demand']), 'raw')
],
'local_op': _urban_nature_balance_percapita_op,
'target_raster_path':
per_cap_urban_nature_balance_pop_group_path,
'datatype_target': gdal.GDT_Float32,
'nodata_target': FLOAT32_NODATA
},
'urban_nature_supply_path':
urban_nature_supply_percapita_to_group_path,
'urban_nature_demand': float(args['urban_nature_demand']),
'target_path':
per_cap_urban_nature_balance_pop_group_path},
task_name=(
f'Calculate per-capita urban nature balance - {pop_group}'),
f'Calculate per-capita urban nature balance-{pop_group}'),
target_path_list=[
per_cap_urban_nature_balance_pop_group_path],
dependent_task_list=[
Expand Down Expand Up @@ -1409,20 +1412,17 @@ def decay_func(dist_array):
])

per_capita_urban_nature_balance_task = graph.add_task(
pygeoprocessing.raster_calculator,
_calculate_urban_nature_balance_percapita,
kwargs={
'base_raster_path_band_const_list': [
(file_registry['urban_nature_supply_percapita'], 1),
(float(args['urban_nature_demand']), 'raw')
],
'local_op': _urban_nature_balance_percapita_op,
'target_raster_path':
file_registry['urban_nature_balance_percapita'],
'datatype_target': gdal.GDT_Float32,
'nodata_target': FLOAT32_NODATA
},
task_name='Calculate per-capita urban nature balance',
target_path_list=[file_registry['urban_nature_balance_percapita']],
'urban_nature_supply_path':
file_registry['urban_nature_supply_percapita'],
'urban_nature_demand': float(args['urban_nature_demand']),
'target_path':
file_registry['urban_nature_balance_percapita']},
task_name=(
'Calculate per-capita urban nature balance}'),
target_path_list=[
file_registry['urban_nature_balance_percapita']],
dependent_task_list=[
urban_nature_supply_percapita_task,
])
Expand Down Expand Up @@ -1469,20 +1469,17 @@ def decay_func(dist_array):
RADIUS_OPT_URBAN_NATURE):
# This is "SUP_DEMi_cap" from the user's guide
per_capita_urban_nature_balance_task = graph.add_task(
pygeoprocessing.raster_calculator,
_calculate_urban_nature_balance_percapita,
kwargs={
'base_raster_path_band_const_list': [
(file_registry['urban_nature_supply_percapita'], 1),
(float(args['urban_nature_demand']), 'raw')
],
'local_op': _urban_nature_balance_percapita_op,
'target_raster_path':
file_registry['urban_nature_balance_percapita'],
'datatype_target': gdal.GDT_Float32,
'nodata_target': FLOAT32_NODATA
},
task_name='Calculate per-capita urban nature balance',
target_path_list=[file_registry['urban_nature_balance_percapita']],
'urban_nature_supply_path':
file_registry['urban_nature_supply_percapita'],
'urban_nature_demand': float(args['urban_nature_demand']),
'target_path':
file_registry['urban_nature_balance_percapita']},
task_name=(
'Calculate per-capita urban nature balance'),
target_path_list=[
file_registry['urban_nature_balance_percapita']],
dependent_task_list=[
urban_nature_supply_percapita_task,
])
Expand Down Expand Up @@ -1976,14 +1973,18 @@ def _get_zonal_stats(raster_path):
feature_id]['sum']
group_sup_dem_in_region = urban_nature_sup_dem_stats[
feature_id]['sum']
group_oversupply_in_region = oversupply_stats[feature_id]['sum']
group_undersupply_in_region = undersupply_stats[feature_id]['sum']
stats_by_feature[feature_id][f'SUP_DEMadm_cap_{groupname}'] = (
group_sup_dem_in_region / group_population_in_region)
stats_by_feature[feature_id][f'Pund_adm_{groupname}'] = (
undersupply_stats[feature_id]['sum'])
group_undersupply_in_region)
stats_by_feature[feature_id][f'Povr_adm_{groupname}'] = (
oversupply_stats[feature_id]['sum'])
group_oversupply_in_region)
sums['supply-demand'][feature_id] += group_sup_dem_in_region
sums['population'][feature_id] += group_population_in_region
sums['oversupply'][feature_id] += group_oversupply_in_region
sums['undersupply'][feature_id] += group_undersupply_in_region

for feature_id in feature_ids:
stats_by_feature[feature_id]['SUP_DEMadm_cap'] = (
Expand Down Expand Up @@ -2122,28 +2123,44 @@ def _write_supply_demand_vector(source_aoi_vector_path, feature_attrs,
target_vector = None


def _urban_nature_balance_percapita_op(urban_nature_supply, urban_nature_demand):
"""Calculate the per-capita urban nature balance.
def _calculate_urban_nature_balance_percapita(
urban_nature_supply_path, urban_nature_demand, target_path):
supply_nodata = pygeoprocessing.get_raster_info(
urban_nature_supply_path)['nodata'][0]

This is the amount of urban nature that each pixel has above (positive
values) or below (negative values) the user-defined ``urban_nature_demand``
value.
def _urban_nature_balance_percapita_op(urban_nature_supply,
urban_nature_demand):
"""Calculate the per-capita urban nature balance.
Args:
urban_nature_supply (numpy.array): The supply of urban nature 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.
urban_nature_demand (float): The policy-defined urban nature requirement,
in square meters per person.
This is the amount of urban nature that each pixel has above (positive
values) or below (negative values) the user-defined
``urban_nature_demand`` value.
Returns:
A ``numpy.array`` of the calculated urban nature budget.
"""
balance = numpy.full(
urban_nature_supply.shape, FLOAT32_NODATA, dtype=numpy.float32)
valid_pixels = ~numpy.isclose(urban_nature_supply, FLOAT32_NODATA)
balance[valid_pixels] = urban_nature_supply[valid_pixels] - urban_nature_demand
return balance
Args:
urban_nature_supply (numpy.array): The supply of urban nature
available to each person in the population. This is ``Ai`` in
the User's Guide.
urban_nature_demand (float): The policy-defined urban nature
requirement, in square meters per person.
Returns:
A ``numpy.array`` of the calculated urban nature budget.
"""
balance = numpy.full(
urban_nature_supply.shape, FLOAT32_NODATA, dtype=numpy.float32)
valid_pixels = ~pygeoprocessing.array_equals_nodata(
urban_nature_supply, supply_nodata)
balance[valid_pixels] = (
urban_nature_supply[valid_pixels] - urban_nature_demand)
return balance

pygeoprocessing.raster_calculator(
base_raster_path_band_const_list=[
(urban_nature_supply_path, 1), (urban_nature_demand, 'raw')],
local_op=_urban_nature_balance_percapita_op,
target_raster_path=target_path,
datatype_target=gdal.GDT_Float32,
nodata_target=FLOAT32_NODATA)


def _urban_nature_balance_totalpop_op(urban_nature_balance, population):
Expand Down
25 changes: 15 additions & 10 deletions tests/test_urban_nature_access.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,14 +217,19 @@ def test_urban_nature_balance(self):
[nodata, 100.5],
[75, 100]], dtype=numpy.float32)
urban_nature_demand = 50
supply_path = os.path.join(self.workspace_dir, 'supply.path')
target_path = os.path.join(self.workspace_dir, 'target.path')

population = numpy.array([
[50, 100],
[40.75, nodata]], dtype=numpy.float32)
pygeoprocessing.numpy_array_to_raster(
urban_nature_supply_percapita, nodata, _DEFAULT_PIXEL_SIZE,
_DEFAULT_ORIGIN, _DEFAULT_SRS.ExportToWkt(), supply_path)

urban_nature_access._calculate_urban_nature_balance_percapita(
supply_path, urban_nature_demand, target_path)

urban_nature_budget = pygeoprocessing.raster_to_numpy_array(
target_path)

urban_nature_budget = (
urban_nature_access._urban_nature_balance_percapita_op(
urban_nature_supply_percapita, urban_nature_demand))
expected_urban_nature_budget = numpy.array([
[nodata, 50.5],
[25, 50]], dtype=numpy.float32)
Expand Down Expand Up @@ -569,10 +574,10 @@ def test_radii_by_pop_group(self):
'pop_female': attributes[0]['pop_female'],
'pop_male': attributes[0]['pop_male'],
'adm_unit_id': 0,
'Pund_adm': 0,
'Pund_adm': 3991.8271484375,
'Pund_adm_female': 2235.423095703125,
'Pund_adm_male': 1756.404052734375,
'Povr_adm': 0,
'Povr_adm': 1084.1727294921875,
'Povr_adm_female': 607.13671875,
'Povr_adm_male': 477.0360107421875,
'SUP_DEMadm_cap': -17.907799109781322,
Expand Down Expand Up @@ -645,10 +650,10 @@ def test_radii_by_pop_group_exponential_kernal(self):
'pop_female': attributes[0]['pop_female'],
'pop_male': attributes[0]['pop_male'],
'adm_unit_id': 0,
'Pund_adm': 0,
'Pund_adm': 4801.7900390625,
'Pund_adm_female': 2689.00244140625,
'Pund_adm_male': 2112.78759765625,
'Povr_adm': 0,
'Povr_adm': 274.2098693847656,
'Povr_adm_female': 153.55752563476562,
'Povr_adm_male': 120.65234375,
'SUP_DEMadm_cap': -17.907799109781322,
Expand Down

0 comments on commit b5a98a4

Please sign in to comment.