I'd like to count the number of raster pixels that either intersect a boundary or are contained inside the boundary and think ST_Intersects might be a good tool for this.

When using ST_Intersects on an ingested shapefile and a raster does it return true for raster pixels residing completely inside the polygon boundary or only if they overlap the border?

This is the current query I'm using but am not sure if it's picking up the pixels inside the boundary.

SELECT COUNT(*) AS intersections FROM ${rastTable} r, boundaries p WHERE p.boundaryname = '${someBoundary}' AND ST_Intersects(r.rast, p.geom);

Maybe I need to union this with an ST_Contains or is it ok as is? On further inspection I'm not even sure this is counting raster pixels. It could be counting raster tiles which isn't what I'm looking for.


After further investigation I'm pretty sure the above query was counting raster tile intersections and not actual pixels. Here is my revised query which counts tile intersections, pixels (excluding nodatas), and pixels (including nodatas).

SELECT count(*) AS tile_intersections, SUM(ST_Count(r.rast)) as pixels_withdata, SUM(ST_Count(r.rast, false)) as total_pixels
FROM ${rastTable} r, boundaries p
WHERE p.boundaryName = ${someBoundary} AND ST_Intersects(r.rast, p.geom);

So now I just need to verify that this is including pixels contained in the boundary.

EDIT 2 (Some progress)

Messed around with this some more, the above query is also incorrect as it counts all pixels in the tiles intersecting the boundary even if they fall outside the boundary. My next attempt to solve this makes use of ST_Clip. The query does some additional work like output the result as a tiff so that I can see/debug what's happening along with using ST_SummaryStats to get other metrics I was interested in.

SELECT p.gid AS gid, 
ST_AsTIFF(ST_Clip(ST_Union(r.rast), p.geom, true)) AS tiffy,
ST_Count(ST_Clip(ST_Union(r.rast), p.geom, true)) AS pixel_count_excluding_nodata, 
ST_Count(ST_Clip(ST_Union(r.rast), p.geom, true), false) AS pixel_count_including_nodata, 
(ST_SummaryStats(ST_Clip(ST_Union(r.rast), p.geom, true), true)).*
FROM ${rastTable} r, fws_boundaries p
WHERE p.orgname = '${boundary}' 
AND ST_Intersects(r.rast, p.geom) 
AND r.rast_date = '${dateString}' 
GROUP BY p.gid;

The above query is not optimized (I'm sure I shouldn't be clipping 4 times :) and groups by gid because some of the boundary names I use map to multiple polygons.

  • Use OpenJump to display geometries and vectorized raster tiles (ST_DumpAsPolygons(rast)).*. It's very easy. – Pierre Racine Aug 22 '17 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.