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.

Vince
  • 20,017
  • 15
  • 45
  • 64
ziggy
  • 4,515
  • 3
  • 37
  • 83

1 Answers1

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.

DavidP
  • 1,353
  • 7
  • 13
  • awesome - this worked and nice use of the double join lateral! – ziggy Nov 07 '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 07 '19 at 15:51
  • is the clip+union of the buffer and raster necessary? – ziggy Nov 07 '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 08 '19 at 07:01
  • But. if you are OK with marginal errors on raster edges, where complete neighborhood cannot be computed, then use your version without union/clip/buffer. – DavidP Nov 08 '19 at 07:08
  • ahhh okay now I understand what you mean about the raster edges – ziggy Nov 08 '19 at 16: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
  • If the query is already going to be using buffers to select the appropriate tiles, I think you might as well at that point just use ST_Clip, and ST_SummarizeStatsAgg, which seems to run faster than this approach- – chrismarx Oct 05 '20 at 19:54