VeniVidiVici
/
invest
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

clean up

This commit is contained in:
Emily Soth 2023-11-01 13:45:15 -07:00
parent ec6e3b8118
commit 0f24212f0a
20 changed files with 151 additions and 304 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/natcap/invest/annual_water_yield.py
View File

@ -735,7 +735,7 @@ def execute(args):
calculate_pet_task = graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=pet_op,
op=numpy.multiply, # PET = ET0 * KC
rasters=[eto_path, tmp_Kc_raster_path],
target_path=tmp_pet_path,
target_dtype=numpy.float32,
@ -781,7 +781,7 @@ def execute(args):
calculate_aet_task = graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=aet_op,
op=numpy.multiply, # AET = fractp * precip
rasters=[fractp_path, precip_path],
target_path=aet_path,
target_nodata=nodata_dict['out_nodata'],
@ -870,12 +870,7 @@ def execute(args):
graph.join()
def pet_op(eto, kc): return eto * kc
def aet_op(fractp, precip): return fractp * precip
# wyield equation to pass to raster_map
def wyield_op(fractp, precip): return (1 - fractp) * precip

11
src/natcap/invest/carbon.py
View File

@ -451,15 +451,13 @@ def execute(args):
continue
output_key = 'delta_cur_' + scenario_type
LOGGER.info("Calculate sequestration scenario '%s'", output_key)
storage_path_list = [
file_registry['tot_c_cur'],
file_registry['tot_c_' + scenario_type]]
diff_rasters_task = graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=diff_op,
rasters=storage_path_list,
op=numpy.subtract, # delta = scenario C - current C
rasters=[file_registry['tot_c_' + scenario_type],
file_registry['tot_c_cur']],
target_path=file_registry[output_key],
target_dtype=numpy.float32,
target_nodata=_CARBON_NODATA),
@ -518,9 +516,10 @@ def execute(args):
"Can't remove temporary file: %s\nOriginal Exception:\n%s",
file_registry[tmp_filename_key], os_error)
# element-wise sum function to pass to raster_map
def sum_op(*xs): return numpy.sum(xs, axis=0)
def diff_op(base, future): return future - base
def _accumulate_totals(raster_path):
"""Sum all non-nodata pixels in `raster_path` and return result."""

29
src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
View File

@ -1514,16 +1514,11 @@ def _calculate_npv(
prices_by_year[year] / (
(1 + discount_rate) ** years_since_baseline))
def _npv_op(*sequestration_arrays):
return numpy.sum(sequestration_arrays, axis=0) * valuation_factor
raster_paths = [
path for (_, path) in net_sequestration_rasters.items() if
year <= target_raster_year]
pygeoprocessing.raster_map(
op=_npv_op,
rasters=raster_paths,
op=lambda *seq_arrays: numpy.sum(seq, axis=0) * valuation_factor,
rasters=[
path for year, path in net_sequestration_rasters.items() if
year <= target_raster_year],
target_path=target_raster_path,
target_dtype=numpy.float32)
@ -1552,11 +1547,8 @@ def _calculate_stocks_after_baseline_period(
``None``.
"""
def _calculate_accumulation_over_years(baseline_matrix, accum_matrix):
return baseline_matrix + (accum_matrix * n_years)
pygeoprocessing.raster_map(
op=_calculate_accumulation_over_years,
op=lambda baseline, accum: baseline + (accum * n_years),
rasters=[baseline_stock_raster_path, yearly_accumulation_raster_path],
target_path=target_raster_path,
target_dtype=numpy.float32)
@ -2151,16 +2143,9 @@ def _reclassify_disturbance_magnitude(
``None``
"""
def _reclassify_disturbance(
landuse_transition_from_matrix, landuse_transition_to_matrix):
"""Pygeoprocessing op to reclassify disturbances."""
return disturbance_magnitude_matrix[
landuse_transition_from_matrix,
landuse_transition_to_matrix].toarray().flatten()
pygeoprocessing.raster_map(
op=_reclassify_disturbance,
op=lambda _from, _to: (
disturbance_magnitude_matrix[_from, _to].toarray().flatten()),
rasters=[landuse_transition_from_raster, landuse_transition_to_raster],
target_path=target_raster_path,
target_dtype=numpy.float32)

2
src/natcap/invest/coastal_vulnerability.py
View File

@ -3137,7 +3137,7 @@ def _aggregate_raster_values_in_radius(
projection matching base_point_vector_path.
sample_distance (float): radius around each point to search
for valid pixels.
kernel_path (string): path to create the kernel raster.
work_dir (string): path to directory for intermediate files.
target_pickle_path (string): path to pickle file storing dict
keyed by point id.
aggregation_op (function): takes a signal array and a mask array

5
src/natcap/invest/crop_production_regression.py
View File

@ -869,9 +869,8 @@ def _x_yield_op(
return result
def _min_op(y_n, y_p, y_k):
"""Calculate the min of the three inputs and multiply by Ymax."""
return numpy.min([y_n, y_k, y_p], axis=0)
"""equation for raster_map: calculate min of inputs and multiply by Ymax."""
def _min_op(y_n, y_p, y_k): return numpy.min([y_n, y_k, y_p], axis=0)
def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):

109
src/natcap/invest/ndr/ndr.py
View File

@ -733,7 +733,7 @@ def execute(args):
s_bar_task = task_graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=_bar_op,
op=numpy.divide, # s_bar = s_accum / flow_accum
rasters=[f_reg['s_accumulation_path'], f_reg['flow_accumulation_path']],
target_path=f_reg['s_bar_path'],
target_dtype=numpy.float32,
@ -890,7 +890,7 @@ def execute(args):
surface_export_task = task_graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=_calculate_export_op,
op=numpy.multiply, # export = load * ndr
rasters=[surface_load_path, ndr_path],
target_path=surface_export_path,
target_dtype=numpy.float32,
@ -930,7 +930,7 @@ def execute(args):
subsurface_export_task = task_graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=_calculate_export_op,
op=numpy.multiply, # export = load * ndr
rasters=[f_reg['sub_load_n_path'], f_reg['sub_ndr_n_path']],
target_path=f_reg['n_subsurface_export_path'],
target_dtype=numpy.float32,
@ -1030,6 +1030,23 @@ def execute(args):
LOGGER.info(r' (_") (_/(__)_) (__) (__) ')
# raster_map equation: Multiply a series of arrays element-wise
def _mult_op(*array_list): return numpy.prod(numpy.stack(array_list), axis=0)
# raster_map equation: Sum a list of arrays element-wise
def _sum_op(*array_list): return numpy.sum(array_list, axis=0)
# raster_map equation: calculate inverse of S factor
def _inverse_op(base_val): return numpy.where(base_val == 0, 0, 1 / base_val)
# raster_map equation: rescale and threshold slope between 0.005 and 1
def _slope_proportion_and_threshold_op(slope):
slope_fraction = slope / 100
slope_fraction[slope_fraction < 0.005] = 0.005
slope_fraction[slope_fraction > 1] = 1
return slope_fraction
def _create_mask_raster(source_raster_path, source_vector_path,
target_raster_path):
"""Create a mask raster from a vector.
@ -1095,14 +1112,6 @@ def _mask_raster(source_raster_path, mask_raster_path,
target_masked_raster_path, target_dtype, nodata)
def _slope_proportion_and_threshold_op(slope):
"""Rescale and threshold slope between 0.005 and 1.0."""
slope_fraction = slope / 100
slope_fraction[slope_fraction < 0.005] = 0.005
slope_fraction[slope_fraction > 1] = 1
return slope_fraction
def _add_fields_to_shapefile(field_pickle_map, target_vector_path):
"""Add fields and values to an OGR layer open for writing.
@ -1217,14 +1226,8 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
if value_count > 0:
value_mean /= value_count
def _normalize_raster_op(array):
"""Divide values by mean."""
if value_mean == 0:
return array
return array / value_mean
pygeoprocessing.raster_map(
op=_normalize_raster_op,
op=lambda array: array if value_mean == 0 else array / value_mean,
rasters=[base_raster_path_band[0]],
target_path=target_normalized_raster_path,
target_dtype=numpy.float32)
@ -1268,14 +1271,6 @@ def _calculate_load(lulc_raster_path, lucode_to_load, target_load_raster):
target_nodata=_TARGET_NODATA)
#Multiply a series of arrays element-wise
def _mult_op(*array_list): return numpy.prod(numpy.stack(array_list), axis=0)
# Sum a list of arrays element-wise
def _sum_op(*array_list): return numpy.sum(array_list, axis=0)
def _map_surface_load(
modified_load_path, lulc_raster_path, lucode_to_subsurface_proportion,
target_surface_load_path):
@ -1334,14 +1329,10 @@ def _map_subsurface_load(
keys = sorted(numpy.array(list(proportion_subsurface_map)))
subsurface_permeance_values = numpy.array(
[proportion_subsurface_map[x] for x in keys])
def _map_subsurface_load_op(lucode_array, modified_load_array):
"""Convert unit load to total load & handle nodata."""
index = numpy.digitize(lucode_array.ravel(), keys, right=True)
return modified_load_array * subsurface_permeance_values[index]
pygeoprocessing.raster_map(
op=_map_subsurface_load_op,
op=lambda lulc, modified_load: (
modified_load * subsurface_permeance_values[
numpy.digitize(lulc.ravel(), keys, right=True)]),
rasters=[lulc_raster_path, modified_load_path],
target_path=target_sub_load_path,
target_dtype=numpy.float32,
@ -1366,47 +1357,29 @@ def _map_lulc_to_val_mask_stream(
"""
lucodes = sorted(lucodes_to_vals.keys())
values = numpy.array([lucodes_to_vals[x] for x in lucodes])
def _map_eff_op(lucode_array, stream_array):
"""Map efficiency from LULC and handle nodata/streams."""
index = numpy.digitize(lucode_array.ravel(), lucodes, right=True)
return values[index] * (1 - stream_array)
pygeoprocessing.raster_map(
op=_map_eff_op,
op=lambda lulc, stream: (
values[numpy.digitize(lulc.ravel(), lucodes, right=True)] *
(1 - stream)),
rasters=[lulc_raster_path, stream_path],
target_path=target_eff_path,
target_dtype=numpy.float32,
target_nodata=_TARGET_NODATA)
def _bar_op(s_accumulation, flow_accumulation):
"""Calculate bar operation of s_accum / flow_accum."""
return s_accumulation / flow_accumulation
def d_up_calculation(s_bar_path, flow_accum_path, target_d_up_path):
"""Calculate d_up = s_bar * sqrt(upslope area)."""
cell_area_m2 = abs(numpy.prod(pygeoprocessing.get_raster_info(
s_bar_path)['pixel_size']))
def _d_up_op(s_bar, flow_accumulation):
"""Calculate d_up index."""
return s_bar * numpy.sqrt(flow_accumulation * cell_area_m2)
pygeoprocessing.raster_map(
op=_d_up_op,
op=lambda s_bar, flow_accum: (
s_bar * numpy.sqrt(flow_accum * cell_area_m2)),
rasters=[s_bar_path, flow_accum_path],
target_path=target_d_up_path,
target_dtype=numpy.float32,
target_nodata=_TARGET_NODATA)
def _inverse_op(base_val):
"""Calculate inverse of S factor."""
return numpy.where(base_val == 0, 0, 1 / base_val)
def calculate_ic(d_up_path, d_dn_path, target_ic_path):
"""Calculate IC as log_10(d_up/d_dn)."""
ic_nodata = float(numpy.finfo(numpy.float32).min)
@ -1439,14 +1412,9 @@ def _calculate_ndr(
ic_factor_raster = None
ic_0_param = (ic_min + ic_max) / 2
def _calculate_ndr_op(effective_retention_array, ic_array):
"""Calculate NDR."""
return (
(1 - effective_retention_array) /
(1 + numpy.exp((ic_0_param - ic_array) / k_param)))
pygeoprocessing.raster_map(
op=_calculate_ndr_op,
op=lambda eff, ic: ((1 - eff) /
(1 + numpy.exp((ic_0_param - ic) / k_param))),
rasters=[effective_retention_path, ic_factor_path],
target_path=target_ndr_path,
target_dtype=numpy.float32,
@ -1456,25 +1424,16 @@ def _calculate_ndr(
def _calculate_sub_ndr(
eff_sub, crit_len_sub, dist_to_channel_path, target_sub_ndr_path):
"""Calculate subsurface: subndr = eff_sub(1-e^(-5*l/crit_len)."""
def _sub_ndr_op(dist_to_channel_array):
"""Calculate subsurface NDR."""
return 1 - eff_sub * (
1 - numpy.exp(-5 * dist_to_channel_array / crit_len_sub))
pygeoprocessing.raster_map(
op=_sub_ndr_op,
op=lambda dist_to_channel: (
1 - eff_sub *
(1 - numpy.exp(-5 * dist_to_channel / crit_len_sub))),
rasters=[dist_to_channel_path],
target_path=target_sub_ndr_path,
target_dtype=numpy.float32,
target_nodata=_TARGET_NODATA)
def _calculate_export_op(load_array, ndr_array):
"""Multiply load by NDR."""
return load_array * ndr_array
def _aggregate_and_pickle_total(
base_raster_path_band, aggregate_vector_path, target_pickle_path):
"""Aggregate base raster path to vector path FIDs and pickle result.

73
src/natcap/invest/pollination.py
View File

@ -1067,6 +1067,25 @@ def execute(args):
task_graph.join()
def pollinator_supply_op(
foraged_flowers_array, floral_resources_array,
convolve_ps_array):
"""raster_map equation: calculate (RA*fa)/FR * convolve(PS)."""
return numpy.where(
floral_resources_array == 0, 0,
foraged_flowers_array / floral_resources_array * convolve_ps_array)
def on_farm_pollinator_abundance_op(h, pat):
"""raster_map equation: return (pat * (1 - h)) / (h * (1 - 2*pat)+pat))"""
return (pat * (1 - h)) / (h * (1 - 2 * pat) + pat)
# raster_map equation: return min(mp_array+FP_array, 1)
def pyt_op(mp_array, FP_array): return numpy.minimum(mp_array + FP_array, 1)
# raster_map equation: return max(0,PYT_array-mp_array)
def pyw_op(mp_array, PYT_array): return numpy.maximum(PYT_array - mp_array, 0)
def _rasterize_vector_onto_base(
base_raster_path, base_vector_path, attribute_id,
target_raster_path, filter_string=None):
@ -1366,6 +1385,19 @@ def _parse_scenario_variables(args):
return result
def _sum_arrays(*array_list):
"""Calculate sum of array_list and account for nodata."""
valid_mask = numpy.zeros(array_list[0].shape, dtype=bool)
result = numpy.empty_like(array_list[0])
result[:] = 0
for array in array_list:
local_valid_mask = ~pygeoprocessing.array_equals_nodata(array, _INDEX_NODATA)
result[local_valid_mask] += array[local_valid_mask]
valid_mask |= local_valid_mask
result[~valid_mask] = _INDEX_NODATA
return result
def _calculate_habitat_nesting_index(
substrate_path_map, species_substrate_index_map,
target_habitat_nesting_index_path):
@ -1394,35 +1426,13 @@ def _calculate_habitat_nesting_index(
def _multiply_by_scalar(raster_path, scalar, target_path):
"""Multiply a raster by a scalar and write out result."""
pygeoprocessing.raster_map(
op=lambda array: array * scalar,
rasters=[raster_path],
target_path=target_path)
def _sum_arrays(*array_list):
"""Calculate sum of array_list and account for nodata."""
valid_mask = numpy.zeros(array_list[0].shape, dtype=bool)
result = numpy.empty_like(array_list[0])
result[:] = 0
for array in array_list:
local_valid_mask = ~pygeoprocessing.array_equals_nodata(array, _INDEX_NODATA)
result[local_valid_mask] += array[local_valid_mask]
valid_mask |= local_valid_mask
result[~valid_mask] = _INDEX_NODATA
return result
def pollinator_supply_op(
foraged_flowers_array, floral_resources_array,
convolve_ps_array):
"""Calculating (RA*fa)/FR * convolve(PS)."""
return numpy.where(
floral_resources_array == 0, 0,
foraged_flowers_array / floral_resources_array * convolve_ps_array)
def _calculate_pollinator_supply_index(
habitat_nesting_suitability_path, floral_resources_path,
species_abundance, target_path):
@ -1441,23 +1451,6 @@ def _calculate_pollinator_supply_index(
target_path=target_path)
def on_farm_pollinator_abundance_op(h_array, pat_array):
"""Return (pad * (1 - h)) / (h * (1 - 2*pat)+pat)) tolerate nodata."""
return (
(pat_array * (1 - h_array)) /
(h_array * (1 - 2 * pat_array) + pat_array))
def pyt_op(mp_array, FP_array):
"""Return min(mp_array+FP_array, 1)."""
return numpy.minimum(mp_array + FP_array, 1)
def pyw_op(mp_array, PYT_array):
"""Return max(0,PYT_array-mp_array)."""
return numpy.maximum(PYT_array - mp_array, 0)
@validation.invest_validator
def validate(args, limit_to=None):
"""Validate args to ensure they conform to `execute`'s contract.

8
src/natcap/invest/scenic_quality/scenic_quality.py
View File

@ -1081,13 +1081,9 @@ def _calculate_visual_quality(source_raster_path, working_dir, target_path):
percentile_bins = numpy.array(percentile_values)
LOGGER.info('Mapping percentile breaks %s', percentile_bins)
def _map_percentiles(valuation_matrix):
return numpy.where(
valuation_matrix == 0, 0, numpy.digitize(
valuation_matrix, percentile_bins))
pygeoprocessing.raster_map(
op=_map_percentiles,
op=lambda val_array: numpy.where(
val_array == 0, 0, numpy.digitize(val_array, percentile_bins)),
rasters=[source_raster_path],
target_path=target_path,
target_dtype=numpy.uint8)

155
src/natcap/invest/sdr/sdr.py
View File

@ -583,7 +583,7 @@ def execute(args):
threshold_slope_task = task_graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=threshold_slope,
op=threshold_slope_op,
rasters=[f_reg['slope_path']],
target_path=f_reg['thresholded_slope_path']),
target_path_list=[f_reg['thresholded_slope_path']],
@ -769,9 +769,11 @@ def execute(args):
task_name='calculate sdr')
sed_export_task = task_graph.add_task(
func=_calculate_sed_export,
args=(
f_reg['usle_path'], f_reg['sdr_path'], f_reg['sed_export_path']),
func=pygeoprocessing.raster_map,
kwargs=dict(
op=numpy.multiply, # export = USLE * SDR
rasters=[f_reg['usle_path'], f_reg['sdr_path']],
target_path=f_reg['sed_export_path']),
target_path_list=[f_reg['sed_export_path']],
dependent_task_list=[usle_task, sdr_task],
task_name='calculate sed export')
@ -798,7 +800,7 @@ def execute(args):
avoided_erosion_task = task_graph.add_task(
func=pygeoprocessing.raster_map,
kwargs=dict(
op=_avoided_erosion_op,
op=numpy.subtract, # avoided erosion = rkls - usle
rasters=[f_reg['rkls_path'], f_reg['usle_path']],
target_path=f_reg['avoided_erosion_path']),
dependent_task_list=[rkls_task, usle_task],
@ -844,7 +846,7 @@ def execute(args):
def _avoided_export_op(avoided_erosion, sdr, sed_deposition):
"""Calculate total retention.
"""raster_map equation: calculate total retention.
Args:
avoided_erosion (numpy.array): Avoided erosion values.
@ -859,19 +861,25 @@ def _avoided_export_op(avoided_erosion, sdr, sed_deposition):
# known as a modified USLE)
return avoided_erosion * sdr + sed_deposition
def _avoided_erosion_op(rkls, usle):
"""Calculate avoided erosion.
def add_drainage_op(stream, drainage):
"""raster_map equation: add drainage mask to stream layer.
Args:
rkls (numpy.array): Computed RKLS values.
usle (numpy.array): Computed USLE values.
stream (numpy.array): binary array where 1 indicates
a stream, and 0 is a valid landscape pixel but not a stream.
drainage (numpy.array): binary array where 1 indicates any water
reaching that pixel drains to a stream.
Returns:
A ``numpy.array`` of the same size and shape of the input numpy
arrays containing computed avoided erosion values.
numpy.array combination of stream and drainage
"""
return rkls - usle
return numpy.where(drainage == 1, 1, stream)
# raster_map equation: calculate USLE
def usle_op(rkls, cp_factor): return rkls * cp_factor
# raster_map equation: calculate the inverse ws factor
def inverse_ws_op(w_factor, s_factor): return 1 / (w_factor * s_factor)
def _calculate_what_drains_to_stream(
@ -1003,15 +1011,11 @@ def _calculate_ls_factor(
None
"""
slope_nodata = pygeoprocessing.get_raster_info(slope_path)['nodata'][0]
flow_accumulation_info = pygeoprocessing.get_raster_info(
flow_accumulation_path)
flow_accumulation_nodata = flow_accumulation_info['nodata'][0]
cell_size = abs(flow_accumulation_info['pixel_size'][0])
cell_size = abs(pygeoprocessing.get_raster_info(
flow_accumulation_path)['pixel_size'][0])
cell_area = cell_size ** 2
def ls_factor_function(percent_slope, flow_accumulation, l_max):
def ls_factor_function(percent_slope, flow_accumulation):
"""Calculate the LS factor.
Args:
@ -1023,15 +1027,6 @@ def _calculate_ls_factor(
ls_factor
"""
# avg aspect intermediate output should always have a defined
# nodata value from pygeoprocessing
valid_mask = (
~pygeoprocessing.array_equals_nodata(percent_slope, slope_nodata) &
~pygeoprocessing.array_equals_nodata(
flow_accumulation, flow_accumulation_nodata))
result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
result[:] = _TARGET_NODATA
# Although Desmet & Govers (1996) discusses "upstream contributing
# area", this is not strictly defined. We decided to use the square
# root of the upstream contributing area here as an estimate, which
@ -1040,9 +1035,8 @@ def _calculate_ls_factor(
# We subtract 1 from the flow accumulation because FA includes itself
# in its count of pixels upstream and our LS factor equation wants only
# those pixels that are strictly upstream.
contributing_area = numpy.sqrt(
(flow_accumulation[valid_mask]-1) * cell_area)
slope_in_radians = numpy.arctan(percent_slope[valid_mask] / 100.0)
contributing_area = numpy.sqrt((flow_accumulation - 1) * cell_area)
slope_in_radians = numpy.arctan(percent_slope / 100)
aspect_length = (numpy.fabs(numpy.sin(slope_in_radians)) +
numpy.fabs(numpy.cos(slope_in_radians)))
@ -1050,7 +1044,7 @@ def _calculate_ls_factor(
# From Equation 4 in "Extension and validation of a geographic
# information system ..."
slope_factor = numpy.where(
percent_slope[valid_mask] < 9.0,
percent_slope < 9,
10.8 * numpy.sin(slope_in_radians) + 0.03,
16.8 * numpy.sin(slope_in_radians) - 0.5)
@ -1064,10 +1058,9 @@ def _calculate_ls_factor(
slope_table = numpy.array([1., 3.5, 5., 9.])
m_table = numpy.array([0.2, 0.3, 0.4, 0.5])
# mask where slopes are larger than lookup table
big_slope_mask = percent_slope[valid_mask] > slope_table[-1]
big_slope_mask = percent_slope > slope_table[-1]
m_indexes = numpy.digitize(
percent_slope[valid_mask][~big_slope_mask], slope_table,
right=True)
percent_slope[~big_slope_mask], slope_table, right=True)
m_exp = numpy.empty(big_slope_mask.shape, dtype=numpy.float32)
m_exp[big_slope_mask] = (
beta[big_slope_mask] / (1 + beta[big_slope_mask]))
@ -1082,14 +1075,12 @@ def _calculate_ls_factor(
# threshold L factor to l_max
l_factor[l_factor > l_max] = l_max
result[valid_mask] = l_factor * slope_factor
return result
return l_factor * slope_factor
pygeoprocessing.raster_calculator(
[(path, 1) for path in [slope_path, flow_accumulation_path]] + [
(l_max, 'raw')],
ls_factor_function, target_ls_factor_path, gdal.GDT_Float32,
_TARGET_NODATA)
pygeoprocessing.raster_map(
op=ls_factor_function,
rasters=[slope_path, flow_accumulation_path],
target_path=target_ls_factor_path)
def _calculate_rkls(
@ -1168,8 +1159,8 @@ def _calculate_rkls(
rkls_function, rkls_path, gdal.GDT_Float32, _TARGET_NODATA)
def threshold_slope(slope):
"""Convert slope to m/m and clamp at 0.005 and 1.0.
def threshold_slope_op(slope):
"""raster_map equation: convert slope to m/m and clamp at 0.005 and 1.0.
As desribed in Cavalli et al., 2013.
"""
@ -1179,21 +1170,6 @@ def threshold_slope(slope):
return slope_m
def add_drainage_op(stream, drainage):
"""Add drainage mask to stream layer.
Args:
stream (numpy.array): binary array where 1 indicates
a stream, and 0 is a valid landscape pixel but not a stream.
drainage (numpy.array): binary array where 1 indicates any water
reaching that pixel drains to a stream.
Returns:
numpy.array combination of stream and drainage
"""
return numpy.where(drainage == 1, 1, stream)
def _calculate_w(
lulc_to_c, lulc_path, w_factor_path,
out_thresholded_w_factor_path):
@ -1227,12 +1203,8 @@ def _calculate_w(
(lulc_path, 1), lulc_to_c, w_factor_path, gdal.GDT_Float32,
_TARGET_NODATA, reclass_error_details)
def threshold_w(w_val):
"""Threshold w to 0.001."""
return numpy.where(w_val < 0.001, 0.001, w_val)
pygeoprocessing.raster_map(
op=threshold_w,
op=lambda w_val: numpy.where(w_val < 0.001, 0.001, w_val),
rasters=[w_factor_path],
target_path=out_thresholded_w_factor_path)
@ -1265,11 +1237,6 @@ def _calculate_cp(lulc_to_cp, lulc_path, cp_factor_path):
_TARGET_NODATA, reclass_error_details)
def usle_op(rkls, cp_factor):
"""Calculate USLE."""
return rkls * cp_factor
def _calculate_bar_factor(
flow_direction_path, factor_path, flow_accumulation_path,
accumulation_path, out_bar_path):
@ -1292,9 +1259,6 @@ def _calculate_bar_factor(
None.
"""
flow_accumulation_nodata = pygeoprocessing.get_raster_info(
flow_accumulation_path)['nodata'][0]
LOGGER.debug("doing flow accumulation mfd on %s", factor_path)
# manually setting compression to DEFLATE because we got some LZW
# errors when testing with large data.
@ -1305,12 +1269,8 @@ def _calculate_bar_factor(
'TILED=YES', 'BIGTIFF=YES', 'COMPRESS=DEFLATE',
'PREDICTOR=3']))
def bar_op(base_accumulation, flow_accumulation):
"""Aggregate accumulation from base divided by the flow accum."""
return base_accumulation / flow_accumulation
pygeoprocessing.raster_map(
op=bar_op,
op=lambda base_accum, flow_accum: base_accum / flow_accum,
rasters=[accumulation_path, flow_accumulation_path],
target_path=out_bar_path)
@ -1320,17 +1280,9 @@ def _calculate_d_up(
"""Calculate w_bar * s_bar * sqrt(flow accumulation * cell area)."""
cell_area = abs(
pygeoprocessing.get_raster_info(w_bar_path)['pixel_size'][0])**2
def d_up_op(w_bar, s_bar, flow_accumulation):
"""Calculate the d_up index.
w_bar * s_bar * sqrt(upslope area)
"""
return w_bar * s_bar * numpy.sqrt(flow_accumulation * cell_area)
pygeoprocessing.raster_map(
op=d_up_op,
op=lambda w_bar, s_bar, flow_accum: (
w_bar * s_bar * numpy.sqrt(flow_accum * cell_area)),
rasters=[w_bar_path, s_bar_path, flow_accumulation_path],
target_path=out_d_up_path)
@ -1340,22 +1292,13 @@ def _calculate_d_up_bare(
"""Calculate s_bar * sqrt(flow accumulation * cell area)."""
cell_area = abs(
pygeoprocessing.get_raster_info(s_bar_path)['pixel_size'][0])**2
def d_up_op(s_bar, flow_accumulation):
"""Calculate the bare d_up index: s_bar * sqrt(upslope area)."""
return numpy.sqrt(flow_accumulation * cell_area) * s_bar
pygeoprocessing.raster_map(
op=d_up_op,
op=lambda s_bar, flow_accum: (
numpy.sqrt(flow_accum * cell_area) * s_bar),
rasters=[s_bar_path, flow_accumulation_path],
target_path=out_d_up_bare_path)
def inverse_ws_op(w_factor, s_factor):
"""Calculate the inverse ws factor."""
return 1 / (w_factor * s_factor)
def _calculate_ic(d_up_path, d_dn_path, out_ic_factor_path):
"""Calculate log10(d_up/d_dn)."""
# ic can be positive or negative, so float.min is a reasonable nodata value
@ -1397,18 +1340,6 @@ def _calculate_sdr(
gdal.GDT_Float32, _TARGET_NODATA)
def _calculate_sed_export(usle_path, sdr_path, target_sed_export_path):
"""Calculate USLE * SDR."""
def sed_export_op(usle, sdr):
"""Sediment export."""
return usle * sdr
pygeoprocessing.raster_map(
op=sed_export_op,
rasters=[usle_path, sdr_path],
target_path=target_sed_export_path)
def _calculate_e_prime(usle_path, sdr_path, stream_path, target_e_prime):
"""Calculate USLE * (1-SDR)."""
def e_prime_op(usle, sdr, streams):

29
src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
View File

@ -965,6 +965,18 @@ def execute(args):
LOGGER.info(' //_| //_|')
# raster_map equation: sum the monthly qfis
def qfi_sum_op(*qf_values): return numpy.sum(qf_values)
def _calculate_l_avail(l_path, gamma, target_l_avail_path):
"""l avail = l * gamma."""
pygeoprocessing.raster_map(
op=lambda l: numpy.min(numpy.stack((gamma * l, l)), axis=0),
rasters=[l_path],
target_path=target_l_avail_path)
def _calculate_vri(l_path, target_vri_path):
"""Calculate VRI as li_array / qb_sum.
@ -1006,10 +1018,6 @@ def _calculate_vri(l_path, target_vri_path):
li_nodata)
"""Sum the monthly qfis."""
def qfi_sum_op(*qf_values): return numpy.sum(qf_values)
def _calculate_monthly_quick_flow(precip_path, n_events_path, stream_path,
si_path, qf_monthly_path):
"""Calculate quick flow for a month.
@ -1361,19 +1369,6 @@ def _aggregate_recharge(
aggregate_vector = None
def _calculate_l_avail(l_path, gamma, target_l_avail_path):
"""l avail = l * gamma."""
def l_avail_op(l_array):
"""Calculate equation [8] L_avail = min(gamma*L, L)."""
return numpy.min(numpy.stack((gamma*l_array, l_array)), axis=0)
pygeoprocessing.raster_map(
op=l_avail_op,
rasters=[l_path],
target_path=target_l_avail_path)
@validation.invest_validator
def validate(args, limit_to=None):
"""Validate args to ensure they conform to `execute`'s contract.

2
tests/test_annual_water_yield.py
View File

@ -43,7 +43,7 @@ class AnnualWaterYieldTests(unittest.TestCase):
'biophysical_table_path': os.path.join(
SAMPLE_DATA, 'biophysical_table.csv'),
'seasonality_constant': 5,
'n_workers': 1,
'n_workers': -1,
}
return args

2
tests/test_carbon.py
View File

@ -114,7 +114,7 @@ class CarbonTests(unittest.TestCase):
'lulc_cur_year': 2016,
'lulc_fut_year': 2030,
'discount_rate': -7.1,
'n_workers': 1,
'n_workers': -1,
}
# Create LULC rasters and pools csv in workspace and add them to args.

1
tests/test_coastal_blue_carbon.py
View File

@ -508,6 +508,7 @@ class TestCBC2(unittest.TestCase):
'price_table_path': os.path.join(target_dir,
'price_table.csv'),
'discount_rate': 4,
'n_workers': 1
}
with open(args['price_table_path'], 'w') as price_table:

5
tests/test_crop_production.py
View File

@ -47,7 +47,7 @@ class CropProductionTests(unittest.TestCase):
'aggregate_polygon_path': os.path.join(
SAMPLE_DATA_PATH, 'aggregate_shape.shp'),
'model_data_path': MODEL_DATA_PATH,
'n_workers': 1
'n_workers': '-1'
}
crop_production_percentile.execute(args)
@ -222,7 +222,7 @@ class CropProductionTests(unittest.TestCase):
'nitrogen_fertilization_rate': 29.6,
'phosphorus_fertilization_rate': 8.4,
'potassium_fertilization_rate': 14.2,
'n_workers': 1
'n_workers': '-1'
}
with open(args['landcover_to_crop_table_path'],
@ -307,7 +307,6 @@ class CropProductionTests(unittest.TestCase):
'nitrogen_fertilization_rate': 29.6,
'phosphorus_fertilization_rate': 8.4,
'potassium_fertilization_rate': 14.2,
'n_workers': 1
}
crop_production_regression.execute(args)

2
tests/test_forest_carbon_edge.py
View File

@ -44,7 +44,7 @@ class ForestCarbonEdgeTests(unittest.TestCase):
'tropical_forest_edge_carbon_model_vector_path': os.path.join(
REGRESSION_DATA, 'input', 'core_data',
'forest_carbon_edge_regression_model_parameters.shp'),
'workspace_dir': '/Users/emily/Documents/fcee22',
'workspace_dir': self.workspace_dir,
'n_workers': -1
}
forest_carbon_edge_effect.execute(args)

1
tests/test_pollination.py
View File

@ -73,7 +73,6 @@ class PollinationTests(unittest.TestCase):
'landcover_biophysical_table_simple.csv'),
'farm_vector_path': os.path.join(
REGRESSION_DATA, 'input', 'blueberry_ridge_farm.shp'),
'n_workers': 2
}
# make empty result files to get coverage for removing if necessary
result_files = ['farm_results.shp', 'total_pollinator_yield.tif',

2
tests/test_scenario_proximity.py
View File

@ -38,7 +38,7 @@ class ScenarioProximityTests(unittest.TestCase):
'focal_landcover_codes': '1 2 3 4 5',
'n_fragmentation_steps': '1',
'replacement_lucode': '12',
'n_workers': 2,
'n_workers': '-1',
}
return args

1
tests/test_seasonal_water_yield_regression.py
View File

@ -615,7 +615,6 @@ class SeasonalWaterYieldRegressionTests(unittest.TestCase):
args = SeasonalWaterYieldRegressionTests.generate_base_args(
self.workspace_dir)
args['workspace_dir'] = '/Users/emily/Documents/swypass'
# Ensure the model can pass when a nodata value is not defined.
size = 100
lulc_array = numpy.zeros((size, size), dtype=numpy.int8)

4
tests/test_ucm.py
View File

@ -27,7 +27,7 @@ class UCMTests(unittest.TestCase):
"""UCM: regression: CC Factors."""
import natcap.invest.urban_cooling_model
args = {
'workspace_dir': '/Users/emily/Documents/testucmfail',#os.path.join(self.workspace_dir, 'workspace'),
'workspace_dir': os.path.join(self.workspace_dir, 'workspace'),
'results_suffix': 'test_suffix',
't_ref': 35.0,
't_obs_raster_path': os.path.join(
@ -208,7 +208,7 @@ class UCMTests(unittest.TestCase):
'avg_tmp_an': 1.608697970397692,
'avd_eng_cn': 7239992.744486,
'avg_wbgt_v': 31.91108630952381,
'avg_ltls_v': 28.73463901689708, # changed due to using float64 vs float32
'avg_ltls_v': 28.73463901689708,
'avg_hvls_v': 75.000000000000000,
}
try:

3
tests/test_urban_nature_access.py
View File

@ -819,14 +819,11 @@ class UNATests(unittest.TestCase):
target_path)
weighted_sum_array = pygeoprocessing.raster_to_numpy_array(target_path)
print(list(weighted_sum_array))
weighted_sum_nodata = pygeoprocessing.get_raster_info(
target_path)['nodata'][0]
print(weighted_sum_nodata)
# check that we have the expected number of nodata pixels
nodata_pixels = numpy.isclose(weighted_sum_array, weighted_sum_nodata)
print(list(nodata_pixels))
self.assertEqual(
numpy.count_nonzero(nodata_pixels), 2)