Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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)

share|improve this question

2 Answers 2

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;
share|improve this answer

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

share|improve this answer
    
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

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.