Thanks Pierre, <div><br></div><div>That is very helpful.</div><div><br></div><div>Unfortunately, the change to ST_Intersects returns an error:</div><div><br></div><blockquote class="webkit-indent-blockquote" style="margin: 0 0 0 40px; border: none; padding: 0px;">
<div><p class="p1">ERROR: function st_intersects(raster, geometry, boolean) does not exist</p></div><div><p class="p1">LINE 3: where ST_Intersects(R.rast, A.the_geom, FALSE)</p></div><div><p class="p1"> ^</p>
</div><div><p class="p1">HINT: No function matches the given name and argument types. You might need to add explicit type casts.</p></div></blockquote><div>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.</div>
<div><br></div><div>Cheers,</div><div>JP<br><div><br><div class="gmail_quote">On Thu, Apr 28, 2011 at 3:32 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div lang="FR" link="blue" vlink="purple"><div><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">Hi JP,</span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">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 <span>ST_Intersects</span> is useful. If it’s not tiled then the <span>ST_Intersects</span> is useless unless you want to eliminate points outside the footprint of your raster coverage.</span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">Another factor is the proportion of <span>nodata</span> value. If there are many, <span>ST_Intersects</span> might be slow because it takes them into account. In this case you could just ignore them by adding FALSE to <span>ST_Intersects</span> ( <span>ST_Intersects</span>(<span>R.rast</span>, <span>V.the_geom</span>, FALSE)), and then parse out the 0s with a WHERE clause:</span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">SELECT V.id, (<span>ptest</span>).<span>geom</span></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">FROM (SELECT V.id, <span>ST_Intersection</span>(<span>R.rast</span>, <span>V.the_geom</span>) AS <span>ptest</span></span></p><p class="MsoNormal">
<span lang="EN-US" style="font-size:11.0pt;color:#1F497D"><span> </span><span> </span>FROM raster R, vector V</span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"><span> </span><span> </span>WHERE <span>ST_Intersects</span>(<span>R.rast</span>, <span>V.the_geom</span>, FALSE)) <span>foo</span></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">WHERE (<span>ptest</span>).<span>val</span> != 0</span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">We could say that <span>ST_Intersects</span> is optimized for the case where the proportion of <span>nodata</span> values is low, which is generally the case.</span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">You should not have to do the same query many time on different <span>rasters</span>. Just load all your <span>rasters</span> in the same table, tiled or not, and use <span>ST_Intersects</span>. Just make sure they are all in the same SRID.</span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">Hope this help.</span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D">Pierre</span></p><p class="MsoNormal"><span lang="EN-US" style="font-size:11.0pt;color:#1F497D"> </span></p><div style="border:none;border-left:solid blue 1.5pt;padding:0cm 0cm 0cm 4.0pt">
<div><div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm"><p class="MsoNormal"><b><span style="font-size:10.0pt">De :</span></b><span style="font-size:10.0pt"> <a href="mailto:postgis-users-bounces@postgis.refractions.net" target="_blank">postgis-users-bounces@postgis.refractions.net</a> [mailto:<a href="mailto:postgis-users-bounces@postgis.refractions.net" target="_blank">postgis-users-bounces@postgis.refractions.net</a>] <b>De la part de</b> JP Glutting<br>
<b>Envoyé :</b> 28 avril 2011 07:39<br><b>À :</b> <a href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a><br><b>Objet :</b> [postgis-users] ST_Intersection(rast, the_geom) and ST_Intersects(rast, the_geom)</span></p>
</div></div><div><div></div><div class="h5"><p class="MsoNormal"> </p><div><p class="MsoNormal">Hi, </p><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">I am using postgis-2.0SVN (revision 7065, I think - I downloaded it on Monday).</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">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:</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">select V.id, ST_Intersection(R.rast, V.the_geom) as pts</p></div><div><p class="MsoNormal">from raster R, vector V</p></div><div><p class="MsoNormal">where ST_Intersects(R.rast, V.the_geom)</p>
</div><div><p class="MsoNormal"> </p></div><div><div><p class="MsoNormal">All I really need is:</p></div><div><p class="MsoNormal"> </p></div></div><blockquote style="margin-left:30.0pt;margin-right:0cm"><div><div><p class="MsoNormal">
<span>select V.id</span></p></div></div><div><div><p class="MsoNormal"><span>from raster R, vector V</span></p></div></div><div><div><div><p class="MsoNormal"><span>where ST_Intersects(R.rast, V.the_geom)</span></p></div>
</div></div></blockquote><div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">and I was hoping that it would be a faster way to do this.</p></div></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">
However, when I run a test on my data, however, with "LIMIT 10", the equivalent code takes 24,926 ms. </p></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">Doing something more complicated, like this (because ST_Intersection seems to return empty geometries for the NoData areas):</p>
</div><div><p class="MsoNormal"> </p></div><blockquote style="margin-left:30.0pt;margin-right:0cm"><div><p class="MsoNormal"><span>select V.id, ST_IsEmpty((ST_Intersection(R.rast, V.the_geom)).geom) as ptest</span></p></div>
<div><p class="MsoNormal"><span>from raster R, vector V</span></p></div><div><p class="MsoNormal"><span>LIMIT 10;</span></p></div></blockquote><blockquote style="margin-left:30.0pt;margin-right:0cm"><div><p class="MsoNormal">
-- then I filter these results for ptest = FALSE</p></div></blockquote><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">takes 589 ms.</p></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">I tried scaling this up to LIMIT 100 (359,304 and 6,444 ms, respectively) and 1000 (1,980,397 and 61,628 ms). </p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">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).</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">JP</p></div></div></div></div></div></div></div><br>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
<a href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
<br></blockquote></div><br></div></div>