Hello,<div><br></div><div>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.</div>
<div><br></div><div>Thanks.</div><div><br></div><div>bw</div><div><br></div><div>#findPopulationOGR.py - read a polygon from the database</div><div># and calculate the population within it</div><div><br></div><div>from osgeo import gdal</div>
<div>from osgeo import ogr</div><div>from osgeo import osr</div><div>from osgeo.gdalconst import *</div><div>from osgeo.gdal_array import *</div><div>import sys, time</div><div>import numpy as np</div><div><br></div><div>
populationFile = "pop4269.tif"</div><div>band = 1</div><div># Zone based</div><div>#segment_id = 147231</div><div># Polygon based</div><div>segment_id = 147810</div><div><br></div><div>sql = "SELECT AsBinary(ST_MemUnion(COALESCE (ppoly.geom, pl.geom))) AS wkb_geometry " </div>
<div>sql = sql + "FROM iris.product_segment ps "</div><div>sql = sql + "JOIN iris.product_segment_to_product_location pspl ON <a href="http://ps.id" target="_blank">ps.id</a> = pspl.product_segment_id "</div>
<div>sql = sql + "JOIN iris.product_location pl ON pspl.product_location_id = <a href="http://pl.id" target="_blank">pl.id</a> "</div>
<div>sql = sql + "JOIN iris.product_location_type plt ON pl.product_location_type_id = <a href="http://plt.id" target="_blank">plt.id</a> "</div><div>sql = sql + "LEFT JOIN iris.product_polygon ppoly ON <a href="http://ps.id" target="_blank">ps.id</a> = ppoly.product_segment_id "</div>
<div>sql = sql + " WHERE <a href="http://ps.id" target="_blank">ps.id</a>=" + `segment_id`</div><div><br></div><div>conn = ogr.Open("PG: host='localhost' dbname=mydb' ")</div><div>layer = conn.ExecuteSQL(sql)</div>
<div>feature = layer.GetFeature(0)</div><div>geom = feature.GetGeometryRef()</div><div><br></div><div>#create memory layer for polygon</div><div>srWGS84 = osr.SpatialReference()</div><div>srWGS84.ImportFromProj4('+proj=longlat +datum=WGS84')</div>
<div>polyds = ogr.GetDriverByName('Memory').CreateDataSource('')</div><div>polylyr = polyds.CreateLayer('poly',geom_type=ogr.wkbPolygon, srs=srWGS84)</div><div>feat = ogr.Feature(polylyr.GetLayerDefn())</div>
<div>feat.SetGeometryDirectly(geom)</div><div>polylyr.CreateFeature(feat)</div><div><br></div><div>#check polygon layer extent</div><div>extent = polylyr.GetExtent()</div><div><br></div><div><br></div><div>print '========Polygon Bounding Box========'</div>
<div>print 'UL: ', extent[0], extent[3]</div><div>print 'LR: ', extent[1], extent[2]</div><div><br></div><div>XUL = extent[0]</div><div>YUL = extent[3]</div><div>XLR = extent[1]</div><div>YLR = extent[2]</div>
<div><br></div><div>gdal.AllRegister()</div><div><br></div><div>ds = gdal.Open( populationFile, GA_ReadOnly )</div><div><br></div><div>if ds is None:</div><div> print('Cannot open %s' % populationFile)</div><div>
sys.exit(1)</div><div> </div><div>cols = ds.RasterXSize</div><div>rows = ds.RasterYSize</div><div>proj = ds.GetProjection()</div><div>geomatrix=ds.GetGeoTransform()</div><div>(success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)</div>
<div>xUL = int(inv_geometrix[0] + inv_geometrix[1] * XUL + inv_geometrix[2] * YUL)</div><div>yUL = int(inv_geometrix[3] + inv_geometrix[4] * XUL + inv_geometrix[5] * YUL)</div><div>xLR = int(inv_geometrix[0] + inv_geometrix[1] * XLR + inv_geometrix[2] * YLR)</div>
<div>yLR = int(inv_geometrix[3] + inv_geometrix[4] * XLR + inv_geometrix[5] * YLR)</div><div>xdif = xLR - xUL</div><div>ydif = yUL - yLR</div><div><br></div><div>#print cols, rows, proj, xUL, yUL, xLR, yLR, xdif, ydif</div>
<div><br></div><div>data = ds.GetRasterBand(band).ReadAsArray(xUL, yLR, xdif, ydif)</div><div><br></div><div>extentPop = np.sum(data)</div><div><br></div><div>print "Print total population in extent = %#10d" % (extentPop )</div>
<div><br></div><div>startTime = time.time()</div><div><br></div><div>for i in range(data.shape[0]):</div><div> for j in range(data.shape[1]):</div><div> X = i + xUL</div><div> Y = j + yLR</div><div> lon = geomatrix[0] + X*geomatrix[1] + Y*geomatrix[2] </div>
<div> lat = geomatrix[3] + X*geomatrix[4] + Y*geomatrix[5]</div><div> qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % (lon, lat))</div><div> if geom.Contains(qpt) == False:</div><div> data[i,j] = 0</div>
<div> </div><div>polyPop = np.sum(data)</div><div><br></div><div>endTime = time.time()</div><div> </div><div>print "Print total population in polygon = %#10d" % (polyPop )</div><div><br></div>
<div>print "The calc took " + str(endTime-startTime) + " seconds"</div><div><br></div><div>ds = None</div><div>conn = None</div><div>polylyr = None</div><div><br></div><div><div>$ python findPopulationOGR.py </div>
<div>========Polygon Bounding Box========</div><div>UL: -108.891024 49.000022</div><div>LR: -104.041596 46.540403</div><div>Print total population in extent = 101602</div><div>Print total population in polygon = 21052</div>
<div>The calc took 604.560879946 seconds</div></div><div> <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>