Frank,<div><br></div><div>Thanks for tip. I was able to rasterize my polygon and used that to calculate my sum. Very fast. </div><div><br></div><div>bw<br><br><div class="gmail_quote">On Fri, Dec 3, 2010 at 2:28 PM, Frank Warmerdam <span dir="ltr"><<a href="mailto:warmerdam@pobox.com">warmerdam@pobox.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div class="im">Brian Walawender wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hello,<br>
<br>
I am looking for some assistance in optimizing a section of code for faster performance. Here is my problem, I have GeoTiff that contains 1 km x 1 km population data over an area roughly the size of the continental US (7020 x 3000). I am trying to calculate the population within a polygon. I am doing this by finding the extent of the polygon and reading in that section of the GeoTiff using the ReadAsArray function. At this point I can quickly calculate the population within the extent by using the numpy sum. However, the only way I can figure to sum the points within the polygon is to iterate over the array checking each point using ogr. If the polygon is very large, this can take an extended period of time. Is there a faster way to calculate the population within a large polygon? My code and some sample output is below.<br>
</blockquote>
<br></div>
Brian,<br>
<br>
I think you would be better off rasterizing your mask polygon and operating<br>
on the masked raster. It might be sufficient to set all cells outside the<br>
region to 0 in your target raster if you are just summing. Or you might<br>
generate the mask separately and just check it instead of doing an<br>
expensive point in polygon test.<br>
<br>
The test script for the rasterize function is available at the following<br>
url and might provide a hint of how to use it.<br>
<br>
<a href="http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py" target="_blank">http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py</a><br>
<br>
You can also find more about the rasterize api from it's C++ docs at:<br>
<br>
<a href="http://www.gdal.org/gdal__alg_8h.html#50caf4bc34703f0bcf515ecbe5061a0a" target="_blank">http://www.gdal.org/gdal__alg_8h.html#50caf4bc34703f0bcf515ecbe5061a0a</a><br>
<br>
Good luck,<br><font color="#888888">
-- <br>
---------------------------------------+--------------------------------------<br>
I set the clouds in motion - turn up | Frank Warmerdam, <a href="mailto:warmerdam@pobox.com" target="_blank">warmerdam@pobox.com</a><br>
light and sound - activate the windows | <a href="http://pobox.com/~warmerdam" target="_blank">http://pobox.com/~warmerdam</a><br>
and watch the world go round - Rush | Geospatial Programmer for Rent<br>
<br>
</font></blockquote></div><br><br clear="all"><br>-- <br><span style="font-family:arial, sans-serif;font-size:13px;border-collapse:collapse"><p style="margin-top:0px;margin-right:0px;margin-bottom:0px;margin-left:0px"><span style="font-size:10pt;color:gray">Brian Walawender</span></p>
<p style="margin-top:0px;margin-right:0px;margin-bottom:0px;margin-left:0px"><span style="font-size:10pt;color:gray">Technique Development Meteorologist</span></p><p style="margin-top:0px;margin-right:0px;margin-bottom:0px;margin-left:0px">
<span style="font-size:10pt;color:gray">Scientific Services Division - Central Region Headquarters</span></p><p style="margin-top:0px;margin-right:0px;margin-bottom:0px;margin-left:0px"><span style="font-size:10pt;color:gray">816-268-3114 - Office</span></p>
<p style="margin-top:0px;margin-right:0px;margin-bottom:0px;margin-left:0px"><span style="font-size:10pt;color:gray">816-805-6497 - Cell</span></p></span><div><br></div><br>
</div>