[gdal-dev] GDAL/OGR Python Question

[Frank Warmerdam]

Brian Walawender wrote:
> Hello,
> 
> 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.

Brian,

I think you would be better off rasterizing your mask polygon and operating
on the masked raster.  It might be sufficient to set all cells outside the
region to 0 in your target raster if you are just summing.  Or you might
generate the mask separately and just check it instead of doing an
expensive point in polygon test.

The test script for the rasterize function is available at the following
url and might provide a hint of how to use it.

   http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

You can also find more about the rasterize api from it's C++ docs at:

   http://www.gdal.org/gdal__alg_8h.html#50caf4bc34703f0bcf515ecbe5061a0a

Good luck,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent