Handling nan nodata values in Cython modules

[davemfish]

This came up originally in #1705 , but is common to several modules.

It's common for our raster processing routines to mask out nodata pixels before applying some calculation. When this happens in the context of raster_calculator, we typically create a mask using numpy.isclose(array, nodata_value, equal_nan=True). But in the context of a Cython module, we typically define our own is_close function, that looks something like this:

cdef int is_close(double x, double y):
    return abs(x-y) <= (1e-8+1e-05*abs(y))

This function does not correctly handle nan values. And nan is a valid GDAL raster value, and nodata value. is_close could be modified to handle nan values like this:

cdef int is_close(double x, double y):
    if isnan(x) and isnan(y):
        return 1
    return abs(x-y) <= (1e-8+1e-05*abs(y))

That's a solution we implemented in pygeoprocessing.routing.routing.pyx

But there's reason to be concerned about the performance of adding this extra check to a function that is called once for every pixel in a raster stack. I'm not sure if we evaluated that concern when implementing the pygeoprocessing function.

We could consider an alternative approach: Pre-processing user input data to ensure that it does not include numpy.nan values, and to reclassify them to something if it does. I don't know how that would compare performance-wise.