[gdal-dev] best practice to evaluate (fast) point in polygon and attribute overlay with python and ogr?

[G. Allegri]

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