Hi,<br>I am using the Python OGR bindings to create a shapefile with a regularly spaced grid over some multipolygons. I only want to have grid cells when the centroid of the cell falls within one of my multipolygons. My code is simple and rustic, but seems to be working. However, it is taking a long time to run, I suspect due to the large number of nested loops. I am attaching the code. The polygons in world.shp are fairly complex, but I was just wondering whether there might be something really obvious here. As for versions, OGR/GDAL is 1.5.3 and it is compiled with GEOS <a href="http://2.2.3.">2.2.3.</a> <br>
<br>#!/usr/bin/python <br>from osgeo import ogr<br>from osgeo import osr<br><br>g_world = ogr.Open(&quot;world.shp&quot;)<br>Lworld = g_world.GetLayer(0)<br>feat_world = Lworld.GetFeature(0)<br>geom_world = feat_world.GetGeometryRef()<br>
<br>t_srs = osr.SpatialReference()<br>t_srs.SetFromUserInput(&#39;WGS84&#39;)<br>drv = ogr.GetDriverByName(&#39;ESRI Shapefile&#39;)<br>ds = drv.CreateDataSource(&#39;/tmp/output_grid.shp&#39;)<br>layer = ds.CreateLayer(ds.GetName(), geom_type = ogr.wkbPolygon, srs = t_srs)<br>
layer.CreateField(ogr.FieldDefn(&#39;id&#39;, ogr.OFTInteger))<br><br>x_res = 0.5<br>y_res = 0.5<br>gid = 0<br>geom_coll = []<br>for xc in xrange(-18, 52, 1):<br>&nbsp;&nbsp;&nbsp; for yc in xrange(-36,36,1):<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; wkt = &quot;POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))&quot;% \<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; ( xc-x_res, yc-y_res, xc+x_res, yc-y_res, xc+x_res, yc+y_res,\<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; xc-x_res,yc+y_res,xc-x_res, yc-y_res)<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; geom = ogr.Geometry(type=ogr.wkbPolygon)<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; geom = ogr.CreateGeometryFromWkt(wkt)<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; geom.AssignSpatialReference ( t_srs )<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; geom_coll.append ( geom )<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; <br>for iw in xrange ( geom_world.GetGeometryCount() ):<br>&nbsp;&nbsp;&nbsp; gw = geom_world.GetGeometryRef ( iw)<br>&nbsp;&nbsp;&nbsp; for geom in geom_coll:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; if gw.Contains ( geom.Centroid()):<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; gid += 1<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; feat = ogr.Feature(feature_def=layer.GetLayerDefn())<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; feat.SetGeometryDirectly(geom)<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; feat.SetFID( 1 )<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; feat.SetField(&#39;id&#39;, gid)<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; layer.CreateFeature(feat)<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; feat.Destroy()<br>&nbsp;&nbsp;&nbsp; <br>ds.Destroy()<br><br clear="all"><br>-- <br>Department of Geography, University College London<br>Gower Street, London WC1E 6BT, UK<br>