[postgis-users] ST_Intersection(rast, the_geom) and ST_Intersects(rast, the_geom)

[Pierre Racine]

Hi JP,

It all depends on the size of your raster coverage. If it is big, it should be tiled (loaded using the raster2pgsal.py -k option and tiled as small as possible 10x10 or 25x25) and in this case the ST_Intersects is useful. If it's not tiled then the ST_Intersects is useless unless you want to eliminate points outside the footprint of your raster coverage.

Another factor is the proportion of nodata value. If there are many, ST_Intersects might be slow because it takes them into account. In this case you could just ignore them by adding FALSE to ST_Intersects ( ST_Intersects(R.rast, V.the_geom, FALSE)), and then parse out the 0s with a WHERE clause:

SELECT V.id, (ptest).geom
FROM (SELECT V.id, ST_Intersection(R.rast, V.the_geom) AS ptest
               FROM raster R, vector V
               WHERE ST_Intersects(R.rast, V.the_geom, FALSE)) foo
WHERE (ptest).val != 0

We could say that ST_Intersects is optimized for the case where the proportion of nodata values is low, which is generally the case.

You should not have to do the same query many time on different rasters. Just load all your rasters in the same table, tiled or not, and use ST_Intersects. Just make sure they are all in the same SRID.

Hope this help.

Pierre

De : postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-bounces at postgis.refractions.net] De la part de JP Glutting
Envoyé : 28 avril 2011 07:39
À : postgis-users at postgis.refractions.net
Objet : [postgis-users] ST_Intersection(rast, the_geom) and ST_Intersects(rast, the_geom)

Hi,

I am using postgis-2.0SVN (revision 7065, I think - I downloaded it on Monday).

I am doing an analysis in which I need to test 50,000 points in a vector table against a binary raster (with 0 set as NoData and 1 being the area of interest). From what I can find online (Amazon has not sent my copy of Postgis in Action yet), the following is the canonical way of doing this:

select V.id, ST_Intersection(R.rast, V.the_geom) as pts
from raster R, vector V
where ST_Intersects(R.rast, V.the_geom)

All I really need is:

select V.id
from raster R, vector V
where ST_Intersects(R.rast, V.the_geom)

and I was hoping that it would be a faster way to do this.

However, when I run a test on my data, however, with "LIMIT 10", the equivalent code takes 24,926 ms.

Doing something more complicated, like this (because ST_Intersection seems to return empty geometries for the NoData areas):

select V.id, ST_IsEmpty((ST_Intersection(R.rast, V.the_geom)).geom) as ptest
from raster R, vector V
LIMIT 10;
-- then I filter these results for ptest = FALSE

takes 589 ms.

I tried scaling this up to LIMIT 100 (359,304  and 6,444 ms, respectively) and 1000 (1,980,397 and 61,628 ms).

This is a huge difference. Is this the best way to use ST_Intersection/ST_Intersects? Am I doing something wrong here? Would I be better off converting this to a vector? Any tips on optimizing this query would be appreciated (I have to run the same query against several different rasters, plus I want to understand how to do this right).

JP
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20110428/b7f133bc/attachment.html>