1

I'm trying to clip a raster with a shapefile in postgis but want the clipped raster to include pixels that overlap the border of the shapefile which ST_Clip doesn't seem to be including.

Here is a screenshot of the shapefile loaded in qgis:

enter image description here

This query returns the raster clipped using the boundary of the geometry (I'm not sure why it doesn't include the top section).

SELECT ST_AsTIFF(ST_Union(ST_Clip(r.rast, ST_Boundary(foo.area), true))) AS the_tiff
FROM (SELECT p.geom AS area FROM areas p WHERE p.orgname = '${name}') as foo
INNER JOIN ${rastTable} r ON ST_Intersects(r.rast, foo.area);

enter image description here

This query returns the raster clipped using the geometry (without wrapping it in ST_Boundary, notice that the boundary pixels are not included)

SELECT ST_AsTIFF(ST_Union(ST_Clip(r.rast, foo.area, true))) AS the_tiff
FROM (SELECT p.geom AS area FROM areas p WHERE p.orgname = '${name}') as foo
INNER JOIN ${rastTable} r ON ST_Intersects(r.rast, foo.area);

enter image description here

My goal is to have the union of these two rasters (both the interior intersections and the border intersections as I want to run statistics on them). I don't understand why ST_Clip doesn't include pixels in the raster intersecting with the border of the shapefile, is it supposed to? Do I really need to be taking the union or how should I write this query?

Basically, I want the rule to be if any part of the shapefile touches a raster pixel, then include the pixel.

Small update

I can get closer to the desired outcome using ST_buffer(foo.area, .01) but this seems like the wrong way to go about it (and still doesn't fix the upper curved section):

enter image description here

  • Is there a nodata value defined in the raster that would cover the missing part? – Pierre Racine Aug 24 '17 at 14:51
  • @Pierre nope, I've checked the backing raster tiles and while the nodata value is -9999 for the raster, none of the pixels in or around the clipped area have that value. – apricity Aug 24 '17 at 18:57
0

I also found that st_clip is only returning values for the pixels that are entirely covered by the geometry. Seems there should be a way to get st_clip to return full coverage of the polygon but I couldn't find it. To get all the values for intersecting pixels, I didn't use st_clip but st_dumpaspolygons and processed as a vector.

WITH a AS (
-- get rasters that intersect with geometry and dump little raster polygons 
SELECT p.geom AS area 
    , (st_dumpaspolygons(r.rast)).*
FROM ${rastTable} r 
JOIN areas p ON st_intersects(p.geom,r.rast)
) 
-- filter to intersecting polygons 
SELECT geom, val AS 
FROM a 
WHERE st_intersects(area, geom)

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.