[gdal-dev] GDAL/OGR Python Question

[Brian Walawender]

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.

Thanks.

bw

#findPopulationOGR.py - read a polygon from the database
# and calculate the population within it

from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo.gdalconst import *
from osgeo.gdal_array import *
import sys, time
import numpy as np

populationFile = "pop4269.tif"
band = 1
# Zone based
#segment_id = 147231
# Polygon based
segment_id = 147810

sql = "SELECT AsBinary(ST_MemUnion(COALESCE (ppoly.geom, pl.geom))) AS
wkb_geometry "
sql = sql + "FROM iris.product_segment ps "
sql = sql + "JOIN iris.product_segment_to_product_location pspl ON ps.id =
pspl.product_segment_id "
sql = sql + "JOIN iris.product_location pl ON pspl.product_location_id =
pl.id "
sql = sql + "JOIN iris.product_location_type plt ON
pl.product_location_type_id = plt.id "
sql = sql + "LEFT JOIN iris.product_polygon ppoly ON ps.id =
ppoly.product_segment_id "
sql = sql + " WHERE ps.id="  + `segment_id`

conn = ogr.Open("PG: host='localhost' dbname=mydb' ")
layer = conn.ExecuteSQL(sql)
feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()

#create memory layer for polygon
srWGS84 = osr.SpatialReference()
srWGS84.ImportFromProj4('+proj=longlat +datum=WGS84')
polyds = ogr.GetDriverByName('Memory').CreateDataSource('')
polylyr = polyds.CreateLayer('poly',geom_type=ogr.wkbPolygon, srs=srWGS84)
feat = ogr.Feature(polylyr.GetLayerDefn())
feat.SetGeometryDirectly(geom)
polylyr.CreateFeature(feat)

#check polygon layer extent
extent = polylyr.GetExtent()


print '========Polygon Bounding Box========'
print 'UL: ', extent[0], extent[3]
print 'LR: ', extent[1], extent[2]

XUL = extent[0]
YUL = extent[3]
XLR = extent[1]
YLR = extent[2]

gdal.AllRegister()

ds = gdal.Open( populationFile, GA_ReadOnly )

if ds is None:
    print('Cannot open %s' % populationFile)
    sys.exit(1)

cols = ds.RasterXSize
rows = ds.RasterYSize
proj = ds.GetProjection()
geomatrix=ds.GetGeoTransform()
(success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)
xUL = int(inv_geometrix[0] + inv_geometrix[1] * XUL + inv_geometrix[2] *
YUL)
yUL = int(inv_geometrix[3] + inv_geometrix[4] * XUL + inv_geometrix[5] *
YUL)
xLR = int(inv_geometrix[0] + inv_geometrix[1] * XLR + inv_geometrix[2] *
YLR)
yLR = int(inv_geometrix[3] + inv_geometrix[4] * XLR + inv_geometrix[5] *
YLR)
xdif = xLR - xUL
ydif = yUL - yLR

#print cols, rows, proj, xUL, yUL, xLR, yLR, xdif, ydif

data = ds.GetRasterBand(band).ReadAsArray(xUL, yLR, xdif, ydif)

extentPop = np.sum(data)

print "Print total population in extent = %#10d" % (extentPop )

startTime = time.time()

for i in range(data.shape[0]):
    for j in range(data.shape[1]):
        X = i + xUL
        Y = j + yLR
        lon = geomatrix[0] + X*geomatrix[1] + Y*geomatrix[2]
        lat = geomatrix[3] + X*geomatrix[4] + Y*geomatrix[5]
        qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % (lon, lat))
        if geom.Contains(qpt) == False:
            data[i,j] = 0

polyPop = np.sum(data)

endTime = time.time()

print "Print total population in polygon = %#10d" % (polyPop )

print "The calc took " + str(endTime-startTime) + " seconds"

ds = None
conn = None
polylyr = None

$ python findPopulationOGR.py
========Polygon Bounding Box========
UL:  -108.891024 49.000022
LR:  -104.041596 46.540403
Print total population in extent =     101602
Print total population in polygon =      21052
The calc took 604.560879946 seconds


-- 

Brian Walawender

Technique Development Meteorologist

Scientific Services Division - Central Region Headquarters

816-268-3114 - Office

816-805-6497 - Cell
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20101203/6e62f869/attachment.html