3

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.

3

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.

| improve this answer | |
  • awesome - this worked and nice use of the double join lateral! – ziggy Nov 7 '19 at 14:31
  • this query also got me the same answer-- select cornerid,avg(v) from( select cornerid, unnest(st_neighborhood(rast, geom, 1,1)) v from corners_final join dem_slope on st_intersects(geom,rast) where cornerid=1160060) a group by cornerid – ziggy Nov 7 '19 at 15:51
  • is the clip+union of the buffer and raster necessary? – ziggy Nov 7 '19 at 15:52
  • I think it is! If your point lies on edge of the raster (on first/last row/column), then there are no neighboring pixels for that point on the same raster. So you have to query the raster table for all rows/raster that can be intersected by a 3x3 window around that point. This can yield more than one raster, which have to be unioned back to one raster for further processing. The clip is not strictly necessary, but in my case, it speeds up processing by a factor of 2-3. – DavidP Nov 8 '19 at 7:01
  • 1
    I have added some other nieghborhoods to this query and have found that using a cte to pre-process the buffer on the point geom and removing the join laterals led to a massive speed improvement – ziggy Nov 26 '19 at 15:25

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.