Merge branch 'main' into bugfix/1509-una-validation-missing-uniform-search-radius
This commit is contained in:
James Douglass
GitHub
commit
85d8028aad
16
HISTORY.rst
16
HISTORY.rst
|
|
@ -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)
|
||||
-------------------
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@ from osgeo import ogr
|
|||
|
||||
from .. import gettext
|
||||
from .. import spec_utils
|
||||
from .. import urban_nature_access
|
||||
from .. import utils
|
||||
from .. import validation
|
||||
from ..model_metadata import MODEL_METADATA
|
||||
|
|
@ -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=(
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -360,7 +360,11 @@ MODEL_SPEC = {
|
|||
"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",
|
||||
|
|
@ -368,7 +372,11 @@ MODEL_SPEC = {
|
|||
"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",
|
||||
|
|
@ -1320,20 +1328,15 @@ def execute(args):
|
|||
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=[
|
||||
|
|
@ -1411,20 +1414,17 @@ def execute(args):
|
|||
])
|
||||
|
||||
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,
|
||||
])
|
||||
|
|
@ -1471,20 +1471,17 @@ def execute(args):
|
|||
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,
|
||||
])
|
||||
|
|
@ -1978,14 +1975,18 @@ def _supply_demand_vector_for_pop_groups(
|
|||
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'] = (
|
||||
|
|
@ -2124,28 +2125,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):
|
||||
|
|
|
|||
|
|
@ -217,14 +217,19 @@ class UNATests(unittest.TestCase):
|
|||
[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)
|
||||
|
|
@ -569,10 +574,10 @@ class UNATests(unittest.TestCase):
|
|||
'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,
|
||||
|
|
@ -645,10 +650,10 @@ class UNATests(unittest.TestCase):
|
|||
'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,
|
||||
|
|
|
|||
Loading…
Reference in New Issue