<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=iso-8859-1"><meta name=Generator content="Microsoft Word 12 (filtered medium)"><!--[if !mso]><style>v\:* {behavior:url(#default#VML);}
o\:* {behavior:url(#default#VML);}
w\:* {behavior:url(#default#VML);}
.shape {behavior:url(#default#VML);}
</style><![endif]--><style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p
{mso-style-priority:99;
mso-margin-top-alt:auto;
margin-right:0cm;
mso-margin-bottom-alt:auto;
margin-left:0cm;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 70.85pt 70.85pt 70.85pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-US link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Andreas,<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Actually, since your raster coverage is tiled you should have used ST_Intersects(priogrid_land.cell, globshort.rast) from the beginning.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Understand that to do this query the system has to determine the value for each single pixel (they are 7 231 680 000 !!!!). Having precomputed statistics, like the histogram for each tile, would not help as your polygon grid do not necessarily correspond to the tile grid. Does it?<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>The real solution, I think, lies in using a lower resolution (overview) of the coverage to compute the statistics. The resulting numbers would be a bit different but not really far from the one you get at highest resolution. You would have to import your raster using the –l option and launch your query on the overview table instead.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>I’m looking for a solution actually counting pixel values but it is still slower than using ST_Intersection() (to my own surprise!). That means, for now, it is faster to convert each tile to a vector representation (that’a what ST_Intersection() is doing) and compute the areas of those polygons than to count and summarize pixel values. Work in progress…<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Pierre<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></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;font-family:"Tahoma","sans-serif"'>From:</span></b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> Andreas Forĝ Tollefsen [mailto:andreasft@gmail.com] <br><b>Sent:</b> 28 février 2011 04:48<br><b>To:</b> PostGIS Users Discussion<br><b>Cc:</b> Pierre Racine; Paragon Corporation<br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon<o:p></o:p></span></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I have reduced search time. It takes 534402 ms for 100 polygons now by using a && operator in the intersection query.<o:p></o:p></p><div><p class=MsoNormal>Having 64818 polygons this would take about 96 hours to complete.<o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><div><div><p class=MsoNormal>DROP TABLE IF EXISTS globshortrel;<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>SELECT <o:p></o:p></p></div><div><p class=MsoNormal>gid, <o:p></o:p></p></div><div><p class=MsoNormal>((foo.geomval).val),<o:p></o:p></p></div><div><p class=MsoNormal>CAST (SUM(ST_Area((foo.geomval).geom)) AS decimal(8,5)) as area,<o:p></o:p></p></div><div><p class=MsoNormal>CAST(SUM(ST_Area((foo.geomval).geom))/0.25*100 AS decimal(6,4)) as percentarea, <o:p></o:p></p></div><div><p class=MsoNormal>CAST (SUM(ST_Area((foo.geomval).geom))/0.0000077160617284 AS int) AS npixels<o:p></o:p></p></div><div><p class=MsoNormal>INTO globshortrel<o:p></o:p></p></div><div><p class=MsoNormal>FROM (SELECT priogrid_land.gid, ST_Intersection(globshort.rast, priogrid_land.cell) AS geomval FROM globshort, priogrid_land WHERE priogrid_land.cell && globshort.rast) AS foo<o:p></o:p></p></div><div><p class=MsoNormal>WHERE gid >= 139300 AND gid <= 139399<o:p></o:p></p></div><div><p class=MsoNormal>GROUP BY gid, (foo.geomval).val;<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Total query runtime: 534402 ms.<o:p></o:p></p></div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>2011/2/28 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com">andreasft@gmail.com</a>><o:p></o:p></p><p class=MsoNormal style='margin-bottom:12.0pt'>If I limit to only 10 polygons it uses 221285 ms.<o:p></o:p></p><div><p class=MsoNormal>2011/2/28 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>><o:p></o:p></p><div><div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm'><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Had the query running for over 48 hours but now it had crashed. Any idea on how to optimize for speed and stability?<o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><div><div><p class=MsoNormal>DROP TABLE IF EXISTS globshortrel;<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>SELECT <o:p></o:p></p></div><div><p class=MsoNormal>gid, <o:p></o:p></p></div><div><p class=MsoNormal>((foo.geomval).val),<o:p></o:p></p></div><div><p class=MsoNormal>CAST (SUM(ST_Area((foo.geomval).geom)) AS decimal(7,5)) as area, <o:p></o:p></p></div><div><p class=MsoNormal>CAST(SUM(ST_Area((foo.geomval).geom))/0.25*100 AS decimal(6,4)) as percentarea, <o:p></o:p></p></div><div><p class=MsoNormal>CAST (SUM(ST_Area((foo.geomval).geom))/0.0000077160617284 AS int) AS npixels<o:p></o:p></p></div><div><p class=MsoNormal>INTO globshortrel<o:p></o:p></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<o:p></o:p></p></div></div><div><div><p class=MsoNormal>GROUP BY gid, (foo.geomval).val<o:p></o:p></p></div></div><div><p class=MsoNormal>ORDER BY gid;<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>-- SELECT * FROM globshortrel;<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>-- SELECT DISTINCT ON(gid) gid, val, percentarea<o:p></o:p></p></div><div><p class=MsoNormal>-- FROM globshortrel<o:p></o:p></p></div><div><p class=MsoNormal>-- GROUP BY gid, val, percentarea<o:p></o:p></p></div><div><p class=MsoNormal>-- ORDER BY gid, percentarea DESC<o:p></o:p></p></div><div><p class=MsoNormal>ERROR: out of memory<o:p></o:p></p></div><div><p class=MsoNormal>DETAIL: Failed on request of size 1753647.<o:p></o:p></p></div><div><p class=MsoNormal>CONTEXT: SQL function "st_dumpaspolygons" statement 1<o:p></o:p></p></div><div><p class=MsoNormal>PL/pgSQL function "st_intersection" line 9 at RETURN QUERY<o:p></o:p></p></div><div><p class=MsoNormal>SQL function "st_intersection" statement 1<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>********** Error **********<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>ERROR: out of memory<o:p></o:p></p></div><div><p class=MsoNormal>SQL state: 53200<o:p></o:p></p></div><div><p class=MsoNormal>Detail: Failed on request of size 1753647.<o:p></o:p></p></div><div><p class=MsoNormal>Context: SQL function "st_dumpaspolygons" statement 1<o:p></o:p></p></div><div><p class=MsoNormal>PL/pgSQL function "st_intersection" line 9 at RETURN QUERY<o:p></o:p></p></div><div><p class=MsoNormal>SQL function "st_intersection" statement 1<o:p></o:p></p></div><div><div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>2011/2/25 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>><o:p></o:p></p><p class=MsoNormal>I wrote this query. It seems to give me the area, percentage of area and the pixelcount within a 0.5x0.5 cell.<o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p><div><div><p class=MsoNormal>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<o:p></o:p></p></div><div><p class=MsoNormal>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<o:p></o:p></p></div><div><p class=MsoNormal>WHERE gid = 139358<o:p></o:p></p></div><div><div><p class=MsoNormal>GROUP BY gid, ((foo.geomval).val)<o:p></o:p></p></div></div><div><p class=MsoNormal>ORDER BY area, ((foo.geomval).val)<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Result: <o:p></o:p></p></div><div><p class=MsoNormal>gid; area; percentarea; npixels<o:p></o:p></p></div><div><div><p class=MsoNormal>139358;0.000123456790106502;0.0493827160426008;15.9999999978027<o:p></o:p></p></div><div><p class=MsoNormal>139358;0.000563271604818283;0.225308641927313;72.9999999844495<o:p></o:p></p></div><div><p class=MsoNormal>139358;0.0966820987653989;38.6728395061596;12529.9999999957<o:p></o:p></p></div><div><p class=MsoNormal>139358;0.152631172839676;61.0524691358705;19781.0000000221<o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>2011/2/25 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>><o:p></o:p></p><div><div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm'><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Hi Pierre. I did tile it when i processed it.<o:p></o:p></p><div><p class=MsoNormal>I used this gdal2wktraster syntax: <span style='font-size:9.0pt;color:black'>python gdal2wktraster.py -r c:\prio_grid\source\globcover\globshort.tif -t globshort -s 4326 -k 720x360 -I -o globshort.sql</span><o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal style='margin-bottom:12.0pt'><span style='font-size:9.0pt'>The size of my raster is X: 129600 Y: 55800</span><o:p></o:p></p><div><p class=MsoNormal>2011/2/24 Pierre Racine <<a href="mailto:Pierre.Racine@sbf.ulaval.ca" target="_blank">Pierre.Racine@sbf.ulaval.ca</a>><o:p></o:p></p><div><div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm'><p class=MsoNormal><o:p> </o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:11.0pt;color:#1F497D'> </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:11.0pt;color:#1F497D'> </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:11.0pt;color:#1F497D'>Pierre</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:11.0pt;color:#1F497D'> </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:11.0pt;color:#1F497D'> </span><o:p></o:p></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 style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><div><div><p class=MsoNormal><br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon<o:p></o:p></p></div></div></div></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Great.<o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Thank you so much. I should have noticed that these were unioned numbers. <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Btw. are there any way of getting the count of pixels within a polygon without actually aggregating them or union them?<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;margin-bottom:12.0pt'>Andreas<o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>2011/2/24 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>><o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>Andreas,</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>So </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>SELECT gid, SUM(ST_Area((foo.geomval).geom))/ [put your pixel area size here] as ct<o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>WHERE gid >= 139358 AND gid <= 139365<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>GROUP BY gid<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>ORDER BY gid<o:p></o:p></p></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div class=MsoNormal align=center style='text-align:center'><hr size=2 width="100%" align=center></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt'><br><b>Subject:</b> Re: [postgis-users] ST_Value from Polygon</span><o:p></o:p></p></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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. <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>What i get is different, and it seems that the query is not providing me with the number of pixels within the grid cell.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Any idea why this is so different?<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>SELECT gid, count((foo.geomval).val) as ct<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>WHERE gid >= 139358 AND gid <= 139365<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>GROUP BY gid<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>ORDER BY gid<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Result:<o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139358;632<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139359;1030<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139360;912<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139361;731<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139362;760<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139363;1230<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139364;1314<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>139365;1014<o:p></o:p></p></div></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>The attached image shows the raster pixels within one cell.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>2011/2/24 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Thanks! <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>That solved it.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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. <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='color:#888888'>Andreas</span><o:p></o:p></p></div><div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>2011/2/23 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>><o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>Andrea,</span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>Try </span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='color:blue'>WHERE gid > 151000 AND gid < 151010</span><o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='color:blue'>GROUP BY gid, (foo.geomval).val</span><o:p></o:p></p></div></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='color:blue'>ORDER BY gid, ct DESC</span><o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div class=MsoNormal align=center style='text-align:center'><hr size=2 width="100%" align=center></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p></div><p class=MsoNormal style='mso-margin-top-alt:auto;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><o:p></o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Hi. Thanks Regina and Leo, <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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. <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>So my result will be the rastervalue with the highest count within each of the grid cells.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>However, as far as i know, there is now Max(COUNT) function in postgresql.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Any idea how i can modify the below query to only return the rastervalue within the grid cell occuring most frequently?<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Consequently i want only one row for each gid, and the maximum occuring rastervalue.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>SELECT gid, (foo.geomval).val, COUNT((foo.geomval).val) AS ct<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>WHERE gid > 151000 AND gid < 151010<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>GROUP BY gid, (foo.geomval).val;<o:p></o:p></p></div></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>gid; val; ct<o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;14;381<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;150;9<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;50;7<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;140;91<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;40;1<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;70;2<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;130;4<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;200;48<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;100;3<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;;0<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;190;1<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;20;203<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;11;111<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;210;16<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>151001;30;105<o:p></o:p></p></div></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>2011/2/23 Paragon Corporation <<a href="mailto:lr@pcorp.us" target="_blank">lr@pcorp.us</a>><o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><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><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='color:blue'><a href="http://www.postgis.us/presentations" target="_blank">http://www.postgis.us/presentations</a></span><o:p></o:p></p><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>Hope that helps,</span><o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'><span style='font-size:10.0pt;color:blue'>Regina and Leo</span><o:p></o:p></p></div><div class=MsoNormal align=center style='text-align:center'><hr size=2 width="100%" align=center></div><p class=MsoNormal style='mso-margin-top-alt:auto;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><o:p></o:p></p><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Hi all, <o:p></o:p></p><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>I am working with a large raster dataset that i want to aggregate into vector grids.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>I have been looking at the ST_Value function, but this is not usable together with the cell polygon.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>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.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Hence i need to calculate the area of each raster class within each cell and select the largest class.<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>Any idea? So far i have only come this far:<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div><div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>DROP TABLE IF EXISTS globshortpoly;<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>SELECT priogrid_land.cell, ST_Value(rast, ST_Centroid(cell))<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>INTO globshortpoly<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>FROM priogrid_land, globshort<o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>WHERE rast && priogrid_land.cell <o:p></o:p></p></div><div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'>LIMIT 1000<o:p></o:p></p></div></div></div></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;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><o:p></o:p></p></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div></div></div></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;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><o:p></o:p></p></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div></div></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div></div></div></div></div><p class=MsoNormal style='mso-margin-top-alt:auto;mso-margin-bottom-alt:auto'> <o:p></o:p></p></div></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><o:p></o:p></p></blockquote></div></div></div><p class=MsoNormal><o:p> </o:p></p></div></blockquote></div></div></div><p class=MsoNormal><o:p> </o:p></p></div></div></div><p class=MsoNormal><o:p> </o:p></p></div></div></div></blockquote></div></div></div><p class=MsoNormal><o:p> </o:p></p></div><p class=MsoNormal><o:p> </o:p></p></div></div></div></div></body></html>