Updating precision notes and assertion.

This commit is contained in:
James Douglass 2021-05-18 08:59:56 -05:00
parent a0d8c4f483
commit ded9ac6cdb
2 changed files with 8 additions and 4 deletions

View File

@ -407,8 +407,8 @@ def _accumulate_totals(raster_path):
raster_sum = 0.0
for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
# The float64 dtype in the sum is needed to reduce numerical error in
# the sum. Users calculated the sum by hand, noticed a difference and
# wrote to us about it on the forum, so hence this fix.
# the sum. Users calculated the sum with ArcGIS zonal statistics,
# noticed a difference and wrote to us about it on the forum.
raster_sum += numpy.sum(
block[~numpy.isclose(block, nodata)], dtype=numpy.float64)
return raster_sum

View File

@ -344,8 +344,9 @@ class CarbonValidationTests(unittest.TestCase):
from natcap.invest import carbon
big_float32_array = numpy.random.default_rng(seed=1).random(
(10000, 10000), dtype=numpy.float32)
(1000, 1000), dtype=numpy.float32)
# Throw in some nodata values for good measure.
nodata = numpy.finfo(numpy.float32).min
big_float32_array[1:15] = nodata
@ -357,6 +358,9 @@ class CarbonValidationTests(unittest.TestCase):
big_float32_array, float(nodata), (2, -2), (2, -2), wkt,
raster_path)
# Verify better-than-float32 precision on raster summation.
# Using a numpy float32 in numpy.sum will pass up to rtol=1e-9.
numpy.testing.assert_allclose(
carbon._accumulate_totals(raster_path),
49935086.309929)
492919.73994,
rtol=1e-12) # Note better precision