Skip to content

Maptotal: Wrong results for large maps  #336

@dk-forestry

Description

@dk-forestry

Hi,
maptotal seems to give wrong results when it is used for large maps (>1 Mio cells) on pcraster 4.3.1, Python 3.9 and Windows 10. The following code should reproduce this behaviour:

from pcraster import *

def do_tests(map_dimensions, map_values):
    results = {}

    for map_dimension in map_dimensions:
        setclone(map_dimension, map_dimension, 1, 0, 0)

        for map_value in map_values:
            map_of_value = ifthen(1, scalar(map_value))
            array_of_value = pcr2numpy(map_of_value, 0)  # convert to numpy array to make sure that the problem is not with ifthen()

            results[f'{map_dimension}x{map_value}'] = {
                'observed_from_maptotal': float(maptotal(map_of_value)),
                'observed_from_array': array_of_value.sum(),
                'expected': map_value * map_dimension ** 2,
            }

    return results


map_dimensions = [10, 100, 500, 1000, 2000, 3000, 4000, 5000]
map_values = [0.01, 1]

results = do_tests(map_dimensions, map_values)

for map_dimensions_x_map_value, map_value_dict in results.items():
    difference = map_value_dict['observed_from_maptotal'] / map_value_dict['expected']
    if round(difference, 2) != 1.00:
        print(
            f'Combination {map_dimensions_x_map_value}, Difference: {difference:.3f}, Expected: {map_value_dict["expected"]}, '
            f'Sum of array: {map_value_dict["observed_from_array"]:.3f}, '
            f'Maptotal: {map_value_dict["observed_from_maptotal"]:.3f}'
        )

Some examples of outputs:

For a map of size 4000x4000 with a constant value of 0.01, the result of maptotal is only 96.3% of what is expected:
Combination 4000x0.01, Difference 0.963, Expected: 160000.0, Sum of array: 160002.875, Maptotal: 154039.000

For a map of size 5000x5000 with a constant value of 1, the result of maptotal is only 67.1% of what is expected:
Combination 5000x1, Difference 0.671, Expected: 25000000, Sum of array: 25000000.000, Maptotal: 16777216.000

Is this indeed a bug? Or is there a hard limit for map sizes that I crossed? Or did I simply misuse pcraster commands? :)

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

Relationships

None yet

Development

No branches or pull requests

Issue actions