Proposal: raster_reduce for raster operations that don't produce a raster

[emlys]

It could be useful to have an equivalent of raster_calculator for operations that don't produce a raster output. If raster_calculator is like map over the blocks of a raster, then this would be like functools.reduce over the blocks of a raster.

Aggregations are probably the most common use case. For the Morgan Stanley project I needed to produce histogram data for a raster. raster_calculator seems like the easiest approach, but it requires writing to a target raster which was totally irrelevant for my use case. I ended up with:

def raster_histogram(raster_path, bins=10):
    """Return histogram data for the given raster divided into bins."""

    histogram = numpy.zeros(bins)

    def add_to_histogram(array):
        nonlocal histogram
        histogram += numpy.histogram(array, bins=bins)[0]
        return array

    pygeoprocessing.raster_calculator(
        [(raster_path, 1)],
        add_to_histogram,
        'foo.tif',
        gdal.GDT_Int32,
        TARGET_NODATA
    )
    return histogram

Given a proposed new function raster_reduce(raster_path, function, initializer) that behaves like functools.reduce, that could be rewritten as

def raster_histogram(raster_path, bins=10):
    """Return histogram data for the given raster divided into bins."""

    def add_to_histogram(array, histogram):
        return histogram + numpy.histogram(array, bins=bins)[0]

    return pygeoprocessing.raster_reduce(
        [(raster_path, 1)],
        add_to_histogram,
        numpy.zeros(bins))

I'm not sure if there are any use cases for this function in invest, but I imagine it would be useful for data analysis.

[emlys]

This would be easy to implement with iterblocks too, like

def raster_reduce(raster_path, function, initializer):
    aggregator = initializer
    for (_, block) in iterblocks(raster_path):
        aggregator = function(aggregator, block)
    return aggregator

The main value of this wrapper function would be so we don't have to think about blocks at all. And I wonder if there's some boilerplate code in raster_calculator that would be useful here too.

[phargogh]

👍 Great idea!

As for other use cases, a few come to mind:

• pygeoprocessing.zonal_statistics seems like it could be rewritten in terms of a reduce function
• There are places in InVEST where we're doing some kind of iterblocks-based accumulation, for example:
  • Crop production (example, another, and another), returning float numbers
  • Scenic Quality, returning a dict of viewpoint coordinates that are over valid DEM pixel values
  • Carbon, for accumulating carbon totals

[emlys]

In our meeting on 12/9 we approved this without a design doc.

[emlys]

This was implemented in #292