3

We need to perform raster calculations using more than two bands in PostGIS. At the moment, ST_MapAlgebra has a two-band form but doesn't support more than that.

I'm open to any solution others would like to offer, but am currently favoring writing a custom plpython function that might be called as follows

SELECT
    ST_FromGDALRaster(
        some_custom_function(
            ST_AsGDALRaster(
                some_raster_column,
                'GTiff'
            )
        )
    ) new_raster_column
FROM some_table

where the custom function accepts and returns a bytea PostgreSQL datatype. Inside the function I am able to import gdal and use following

vsipath = '/vsimem/buffer'
gdal.FileFromMemBuffer(vsipath, gdal_rast)
rast = gdal.Open(vsipath)

gdal_rast above is the parameter being passed into the custom function resulting from the ST_AsGDALRaster function in the select clause. From there, I can turn the bands into numpy arrays and process as desired.

My problem comes in returning the processed raster back from the custom function to the database. I have tried the MEM driver to create an in-memory raster, but when I return that raster to the ST_FromGDALRaster function from the custom function I receive this error:

ERROR: RASTER_fromGDALRaster: Could not open bytea with GDAL.
Check that the bytea is of a GDAL supported format

Can the plpython function approach be salvaged? Or is there a better solution that avoids plpython altogether?

5
  • If you functions are linear, can you not call ST_MapAlgebra twice to combine 3 rasters. I realize that this does not directly answer the question. Jan 14, 2016 at 9:38
  • That is looking like the only solution that will work, I will try that approach.
    – dadbot
    Jan 14, 2016 at 19:41
  • Our data scientist just let me know that the functions are unfortunately non-linear so nesting ST_MapAlgebra functions won't work. Please let me know if there are other ideas out there.
    – dadbot
    Jan 14, 2016 at 20:46
  • 1
    Are you sure postgis only supports 2 bands? The documentation seems to indicate otherwise. postgis.net/docs/RT_ST_MapAlgebra.html
    – tilt
    Jan 15, 2016 at 20:27
  • I was looking at the expression version of ST_MapAlgebra, looks like the rastbandarg of the callback version will meet our needs. Missed that, thanks for calling it to my attention!
    – dadbot
    Jan 18, 2016 at 17:20

1 Answer 1

1

tilt pointed out in the comments above that the rastbandarg version of ST_MapAlgebra's callback form accepts more than two raster bands. This approach solves the problem.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.