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
