2

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?

  • 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. – John Powell Jan 14 '16 at 9:38
  • That is looking like the only solution that will work, I will try that approach. – dadbot Jan 14 '16 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 '16 at 20:46
  • Are you sure postgis only supports 2 bands? The documentation seems to indicate otherwise. postgis.net/docs/RT_ST_MapAlgebra.html – tilt Jan 15 '16 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 '16 at 17:20
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, privacy policy and cookie policy

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