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

[G. Allegri]

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 <giohappy at gmail.com>:
> 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
>