# Using custom python function with raster calculator

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.

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.

Update:

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_)
``````

The `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.

• Can you post the code you have now? – GBG Apr 24 '18 at 15:38
• How big are your rasters? nrows, ncols, datatype? If they're not too huge, you can just use gdal.Open('path_to_raster').ReadAsArray() in Python to get each as a numpy array, then apply your functions. The only tricky thing is writing this back to disk, and it's not that tricky. – Jon Apr 27 '18 at 17:30
• Data sets vary in size, but the largest one is roughly 12000 by 50000. I tried ReadAsArray() and it through a memory error. Average size is around 3000 by 3000. – Mike Bannister Apr 27 '18 at 17:45