PostGIS raster MapAlgebra callback optimization

Question

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 
  immutable
as $$
  declare
    roaddist float;
    protected_area float;
    citydist float;
    return_val float;
  begin

    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;
  end;
$$;

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.

Answer 1

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)

Answer 2

Or

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.