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?