If I limit to only 10 polygons it uses 221285 ms.<br><br><div class="gmail_quote">2011/2/28 Andreas Forĝ Tollefsen <span dir="ltr"><<a href="mailto:andreasft@gmail.com">andreasft@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Had the query running for over 48 hours but now it had crashed. Any idea on how to optimize for speed and stability?<div><br></div><div><div>DROP TABLE IF EXISTS globshortrel;</div><div><br></div><div>SELECT </div><div>gid, </div>
<div>((foo.geomval).val),</div><div>CAST (SUM(ST_Area((foo.geomval).geom)) AS decimal(7,5)) as area, </div><div>CAST(SUM(ST_Area((foo.geomval).geom))/0.25*100 AS decimal(6,4)) as percentarea, </div><div>CAST (SUM(ST_Area((foo.geomval).geom))/0.0000077160617284 AS int) AS npixels</div>
<div>INTO globshortrel</div><div class="im"><div>FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo</div>
</div><div class="im"><div>GROUP BY gid, (foo.geomval).val</div>
</div><div>ORDER BY gid;</div><div><br></div><div><br></div><div>-- SELECT * FROM globshortrel;</div><div><br></div><div>-- SELECT DISTINCT ON(gid) gid, val, percentarea</div><div>-- FROM globshortrel</div><div>-- GROUP BY gid, val, percentarea</div>
<div>-- ORDER BY gid, percentarea DESC</div><div>ERROR: out of memory</div><div>DETAIL: Failed on request of size 1753647.</div><div>CONTEXT: SQL function "st_dumpaspolygons" statement 1</div><div>PL/pgSQL function "st_intersection" line 9 at RETURN QUERY</div>
<div>SQL function "st_intersection" statement 1</div><div><br></div><div>********** Error **********</div><div><br></div><div>ERROR: out of memory</div><div>SQL state: 53200</div><div>Detail: Failed on request of size 1753647.</div>
<div>Context: SQL function "st_dumpaspolygons" statement 1</div><div>PL/pgSQL function "st_intersection" line 9 at RETURN QUERY</div><div>SQL function "st_intersection" statement 1</div><div>
<div></div><div class="h5"><br><div class="gmail_quote">
2011/2/25 Andreas Forĝ Tollefsen <span dir="ltr"><<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I wrote this query. It seems to give me the area, percentage of area and the pixelcount within a 0.5x0.5 cell.<div><br><div><div>SELECT gid, SUM(ST_Area((foo.geomval).geom)) as area, SUM(ST_Area((foo.geomval).geom))/0.25*100 as percentarea, SUM(ST_Area((foo.geomval).geom))/7.716049382716049382716049382716e-6 AS Npixels</div>
<div>FROM (SELECT onecell.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(onecell.rast, priogrid_land.cell) AS geomval FROM onecell, priogrid_land) AS foo</div><div>WHERE gid = 139358</div><div><div>
GROUP BY gid, ((foo.geomval).val)</div>
</div><div>ORDER BY area, ((foo.geomval).val)</div><div><br></div><div>Result: </div><div>gid; area; percentarea; npixels</div><div><div>139358;0.000123456790106502;0.0493827160426008;15.9999999978027</div><div>139358;0.000563271604818283;0.225308641927313;72.9999999844495</div>
<div>139358;0.0966820987653989;38.6728395061596;12529.9999999957</div><div>139358;0.152631172839676;61.0524691358705;19781.0000000221</div></div><div><br></div><br><div class="gmail_quote">2011/2/25 Andreas Forĝ Tollefsen <span dir="ltr"><<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>></span><div>
<div></div><div><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Pierre. I did tile it when i processed it.<div>I used this gdal2wktraster syntax: <span lang="EN-US" style="font-size:9.0pt;line-height:150%;font-family:"Times New Roman","serif";color:black">python gdal2wktraster.py -r c:\prio_grid\source\globcover\globshort.tif
-t globshort -s 4326 -k 720x360 -I -o globshort.sql</span></div><div><font face="'Times New Roman', serif"><span style="font-size:12px;line-height:18px"><br></span></font></div>
<div><span lang="EN-US" style="color:black"></span><font face="'Times New Roman', serif"><span style="font-size:12px;line-height:18px">The size of my raster is X: 129600 Y: 55800<br>
</span></font><br><div class="gmail_quote">2011/2/24 Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca" target="_blank">Pierre.Racine@sbf.ulaval.ca</a>></span><div><div></div><div>
<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div lang="EN-US" link="blue" vlink="purple"><div><p class="MsoNormal"><span style="font-size:11.0pt;color:#1F497D">You should just divide the area of your polygon by the area of one of your pixel.</span></p><p class="MsoNormal">
<span style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span style="font-size:11.0pt;color:#1F497D">I want to investigate this example a little bit further. What is the size of your raster in pixels? I understand that you did not tile it. Did you?</span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:#1F497D"> </span></p><p class="MsoNormal"><span style="font-size:11.0pt;color:#1F497D">Pierre</span></p><p class="MsoNormal"><span style="font-size:11.0pt;color:#1F497D"> </span></p>
<p class="MsoNormal"><span 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">From:</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>On Behalf Of </b>Andreas Forĝ Tollefsen<br>
<b>Sent:</b> 24 février 2011 09:53<br><b>To:</b> Paragon Corporation<br><b>Cc:</b> PostGIS Users Discussion</span></p><div><div></div><div><br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon</div></div>
<p></p></div></div><div><div></div><div><p class="MsoNormal"> </p><p class="MsoNormal">Great.</p><div><p class="MsoNormal">Thank you so much. I should have noticed that these were unioned numbers. </p></div><div>
<p class="MsoNormal">Btw. are there any way of getting the count of pixels within a polygon without actually aggregating them or union them?</p></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal" style="margin-bottom:12.0pt">
Andreas</p><div><p class="MsoNormal">2011/2/24 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>></p><div><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">Andreas,</span></p>
<p class="MsoNormal"><span style="font-size:10.0pt;color:blue">Sorry should have recognized what you're doing. The intersection returns a polygon which is a union of the clipped raster pixel squares. So you need to use Sum of area instead and then divide by the area of a pixel to get the equivalent of your count.</span></p>
<p class="MsoNormal"> </p><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">So </span></p><p class="MsoNormal"> </p><div><p class="MsoNormal">SELECT gid, SUM(ST_Area((foo.geomval).geom))/ [put your pixel area size here] as ct</p>
</div><div><div><p class="MsoNormal">FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo</p></div><div><p class="MsoNormal">
WHERE gid >= 139358 AND gid <= 139365</p></div><div><p class="MsoNormal">GROUP BY gid</p></div><div><p class="MsoNormal">ORDER BY gid</p></div></div><p class="MsoNormal"> </p><div class="MsoNormal" align="center" style="text-align:center">
<hr size="2" width="100%" align="center"></div><p class="MsoNormal"><b><span style="font-size:10.0pt">From:</span></b><span style="font-size:10.0pt"> Andreas Forĝ Tollefsen [mailto:<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>] <br>
<b>Sent:</b> Thursday, February 24, 2011 8:33 AM<br><b>To:</b> PostGIS Users Discussion<br><b>Cc:</b> Paragon Corporation</span></p><div><div><p class="MsoNormal"><span style="font-size:10.0pt"><br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon</span></p>
</div></div><p class="MsoNormal"> </p><div><div><p class="MsoNormal">I am a bit unsure whether my results are actually correct. According to a total count using the below query, I get very different results between the cells. </p>
<div><p class="MsoNormal">Since the raster does actually cover the whole vector cell, i would assume that the count should be similar in all cells. Meaning, the pixel count should be the same.</p></div><div><p class="MsoNormal">
What i get is different, and it seems that the query is not providing me with the number of pixels within the grid cell.</p></div><div><p class="MsoNormal">Any idea why this is so different?</p></div><div><p class="MsoNormal">
</p></div><div><div><p class="MsoNormal">SELECT gid, count((foo.geomval).val) as ct</p></div><div><p class="MsoNormal">FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo</p>
</div><div><p class="MsoNormal">WHERE gid >= 139358 AND gid <= 139365</p></div><div><p class="MsoNormal">GROUP BY gid</p></div><div><p class="MsoNormal">ORDER BY gid</p></div><div><p class="MsoNormal"> </p></div><div>
<p class="MsoNormal">Result:</p></div><div><div><p class="MsoNormal">139358;632</p></div><div><p class="MsoNormal">139359;1030</p></div><div><p class="MsoNormal">139360;912</p></div><div><p class="MsoNormal">139361;731</p>
</div><div><p class="MsoNormal">139362;760</p></div><div><p class="MsoNormal">139363;1230</p></div><div><p class="MsoNormal">139364;1314</p></div><div><p class="MsoNormal">139365;1014</p></div></div><div><p class="MsoNormal">
</p></div><div><p class="MsoNormal">The attached image shows the raster pixels within one cell.</p></div><div><p class="MsoNormal"> </p></div><p class="MsoNormal"> </p><div><p class="MsoNormal">2011/2/24 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>></p>
<p class="MsoNormal">Thanks! </p><div><p class="MsoNormal">That solved it.</p></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">This will probably take a lot of time. I have 259200 polygons measuring 0.5 x 0.5 decimal degrees while the raster dataset is of global cover and has a pixelsize of 0.00277777777777778x0.00277777777777778. </p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal"><span style="color:#888888">Andreas</span></p></div><div><div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal"> </p><div><p class="MsoNormal">
2011/2/23 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>></p><div><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">Andrea,</span></p><p class="MsoNormal"> </p><p class="MsoNormal">
<span style="font-size:10.0pt;color:blue">Try </span></p><p class="MsoNormal"> </p><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">SELECT DISTINCT ON(gid) gid, </span><span style="color:blue">(foo.geomval).val, COUNT((foo.geomval).val) AS ct</span></p>
<div><div><p class="MsoNormal"><span style="color:blue">FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo</span></p>
</div><div><p class="MsoNormal"><span style="color:blue">WHERE gid > 151000 AND gid < 151010</span></p></div><div><p class="MsoNormal"><span style="color:blue">GROUP BY gid, (foo.geomval).val</span></p></div></div>
<div>
<p class="MsoNormal"><span style="color:blue">ORDER BY gid, ct DESC</span></p></div><div><p class="MsoNormal"> </p></div><div class="MsoNormal" align="center" style="text-align:center"><hr size="2" width="100%" align="center">
</div><div><p class="MsoNormal"><b><span style="font-size:10.0pt">From:</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>On Behalf Of </b>Andreas Forĝ Tollefsen</span></p>
</div><p class="MsoNormal" style="margin-bottom:12.0pt"><b><span style="font-size:10.0pt">Sent:</span></b><span style="font-size:10.0pt"> Wednesday, February 23, 2011 4:05 AM<br><b>To:</b> PostGIS Users Discussion<br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon</span></p>
<div><div><p class="MsoNormal">Hi. Thanks Regina and Leo, </p><div><p class="MsoNormal">I have been testing the raster and geom intersection a bit. I guess what i need is to use the ST_Intersection together with a max(count) function. </p>
<div><p class="MsoNormal">So my result will be the rastervalue with the highest count within each of the grid cells.</p></div><div><p class="MsoNormal">However, as far as i know, there is now Max(COUNT) function in postgresql.</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">Any idea how i can modify the below query to only return the rastervalue within the grid cell occuring most frequently?</p></div><div><p class="MsoNormal">
Consequently i want only one row for each gid, and the maximum occuring rastervalue.</p></div><div><p class="MsoNormal"> </p></div><div><div><p class="MsoNormal">SELECT gid, (foo.geomval).val, COUNT((foo.geomval).val) AS ct</p>
</div><div><p class="MsoNormal">FROM (SELECT globshort.rid, priogrid_land.cell, priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land) AS foo</p></div><div><p class="MsoNormal">
WHERE gid > 151000 AND gid < 151010</p></div><div><p class="MsoNormal">GROUP BY gid, (foo.geomval).val;</p></div></div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">gid; val; ct</p></div><div><div>
<p class="MsoNormal">
151001;14;381</p></div><div><p class="MsoNormal">151001;150;9</p></div><div><p class="MsoNormal">151001;50;7</p></div><div><p class="MsoNormal">151001;140;91</p></div><div><p class="MsoNormal">151001;40;1</p></div><div><p class="MsoNormal">
151001;70;2</p></div><div><p class="MsoNormal">151001;130;4</p></div><div><p class="MsoNormal">151001;200;48</p></div><div><p class="MsoNormal">151001;100;3</p></div><div><p class="MsoNormal">151001;;0</p></div><div><p class="MsoNormal">
151001;190;1</p></div><div><p class="MsoNormal">151001;20;203</p></div><div><p class="MsoNormal">151001;11;111</p></div><div><p class="MsoNormal">151001;210;16</p></div><div><p class="MsoNormal">151001;30;105</p></div></div>
<div><p class="MsoNormal"> </p></div><div><p class="MsoNormal"> </p><div><p class="MsoNormal">2011/2/23 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>></p><div><p class="MsoNormal">
<span style="font-size:10.0pt;color:blue">Have you looked at ST_Intersection. I'm not sure how large your grids are so might still be a bit too slow. </span></p><p class="MsoNormal"> </p><p class="MsoNormal"> </p><p class="MsoNormal">
<span style="color:blue"><a href="http://www.postgis.org/documentation/manual-svn/RT_ST_Intersection.html" target="_blank">http://www.postgis.org/documentation/manual-svn/RT_ST_Intersection.html</a></span></p><p class="MsoNormal">
</p><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">Below is a link to our slides from our North Carolina GIS meeting that may answer some of your questions (shows some Raster examples) as well as the 3D ones people have asked.</span></p>
<p class="MsoNormal"> </p><p class="MsoNormal"><span style="color:blue"><a href="http://www.postgis.us/presentations" target="_blank">http://www.postgis.us/presentations</a></span></p><p class="MsoNormal"> </p><div><p class="MsoNormal">
<span style="font-size:10.0pt;color:blue">Hope that helps,</span></p></div><div><p class="MsoNormal"><span style="font-size:10.0pt;color:blue">Regina and Leo</span></p></div><div class="MsoNormal" align="center" style="text-align:center">
<hr size="2" width="100%" align="center"></div><p class="MsoNormal" style="margin-bottom:12.0pt"><b><span style="font-size:10.0pt">From:</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>On Behalf Of </b>Andreas Forĝ Tollefsen<br>
<b>Sent:</b> Tuesday, February 22, 2011 4:28 AM<br><b>To:</b> PostGIS Users Discussion<br><b>Subject:</b> [postgis-users] ST_Value from Polygon</span></p><div><div><p class="MsoNormal">Hi all, </p><div><p class="MsoNormal">
</p></div><div><p class="MsoNormal">I am working with a large raster dataset that i want to aggregate into vector grids.</p></div><div><p class="MsoNormal">The raster dataset is a landcover dataset, and i want to find which of the raster values are the most dominant within each of the vector grid cells.</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">I have been looking at the ST_Value function, but this is not usable together with the cell polygon.</p></div><div><p class="MsoNormal"> </p></div><div>
<p class="MsoNormal">I have written a script that gives me the raster value of the centroid of each cell, but i want to find which raster class is the largest.</p></div><div><p class="MsoNormal">Hence i need to calculate the area of each raster class within each cell and select the largest class.</p>
</div><div><p class="MsoNormal"> </p></div><div><p class="MsoNormal">Any idea? So far i have only come this far:</p></div><div><p class="MsoNormal"> </p></div><div><div><p class="MsoNormal">DROP TABLE IF EXISTS globshortpoly;</p>
</div><div><p class="MsoNormal">SELECT priogrid_land.cell, ST_Value(rast, ST_Centroid(cell))</p></div><div><p class="MsoNormal">INTO globshortpoly</p></div><div><p class="MsoNormal">FROM priogrid_land, globshort</p></div>
<div><p class="MsoNormal">WHERE rast && priogrid_land.cell </p></div><div><p class="MsoNormal">LIMIT 1000</p></div></div></div></div></div><p class="MsoNormal" style="margin-bottom:12.0pt"><br>_______________________________________________<br>
postgis-users mailing list<br><a href="mailto:postgis-users@postgis.refractions.net" target="_blank">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></p>
</div><p class="MsoNormal"> </p></div></div></div></div></div><p class="MsoNormal" style="margin-bottom:12.0pt"><br>_______________________________________________<br>postgis-users mailing list<br><a href="mailto:postgis-users@postgis.refractions.net" target="_blank">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></p></div><p class="MsoNormal"> </p></div></div></div></div><p class="MsoNormal">
</p></div></div></div></div></div><p class="MsoNormal"> </p></div></div></div></div></div></div><br>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@postgis.refractions.net" target="_blank">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></div></div><br></div>
</blockquote></div></div></div><br></div></div>
</blockquote></div><br></div></div></div>
</blockquote></div><br>