PostGIS raster calculate average slope of 3x3 neighborhood around a geometry point

Question

I have two tables

1. corners_final (vector point data)
2. dem_slope (1ft slope DEM, calculated in degrees)

With this simple query I am able to extract the slope value that intersects a given point

SELECT ST_Value(rast,geom)
FROM dem_slope 
JOIN corners_final ON ST_Intersects(geom,rast)
WHERE cornerid=1160060

But I want to find the average slope around that point by 3x3 area.

I cannot figure out how to use ST_Neighborhood... Should I be using ST_MapAlgebra?

I do not see a lot of examples on these functions.

Answer 1

It is bit tricky, especially at raster edges so I came up with this:

select avg(value), cornerid 
from corners_final cf
join lateral (select st_union(st_clip(rast, st_buffer(cf.geom, 10))) raster from dem_slope where rast && st_buffer(cf.geom, 10)) slope on true
join lateral (select unnest(st_neighborhood(raster, cf.geom, 1,1)) value) vals on true
where cornerid=1160060
group by cornerid

It is necessary to retrieve all rasters that are within a 3x3 window around your points, to be able to compute the neighborhood later on. This is done by firs join in the query. Please adjust the buffer amount according to your requirements, it should be a bit bigger than 3x3 pixel diameter. So after first join, each point has its own raster, that contains only pixels around them. Then, in second join, st_neighborhood function is used to retrieve values from the desired neighborhood. Because st_neighborhood returns 2d array, unnest is necessary to retrieve individual pixel values. After that, it is only a matter of siple everage of all values for a specified point (cornerid). Group by is added only for convenience, this way you can specify multiple corner id's, or ewen avoid that condition and still be able to retrieve average values per corner id.

I have tested this query on different datasets (municipality points and DEM raster), and it performs pretty good... about 200ms/100 points.

It is also worth notice, that if you have small rasters, then it may be faster to simply union them without clipping...

Hope that helps.