I'm trying to speed up an ST_MapAlgebra callback function that essentially takes a raster with four bands, checks each band for values less than a given number, and where the check is true, return the fourth band.

This is my messy function:

create or replace function 
callback_fn(pixel float[], pos integer[], variadic userargs text[]) 
  returns float 
  language plpgsql 
as $$
    roaddist float;
    protected_area float;
    citydist float;
    return_val float;

    citydist := pixel[1][1][1];
    roaddist := pixel[2][1][1];
    protected_area := pixel[3][1][1];
    return_val := pixel[4][1][1];

    if protected_area = 1 then return null;
    end if;
    if citydist <= 2000 then return return_val;
    elsif roaddist <= 500 then return return_val;
    else return null;
    end if;

So from what I understand, the function above will go through every pixel and either return the "return_val" or it will return NULL. So what I end up with is basically a single band raster with excluded areas removed.

This is exactly what I want, but I was hoping to run this faster. Is there a way to run multiple pixels simultaneously instead of going through each one individually? I'm confused why ST_Reclass can run so quickly to reclassify a whole raster while ST_MapAlgebra seems to take so much longer when (in my mind, at least) the functions are doing similar things - checking the values of a raster and where true, change value, otherwise set null. Is there a way to use ST_Reclass in a callback function to speed up this process?

Additionally, is there a good source for documentation on writing PostGIS raster callback functions? I have been googling a lot over the last few days and have found extremely meager examples.

  • I think the difference is that ST_Reclass is operating on a single raster, whereas MapAlgebra is designed to work on multiple rasters and with offsets from the current pixel, which allows for more complex functions, but which need to be coded up with plpgsql. In your case, you could use ST_Reclass, as you don't seem to be combining bands directly. I agree that the documentation for ST_MapAlgebra is scarce, though the functionality is great. Feb 21 '16 at 8:15

I think to achieve what you want, you'd

1) use ST_Reclass to reclassify your citydist, roaddist .. http://postgis.net/docs/RT_ST_Reclass.html

2) Convert those to geometry with ST_Polygon - http://postgis.net/docs/RT_ST_Polygon.html

3) Then use ST_Clip with your final- http://postgis.net/docs/RT_ST_Clip.html (setting crop to false)



1) ST_Reclass() to get bands of 1 and NULL

2) ST_Intersection(ST_Intersection(ST_Intersection(band1, band2, 'BAND2'), band3, 'BAND2'), band4, 'BAND2')

I guess any call to some standard API function will always be faster than a custom ST_MapAlgebra() function. Also anything which is raster only should be faster than any conversion to vector. Every case can be solve in many ways and all require test and probing.

ST_ReClass() is faster because it has access to low level code and memory that ST_MapAlgebra() can't have access to. ST_MapAlgebra() also has the overhead of having to call a high level function for each pixel.

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.