0

I need to overlay a geometry (Shapefile of some city in UK) over a raster (grids of entire UK Map) and extract cells/pixels from raster that intersect with points of geometry.

Current Implementation: I am able to extract the points using below functions:

ST_PixelAsPoints(ST_MapAlgebra(rast,1,ST_AsRaster(geom,rast,'8BUI',1,0,True),1,'[rast.val]')))

Problem with above, is that it extracts all the points, including those points which just touch the raster during overlay.

Expected Implementation:

I need to extract only those points/pixels which overlap above a certain % on the raster. i.e. I need to extract points based on % of intersection.

In essence , I need to replicate the functionality provided by extract function of R . Also more importantly , I need to replicate the weights option provided by extract function. https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/extract

  • 1
    Welcome to GIS SE. As a new user, please take the Tour. We use a "Focused question / Best answer" model here, which relies on questions that not only state what is needed, but also what has been done toward accomplishing it, and what problem had been encountered. – Vince Feb 5 '18 at 12:49
1

Try something like this (untested):

WITH pixelsvalues AS (
  SELECT ST_PixelAsPolygons(rast) gv
  FROM rasttable, polytable
  WHERE ST_Intersects(rast, geom)
)
SELECT (gv).val,
       ST_Intersection((gv).geom, p.geom) geom
FROM pixelsvalues, polytable p
WHERE ST_Intersects((gv).geom, p.geom) AND ST_Area(ST_Intersection((gv).geom, p.geom))/ST_Area((gv).geom) > 0.05
  • Thanks for answer. I tried with the above code. I need more help with the above code to complete my requirements. – Sachin Feb 8 '18 at 16:42
  • Thanks for answer. I need more help with the above code to complete my requirements. I need to extract points from raster . I am bit lost on how to accomplish it. I tried with ST_PixelAsPoints , but points don't have an area hence it fails . I extracted x and y values from gv and tried to extract points from raster using ST_pixelAsPoint(raster,x,y) -> but this interestingly fails to provide me all the points, especially the points in geometry which are just touching the raster. I get 16 rows on my test data with above query , but they are polygons , I need to get 16 points instead . – Sachin Feb 8 '18 at 17:18

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.