2

I am trying to understand the behaviour of ST_Intersects() in PostGIS.

My goal is to get the mean value of a raster area corresponding to each polygon in a vector layer.

To test my SQL I run initially the script with a vector layer containing only one geometry:

CREATE TABLE mean AS 
SELECT id, (ST_SummaryStats(ST_Clip(rast,1,geom))).mean
FROM vector, raster
WHERE ST_Intersects(geom,rast)

What returns is a two rows table:

1;-146.258943781942
1;-115.587301587302

How it is possible that the same geometry intersects the same raster twice with two different values?

I don't know if it is a coincidence, but the mean value I got doing the same operation with the same two layers with the QGIS Zonal Statistics is -138.788659793814 (the mean between -146.258943781942 and -115.587301587302 is -130... pretty close)

1

Your geometry straddles several tiles. For example try

CREATE TABLE mean AS 
SELECT id, rid, (ST_SummaryStats(ST_Clip(rast,1,geom))).*
FROM vector, raster
WHERE ST_Intersects(geom,rast);

and you will see the rid for each tile in the intersect. You need to summarise the information across the tiles.

You can use ST_SummaryStatsAgg

SELECT gid, (ST_SummaryStatsAgg(ST_Clip(rast, geom, true))).*
FROM raster, vector
WHERE ST_Intersects(rast, geom) /*and gid = 1*/
GROUP BY gid;
-1

Most likely there is one record on the very edge which gets included in the result (as is it true it does intersect).

You probably what to use ST_Within or ST_Contains

http://postgis.net/docs/ST_Within.html

http://postgis.net/docs/ST_Contains.html

  • Why would this happen with ST_SummaryStats() though? Surely that is being calculated on all of the raster data that intersects? – djq Feb 14 '13 at 14:16
  • @Mapperz How is this? My vector layer contains only #1 record, a geometry with unique id "1". The output confirms this since the id "1" is reported in both rows. Do you mean the vector intersects with two records of the raster? – Francesco Feb 14 '13 at 14:51
  • Perhaps it is returning edge id as a vector would have at least two edges – lewis Feb 15 '13 at 19:30

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.