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

[JP Glutting]

Ah, ok, thanks. I will try that - it looks pretty slick.

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

> My mistake : The FALSE comes right after the raster parameter:
>
>
>
> ST_Intersects(R.rast, FALSE, A.the_geom)
>
>
>
> 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 12:28
> *À :* PostGIS Users Discussion
> *Objet :* Re: [postgis-users] ST_Intersection(rast, the_geom) and
> ST_Intersects(rast, the_geom)
>
>
>
> 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
>
>
>
> _______________________________________________
> 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/5e7b34e6/attachment.html>