[Gdal-dev] ogr.py Intersect with wkbGeometryCollection

Thu May 5 18:37:33 EDT 2005

> -----Original Message-----
> From: Frank Warmerdam [mailto:fwarmerdam at gmail.com] 
> Sent: Thursday, May 05, 2005 11:22 AM
> To: Ryan, Adam
> Cc: gdal-dev at remotesensing.org
> Subject: Re: [Gdal-dev] ogr.py Intersect with wkbGeometryCollection
> 
> 
> On 5/5/05, Ryan, Adam <ARyan at co.linn.or.us> wrote:
> > 
> > Hi,
> > 
> > I tried to use Geometry.Intersect with a wkbGeometryCollection type 
> > Geometry and got hung.  Is this not permited?
> 
> Adam,
> 
> I am not aware of a problem with this, but then I doubt I 
> have specifically tried it.  Can you prepare a relatively 
> simple example to trigger this 
> problem?

Frank, thanks for your reply on this.  I'm still on the steep learning curve
here so bear with me.  I've made a little test script and have noted what
works and what doesn't.

The goal for me is to write a generic theme-on-theme spatial query method,
where you would pass an ogrQueryLayer and an ogrSelectLayer with a
selectList (of record indexes obtained from a previous query that likely
involved mapscript and not ogr directly) and a method (for this test I just
do an intersect) and a mode (new, and, or, xor, not - for later) and
optionally an queryList(again, later).

This test is a step towards that goal.  Any advice on how to optimize for
speed would be greatly appreciated (including rewriting using a whole
different approach).

But first, one comment about pymod.ogr.Geometry.Empty().
When I try:
selectPoly.Empty()
I get:
TypeError: Empty() takes exactly 2 arguments (1 given)
If I change ogr.py from:
    def Empty( self, other_geom ):
        return _gdal.OGR_G_Empty( self._o, other_geom._o )
to:
    def Empty( self ): #   , other_geom ):
        return _gdal.OGR_G_Empty( self._o ) #   , other_geom._o )
It seems to work.

And my test script.  Note that I'm using Howard Butler's ogr.py for python24
on Windows XP:
--------------

#!C:\python24\python.exe

from pymod import ogr,gdal

#   queryLayer is the layer I want a result list for
#   selectLayer is the layer that's previously been queried
#   selectList is the list of record numbers from the previous query
#
#   This will return a list of queryLayer indexes, or an empty list  
#
def intersectByList(queryLayer,selectLayer,selectList):

    from pymod import ogr,gdal
    import sys
    import time

    #   Keep track of time so I can try to optimize performance
    starttime = time.clock()

    #   Start a polygon geometry of feature envelopes
    #   that get all unioned together
    selectEnvelopes = ogr.Geometry(ogr.wkbPolygon)

    #   Start a geometryCollection of select Geometries
    selectGeomCol = ogr.Geometry(ogr.wkbGeometryCollection)

    #   Create a couple of geometries that I use when processing
    #   the envelopes
    selectPoly = ogr.Geometry(ogr.wkbPolygon)
    selectRing = ogr.Geometry(ogr.wkbLinearRing)

    #   Step through the passed record list
    for selectRec in selectList:

        #   Get the select Feature and Geometry
        selectFea = selectLayer.GetFeature(selectRec)
        selectGeom = selectFea.GetGeometryRef()

        #   Add the Geometry to our collection
        selectGeomCol.AddGeometry(selectGeom.Clone())

        #   Add the points of the envelope to the selectRing
        #   and add it to the selectPoly 
        selectEnv = selectGeom.GetEnvelope()
        selectRing.AddPoint(selectEnv[0],selectEnv[2])
        selectRing.AddPoint(selectEnv[1],selectEnv[2])
        selectRing.AddPoint(selectEnv[1],selectEnv[3])
        selectRing.AddPoint(selectEnv[0],selectEnv[3])
        selectRing.AddPoint(selectEnv[0],selectEnv[2])
        selectPoly.AddGeometry(selectRing)

        #   Union the existing selectEnvelopes with this polygon
        selectEnvelopes = selectEnvelopes.Union(selectPoly)

        #   CleanUp
        selectPoly.Empty()
        selectRing.Empty()
        selectFea.Destroy()

    #   Check time
    print "Made envelops in: ", time.clock() 

    #   Set a spatial filter on the qlayer with selectEnvelopes
    queryLayer.SetSpatialFilter(selectEnvelopes) 

    #   Check time
    print "SetSpatialFilter in : ", time.clock() 

    #   Start the return list
    queryList = []

    #    Step through the query features and test for intersect
    queryFeature = queryLayer.GetNextFeature()    
    while queryFeature is not None:

        queryGeom = queryFeature.GetGeometryRef()

        #-------------------------------------
        """#***** This doesn't work:

        #   Test for Intersect against select collection
        if queryGeom.Intersect(selectGeomCol):
            queryList.append(queryFeature.GetFID())

        """
        #-------------------------------------
        #   So I do the following ...                       

        #   Step through the list of features in the select collection
        for selectI in range(selectGeomCol.GetGeometryCount()):

            #   Get one select geometry out of the collection
            selectGeom = selectGeomCol.GetGeometryRef(selectI)

            #   Test for Intersect and if true move on to next query Feature
            if queryGeom.Intersect(selectGeom):
                queryList.append(queryFeature.GetFID())
                break
        #-------------------------------------

        queryFeature.Destroy()
        queryFeature = queryLayer.GetNextFeature() 

    #   Check time
    print "Intersected in : ", time.clock() 

    #   Return the list
    return queryList

#===========================    

#   RUN TIME    

#   Open ogr select layer
ogrSDS = ogr.Open("c:/gis/shapefiles/surveys.shp")
ogrSlayer = ogrSDS.GetLayer(0)

slayerSet = [27735, 27776, 27782, 27784, 27794, 27795, 27806, 27810, 27812,
27815, 27816, 27829, 27833, 27835, 27838, 27846, 27850, 27851, 27863, 27870,
27877, 27878, 27883, 27884, 27892, 27898, 27899, 27902, 27903, 27905, 27906,
27908, 27910, 27917, 27925, 27942, 27951, 27952, 27953] 

#   Open ogr query layer
ogrQDS = ogr.Open("c:/gis/shapefiles/taxlots.shp")
ogrQlayer = ogrQDS.GetLayer(0)

#   Get the intersect list
qlayerSet = intersectByList(ogrQlayer,ogrSlayer,slayerSet)

print qlayerSet

ogrSDS.Destroy()
ogrQDS.Destroy() 

-----

Thanks for any help.  The speed is OK if I'm noodling around on my own, but
kind of slow for a weblication.  Like I said, I'm still learning lots daily.
Cheers,

Adam