Here I am, back to office ;-) I've had support by IRC people. As I supposed there aren't obvious improvements to the routine code. The only, real, boost is given by using a spatial index. My polygonal layer is a shapefile, so I've created a quadtree index using the shptree mapserver utility, which produces a .qix index file that OGR can use to gain performance in spatial operations like mine. More informations (with the fresh add of the shptree option) can be found at the shapefile format page of gdal: http://www.gdal.org/ogr/drv_shapefile.html
thanks to christopher schmidt, hobu and the others on IRC. Have a good day! giovanni 2009/7/1 G. Allegri <[email protected]>: > Hi list. > I needed to make a simple routine to create a regular point grid (as > csv) on the base of an input polygons layer and its attributes. > I've compined the two needs: > > - verify point/polygon containment > - extract polygon attribute and attach it to the point feature > > so I've used the <layer>.SetSpatialFilter(<geometry>) method. When I > need to evaluate some thousends of points the performance gets really > low. > This is the simple algoithm (I've omitted various error checking here): > > > def doGrid(self,spacing=10,fout='pointGrid',attr=None): > > spacingm = spacing*1000 # POINTS SPACING IN METERS > point = ogr.Geometry(ogr.wkbPoint) > fout = "%s_%skm.txt" % (fout,spacing) > header = "lat\tlon\tcode" > rowtpl = "%s\t%s\t%s" > try: > f = open(fout,'w') > except: > print "Can't open %s" % s > exit(1) > > if attr: > header += "\tattr_value" > rowtpl += "\t%s" > f.write(header+"\n") > > plat = self._mlat > r = 1 > while (plat < self._Mlat): # start from the lower row > plon = self._mlon > r += 1 > c = 1 > while (plon < self._Mlon): # start from the leftmost > column > code = "r"+str(r)+"c"+str(c) > strpars = (plat,plon,code) > if attr: > point.AddPoint_2D(plon,plat) # I use > this to add the actual > point to the Point geometry > self._inlayer.ResetReading() # reset > the layer pointer > self._inlayer.SetSpatialFilter(point) > # do the filtering. It was > a try and it worked: passing a point geometry is ok, I don't need a > box... > feat = self._inlayer.GetNextFeature() > # I try to get a feature. > If it fails, the point is out of the layer. > if feat: > try: > val = > feat.GetField(attr) # I extract the attribute > strpars += (val,) > rowstring = (rowtpl % > strpars)+"\n" > except ValueError: > print "attribute > doesn't belong to feature" > else: > rowstring = "" > else: > rowstring = (rowtpl % strpars)+"\n" > > f.write(rowstring) > plon += spacingm > c += 1 > plat += spacingm > f.flush() > f.close() > print "Punti totali: %s" % ((c-1)*(r-1)) > > > thanks for any suggestion, > giovanni > _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
