I have a function in a python class that takes two values and returns a third value. I would like to apply the function to two overlapping rasters on a per pixel level and come up with a third raster. It doesn't appear that either the qgis or arcgis raster calculators can accept a custom function. I'm sure this can be done with GDAL, but don't have any experience with it. GDAL is certainly an option if an easier system doesn't exist.
I have access to QGIS, arcmap 10.6, and arcgis pro.
For a little context I'm using an NDVI raster and a vegetation height raster to calculate hydraulic surface roughness based on a table I've already created. Worse case situation will be to write a quick script to convert the table to a raster calculator equation. As the table has 49 "bins" the equation will be quite long and it seems like there should be a more elegant way to do this.
Simplified code for returning n value based on NDVI and vegetation height, a lot of error checking has been removed for clarity:
class RoughnessTable(object): """ Extracts predicted mean and std Manning's n values for a XS based on lidar and NDVI v0.01b """ def __init__(self, table): """ table is the lidar/ndvi/n-value dataframe""" self.table = table def mean(self, l, n): """ Return mean value from roughness table for lidar height 'l' and NDVI 'n' :param l: mean lidar height of xs segment :param n: mean ndvi of xs segment :returns: mean predicted surface roughness (float) """ rt_row = self.table.query('l_bin_min <= @l < l_bin_max & n_bin_min <= @n < n_bin_max') return float(rt_row.mean_)
table passed to
RoughnessTable.__init__() is a pandas dataframe that represents the NDVI/vegetation height bins. A sample is shown below.
I've had success using the raster_calculator from the pygeoprocessing library, which handles opening the raster, and splitting it into 256 x 256 pixel ndarray chunks. I'm still having to iterate over the chunk from raster_calculator, and it seems rather hacky, but it works. It is rather slow. I haven't profiled it yet, but I'm guessing the
self.table.query(...) line is causing a large portion of the problem.