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

[JP Glutting]

Thanks Pierre,

That is very helpful.

Unfortunately, the change to ST_Intersects returns an error:

ERROR:  function st_intersects(raster, geometry, boolean) does not exist

LINE 3:  where ST_Intersects(R.rast, A.the_geom, FALSE)

               ^

HINT:  No function matches the given name and argument types. You might need
to add explicit type casts.

There are only small portions of the rasters that have values of 1. However,
the ST_Intersects statement in the where clause seems to filter out all the
empty geometries, so that is what I wanted. I will run a few tests to see if
it is working correctly.

Cheers,
JP

On Thu, Apr 28, 2011 at 3:32 PM, Pierre Racine
<Pierre.Racine at sbf.ulaval.ca>wrote:

> 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
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20110428/bf46c8b7/attachment.html>