Ok. I was worried there was something genious about your suggestion that I could not grasp :)<div><br></div><div>Andreas<br><br><div class="gmail_quote">2013/6/24 Kim Bisgaard <span dir="ltr"><<a href="mailto:kib@dmi.dk" target="_blank">kib@dmi.dk</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
Hi,<br>
<br>
You are right, I mistook ST_intersects() for ST_intersection() - me butting out of this thread..<br>
<br>
Regards,<br>
Kim<br>
<br>
<br>
<div>On 2013-06-23 10:46, Andreas Forų Tollefsen wrote:<br>
</div>
<blockquote type="cite">
Thanks for your answers.
<div><div class="h5"><div><br>
</div>
<br>
<div>Kim: We do not directly use ST_Intersection in our script, but ST_Clip to clip the raster according to our polygons, and only where they ST_Intersects(), not ST_Intersection(); </div>
<div>I still do not understand your suggestion.</div>
<div><br>
</div>
<div>Best,</div>
<div>Andreas<br>
<br>
<div class="gmail_quote">2013/6/18 Kim Bisgaard <span dir="ltr"><<a href="mailto:kib@dmi.dk" target="_blank">kib@dmi.dk</a>></span><br>
<blockquote class="gmail_quote">
<div>Hi,<br>
<br>
Because ST_Intersection() returns neighbour pixels sharing a same value as only one<br>
polygon. You thus get a bigger area and thus more points.<br>
<br>
Bit me once :-/<br>
<br>
Regards,<br>
Kim
<div>
<div><br>
<br>
<br>
<br>
<br>
<div>On 2013-06-18 14:59, Andreas Forų Tollefsen wrote:<br>
</div>
<blockquote type="cite">Hi Kim,
<div><br>
</div>
<div>Thanks for your answer. However, we want raster as an output, since we want to be able to use the summarystats function.</div>
<div>Please elaborate how you think ST_PixelAsPolygons should solve out issue?</div>
<div>Thanks.</div>
<div><br>
</div>
<div>Andreas<br>
<br>
<div class="gmail_quote">2013/6/18 Kim Bisgaard <span dir="ltr"><<a href="mailto:kib@dmi.dk" target="_blank">kib@dmi.dk</a>></span><br>
<blockquote class="gmail_quote">
<div>Hi,<br>
<br>
Try to use 'ST_PixelAsPolygons(ST_Clip(n.rast, p.cell))'<br>
instead of 'ST_Intersects(n.rast, p.cell)'<br>
<br>
Regards,<br>
Kim
<div>
<div><br>
<br>
<div>On 2013-06-18 11:03, Andreas Forų Tollefsen wrote:<br>
</div>
</div>
</div>
<blockquote type="cite">
<div>
<div>Hi,
<div><br>
</div>
<div>We are working on a raster summarystats script to calculate various statistics for the pixels within fishnet polygons.</div>
<div><br>
</div>
<div>Our raster cell size is 0.0083333333333... x 0.0083333333333... degrees while our quadrat polygons are 0.5 x 0.5 decimal degrees.</div>
<div>This should give us 60x60 raster pixels within each of our polygons. ArcGIS zonal statistics returns a pixel count of 3600 in addition to other statistics.</div>
<div>However, PostGIS returns 3721 pixel count.</div>
<div><br>
</div>
<div>We do not really understand why, but it seems that our query includes some pixels that are outside of the polygon, but still touches the vertices of the polygon and are therefore included in the calculation.</div>
<div>Are there any way of modifying our script to return the same result as ArcGIS?</div>
<div>Thanks!</div>
<div><br>
</div>
<div>Andreas</div>
<div><br>
</div>
<div>script:</div>
<div><br>
</div>
<div>
<div>/* This query makes one raster for each PRIO-GRID cell. Clip and union is the procedure. */</div>
<div>INSERT INTO nightlightsprio (gid, "year", rast)</div>
<div>(SELECT gid, "year", ST_Union(raster) as rast</div>
<div>FROM </div>
<div>(SELECT p.gid, n."year", ST_Clip(n.rast, p.cell) as raster </div>
<div>FROM nightlights n, priogridyear p</div>
<div>WHERE ST_Intersects(n.rast, p.cell)</div>
<div>AND n."year" = p."year"</div>
<div>) </div>
<div>as priorast</div>
<div>GROUP BY gid, "year");</div>
<div><br>
</div>
<div><br>
</div>
<div>/* Default BandNoDataValue is 0. Raster value 0 means no light, not no data. Setting to NULL. This produces correct results. */</div>
<div>UPDATE nightlightsprio2 SET rast = ST_SetBandNoDataValue(rast, 1, NULL);</div>
<div><br>
</div>
<div><br>
</div>
<div>ALTER TABLE nightlightsprio2 ADD COLUMN nightlights_sum double precision,</div>
<div><span></span>ADD COLUMN nightlights_mean double precision,</div>
<div><span></span>ADD COLUMN nightlights_sd double precision,</div>
<div><span></span>ADD COLUMN nightlights_min double precision, </div>
<div><span></span>ADD COLUMN nightlights_max double precision, </div>
<div><span></span>ADD COLUMN nightlights_count integer;</div>
<div><br>
</div>
<div>UPDATE nightlightsprio2 SET nightlights_sum = (ST_SummaryStats(rast)).sum;</div>
<div>UPDATE nightlightsprio2 SET nightlights_mean = (ST_SummaryStats(rast)).mean;</div>
<div>UPDATE nightlightsprio2 SET nightlights_sd = (ST_SummaryStats(rast)).stddev;</div>
<div>UPDATE nightlightsprio2 SET nightlights_min = (ST_SummaryStats(rast)).min;</div>
<div>UPDATE nightlightsprio2 SET nightlights_max = (ST_SummaryStats(rast)).max;</div>
<div>UPDATE nightlightsprio2 SET nightlights_count = (ST_SummaryStats(rast)).count;</div>
</div>
<div><br>
</div>
<br>
<fieldset></fieldset> <br>
</div>
</div>
<pre>_______________________________________________
postgis-users mailing list
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><span>
</span></pre>
<span></span></blockquote>
<span><br>
<pre cols="140">--
Kim Bisgaard
Application Development Division Phone: <a href="tel:%2B45%203915%207562" value="+4539157562" target="_blank">+45 3915 7562</a> (direct)
Danish Meteorological Institute Fax: <a href="tel:%2B45%203915%207460" value="+4539157460" target="_blank">+45 3915 7460</a> (division)</pre>
</span></div>
<br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
<br>
</blockquote>
</div>
<br>
</div>
<br>
<fieldset></fieldset> <br>
<pre>_______________________________________________
postgis-users mailing list
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a>
</pre>
</blockquote>
<br>
<pre cols="140">--
Kim Bisgaard
Application Development Division Phone: <a href="tel:%2B45%203915%207562" value="+4539157562" target="_blank">+45 3915 7562</a> (direct)
Danish Meteorological Institute Fax: <a href="tel:%2B45%203915%207460" value="+4539157460" target="_blank">+45 3915 7460</a> (division)</pre>
</div>
</div>
</div>
<br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
<br>
</blockquote>
</div>
<br>
</div>
<br>
<fieldset></fieldset> <br>
<pre>_______________________________________________
postgis-users mailing list
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a>
</pre>
</div></div></blockquote><div><div class="h5">
<br>
<pre cols="140">--
Kim Bisgaard
Application Development Division Phone: <a href="tel:%2B45%203915%207562" value="+4539157562" target="_blank">+45 3915 7562</a> (direct)
Danish Meteorological Institute Fax: <a href="tel:%2B45%203915%207460" value="+4539157460" target="_blank">+45 3915 7460</a> (division)</pre>
</div></div></div>
<br>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
<br></blockquote></div><br></div>