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

The table passed to RoughnessTable.__init__() is a pandas dataframe that represents the NDVI/vegetation height bins. A sample is shown below.

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.