Abstraction for iterating over (x, y) coords by raster blocks

[emlys]

A very common pattern in our cython raster algorithms looks like:

for offset_dict in pygeoprocessing.iterblocks(raster_path_band, offset_only=True):
    if ctime(NULL)-last_log_time > _LOGGING_PERIOD:
        LOGGER.info(f'{n_processed/n_pixels*100:.1f}% complete')
    last_log_time = ctime(NULL)
    xoff = offset_dict['xoff']
    yoff = offset_dict['yoff']
    win_xsize = offset_dict['win_xsize']
    win_ysize = offset_dict['win_ysize']
    n_processed += win_xsize * win_ysize

    for xi in range(win_xsize):
        for yi in range(win_ysize):
            x = xoff + xi
            y = yoff + yi
           
            ... do something with (x, y)

When using iterblocks to process one block at a time, we need to convert the offset dicts into actual coordinates, and iterate over those one pixel at a time. Implementing this in a reusable function would let us deduplicate a lot of similar code. And it would save us 3 levels of indentation in the subsequent algorithm!

I count at least 5 cases where this could be used just in pygeoprocessing.routing, and I know it could be used some places in invest too.

There is also a similar pattern where we create a buffered array for each block that extends some number of pixels past the edges of the block. That could be abstracted out similarly.

For compatibility with both cython and C++, this would be best implemented as a C++ extension.

[phargogh]

I think this is a great idea! And providing both python and C++ interfaces to this sounds useful as well. To save future us from possibly getting caught with super long iteration times, maybe we could raise a warning when the python interface is used on a raster that's sufficiently large to cause runtime problems?