<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<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:x="urn:schemas-microsoft-com:office:excel" xmlns:p="urn:schemas-microsoft-com:office:powerpoint" xmlns:a="urn:schemas-microsoft-com:office:access" xmlns:dt="uuid:C2F41010-65B3-11d1-A29F-00AA00C14882" xmlns:s="uuid:BDC6E3F0-6DA3-11d1-A2A3-00AA00C14882" xmlns:rs="urn:schemas-microsoft-com:rowset" xmlns:z="#RowsetSchema" xmlns:b="urn:schemas-microsoft-com:office:publisher" xmlns:ss="urn:schemas-microsoft-com:office:spreadsheet" xmlns:c="urn:schemas-microsoft-com:office:component:spreadsheet" xmlns:odc="urn:schemas-microsoft-com:office:odc" xmlns:oa="urn:schemas-microsoft-com:office:activation" xmlns:html="http://www.w3.org/TR/REC-html40" xmlns:q="http://schemas.xmlsoap.org/soap/envelope/" xmlns:rtc="http://microsoft.com/officenet/conferencing" xmlns:D="DAV:" xmlns:Repl="http://schemas.microsoft.com/repl/" xmlns:mt="http://schemas.microsoft.com/sharepoint/soap/meetings/" xmlns:x2="http://schemas.microsoft.com/office/excel/2003/xml" xmlns:ppda="http://www.passport.com/NameSpace.xsd" xmlns:ois="http://schemas.microsoft.com/sharepoint/soap/ois/" xmlns:dir="http://schemas.microsoft.com/sharepoint/soap/directory/" xmlns:ds="http://www.w3.org/2000/09/xmldsig#" xmlns:dsp="http://schemas.microsoft.com/sharepoint/dsp" xmlns:udc="http://schemas.microsoft.com/data/udc" xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:sub="http://schemas.microsoft.com/sharepoint/soap/2002/1/alerts/" xmlns:ec="http://www.w3.org/2001/04/xmlenc#" xmlns:sp="http://schemas.microsoft.com/sharepoint/" xmlns:sps="http://schemas.microsoft.com/sharepoint/soap/" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:udcs="http://schemas.microsoft.com/data/udc/soap" xmlns:udcxf="http://schemas.microsoft.com/data/udc/xmlfile" xmlns:udcp2p="http://schemas.microsoft.com/data/udc/parttopart" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns:st="" xmlns="http://www.w3.org/TR/REC-html40"><head><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;}
span.EmailStyle17
{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'>You should just divide the area of your polygon by the area of one of your pixel.<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 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?<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><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"'> postgis-users-bounces@postgis.refractions.net [mailto:postgis-users-bounces@postgis.refractions.net] <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<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>Great.<o:p></o:p></p><div><p class=MsoNormal>Thank you so much. I should have noticed that these were unioned numbers. <o:p></o:p></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?<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'>Andreas<o:p></o:p></p><div><p class=MsoNormal>2011/2/24 Paragon Corporation <<a href="mailto:lr@pcorp.us">lr@pcorp.us</a>><o:p></o:p></p><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";color:blue'>Andreas,</span><o:p></o:p></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";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> <o:p></o:p></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";color:blue'>So </span><o:p></o:p></p><p class=MsoNormal> <o:p></o:p></p><div><p class=MsoNormal>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>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>WHERE gid >= 139358 AND gid <= 139365<o:p></o:p></p></div><div><p class=MsoNormal>GROUP BY gid<o:p></o:p></p></div><div><p class=MsoNormal>ORDER BY gid<o:p></o:p></p></div></div><p class=MsoNormal><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><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:<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<o:p></o:p></span></p><div><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'><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><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. <o:p></o:p></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.<o:p></o:p></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.<o:p></o:p></p></div><div><p class=MsoNormal>Any idea why this is so different?<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><div><p class=MsoNormal>SELECT gid, count((foo.geomval).val) as ct<o:p></o:p></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<o:p></o:p></p></div><div><p class=MsoNormal>WHERE gid >= 139358 AND gid <= 139365<o:p></o:p></p></div><div><p class=MsoNormal>GROUP BY gid<o:p></o:p></p></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>Result:<o:p></o:p></p></div><div><div><p class=MsoNormal>139358;632<o:p></o:p></p></div><div><p class=MsoNormal>139359;1030<o:p></o:p></p></div><div><p class=MsoNormal>139360;912<o:p></o:p></p></div><div><p class=MsoNormal>139361;731<o:p></o:p></p></div><div><p class=MsoNormal>139362;760<o:p></o:p></p></div><div><p class=MsoNormal>139363;1230<o:p></o:p></p></div><div><p class=MsoNormal>139364;1314<o:p></o:p></p></div><div><p class=MsoNormal>139365;1014<o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>The attached image shows the raster pixels within one cell.<o:p></o:p></p></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/24 Andreas Forĝ Tollefsen <<a href="mailto:andreasft@gmail.com" target="_blank">andreasft@gmail.com</a>><o:p></o:p></p><p class=MsoNormal>Thanks! <o:p></o:p></p><div><p class=MsoNormal>That solved it.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></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. <o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal><span style='color:#888888'>Andreas<o:p></o:p></span></p></div><div><div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>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><span style='font-size:10.0pt;font-family:"Arial","sans-serif";color:blue'>Andrea,</span><o:p></o:p></p><p class=MsoNormal> <o:p></o:p></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";color:blue'>Try </span><o:p></o:p></p><p class=MsoNormal> <o:p></o:p></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";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><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><span style='color:blue'>WHERE gid > 151000 AND gid < 151010</span><o:p></o:p></p></div><div><p class=MsoNormal><span style='color:blue'>GROUP BY gid, (foo.geomval).val</span><o:p></o:p></p></div></div><div><p class=MsoNormal><span style='color:blue'>ORDER BY gid, ct DESC</span><o:p></o:p></p></div><div><p class=MsoNormal> <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><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"'> <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<o:p></o:p></span></p></div><p class=MsoNormal style='margin-bottom:12.0pt'><b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>Sent:</span></b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> 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>Hi. Thanks Regina and Leo, <o:p></o:p></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. <o:p></o:p></p><div><p class=MsoNormal>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>However, as far as i know, there is now Max(COUNT) function in postgresql.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></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?<o:p></o:p></p></div><div><p class=MsoNormal>Consequently i want only one row for each gid, and the maximum occuring rastervalue.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><div><p class=MsoNormal>SELECT gid, (foo.geomval).val, COUNT((foo.geomval).val) AS ct<o:p></o:p></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<o:p></o:p></p></div><div><p class=MsoNormal>WHERE gid > 151000 AND gid < 151010<o:p></o:p></p></div><div><p class=MsoNormal>GROUP BY gid, (foo.geomval).val;<o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>gid; val; ct<o:p></o:p></p></div><div><div><p class=MsoNormal>151001;14;381<o:p></o:p></p></div><div><p class=MsoNormal>151001;150;9<o:p></o:p></p></div><div><p class=MsoNormal>151001;50;7<o:p></o:p></p></div><div><p class=MsoNormal>151001;140;91<o:p></o:p></p></div><div><p class=MsoNormal>151001;40;1<o:p></o:p></p></div><div><p class=MsoNormal>151001;70;2<o:p></o:p></p></div><div><p class=MsoNormal>151001;130;4<o:p></o:p></p></div><div><p class=MsoNormal>151001;200;48<o:p></o:p></p></div><div><p class=MsoNormal>151001;100;3<o:p></o:p></p></div><div><p class=MsoNormal>151001;;0<o:p></o:p></p></div><div><p class=MsoNormal>151001;190;1<o:p></o:p></p></div><div><p class=MsoNormal>151001;20;203<o:p></o:p></p></div><div><p class=MsoNormal>151001;11;111<o:p></o:p></p></div><div><p class=MsoNormal>151001;210;16<o:p></o:p></p></div><div><p class=MsoNormal>151001;30;105<o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>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><span style='font-size:10.0pt;font-family:"Arial","sans-serif";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> <o:p></o:p></p><p class=MsoNormal> <o:p></o:p></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><o:p></o:p></p><p class=MsoNormal> <o:p></o:p></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";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> <o:p></o:p></p><p class=MsoNormal><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> <o:p></o:p></p><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";color:blue'>Hope that helps,</span><o:p></o:p></p></div><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Arial","sans-serif";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='margin-bottom:12.0pt'><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"'> <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>Hi all, <o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>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>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><o:p> </o:p></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.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></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.<o:p></o:p></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.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Any idea? So far i have only come this far:<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><div><p class=MsoNormal>DROP TABLE IF EXISTS globshortpoly;<o:p></o:p></p></div><div><p class=MsoNormal>SELECT priogrid_land.cell, ST_Value(rast, ST_Centroid(cell))<o:p></o:p></p></div><div><p class=MsoNormal>INTO globshortpoly<o:p></o:p></p></div><div><p class=MsoNormal>FROM priogrid_land, globshort<o:p></o:p></p></div><div><p class=MsoNormal>WHERE rast && priogrid_land.cell <o:p></o:p></p></div><div><p class=MsoNormal>LIMIT 1000<o:p></o:p></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><o:p></o:p></p></div><p class=MsoNormal><o:p> </o:p></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><o:p></o:p></p></div><p class=MsoNormal><o:p> </o:p></p></div></div></div></div><p class=MsoNormal><o:p> </o:p></p></div></div></div></div></div><p class=MsoNormal><o:p> </o:p></p></div></div></div></body></html>