[Gdal-dev] ogr find neighbours

[Howard Butler]

Didrik,

In my experience, although I don't know for sure 
if really true, testing for Overlaps/Disjoint is 
much faster than doing Touches, which probably 
has to walk a lot of vertices and may take on 
some precision model stuff to do comparisons.  If 
the two geometries are Disjoint, or not Overlaps, 
then they can't Touch.  Only run the Touch test 
on those pairs that really have a chance to do so.

I could be totally wrong here though, and GEOS 
may already be doing this optimization for us (at 
least testing the bounding boxes to throw out 
obvious cases), but I haven't looked at the 
source code to know  for sure.

Another thing I would add is if you are trying to 
do something like a dissolve, you need to 
bucketize the neighbors into groups so that you 
can go along and union them together.  I have 
some code that uses the Python Cartographic 
Library to do a dissolve.  Maybe some of the 
concepts of what you are doing are similar.  Let 
me know if you are interested.

Howard

At 3:29 PM +0200 5/10/06, Didrik Pinte wrote:
>Content-Type: multipart/signed; micalg=pgp-sha1;
>	protocol="application/pgp-signature";
>	boundary="=-ylVKZCaZfiSwx0EHy8kT"
>
>Hi,
>
>In the process of implementing a model, I need to get a list of
>neighbours for a shapefile containing polygons.
>
>The file contains 11 polygons having lots of points.
>
>The very simple python code i've developped now is just to test if for
>each polygon it touches or overlaps the other polygons from the file.
>
>BUT :
>
>with 11 polygons, I makes 122 (11 * 11) tests of geometry.Touches and
>122 tests of geometry.Overlaps. This is done in 60 seconds and thus is
>too slow.
>
>By hand, I can do it in 15 seconds maximum taking into account the time
>to open the file with GRASS ;-)
>
>Does anyone here have a great idea on how to do this more quickly ?
>
>I've pasted the code below.
>
>Thanks for your help
>
>Didrik
>
>--
>
>def get_neighbours(self,queryLayer):
>  '''
>  Returns a list of references to adjacent cells
>  '''
>  self.neighbours = []
>  cobj = self.geom
>  queryLayer.ResetReading()
>  queryFeat = queryLayer.GetNextFeature()
>  while queryFeat is not None:
>   geometry = queryFeat.GetGeometryRef()
>   found = buffer.Touches(geometry)
>   if found > 0:  
>     self.neighbours.append(queryFeat.Clone())
>   else:
>    found = buffer.Overlaps(geometry)
>    if found > 0:
>     self.neighbours.append(queryFeat.Clone())
>    queryFeat.Destroy()
>    queryFeat = queryLayer.GetNextFeature()
>   queryLayer = None
>   return self.neighbours
>
>
>Content-Type: application/pgp-signature; name=signature.asc
>Content-Description: Ceci est une partie de message
>	numériquement signée
>
>Attachment converted: Macintosh HD:signature 331.asc (    /    ) (003EDC23)
>_______________________________________________
>Gdal-dev mailing list
>Gdal-dev at lists.maptools.org
>http://lists.maptools.org/mailman/listinfo/gdal-dev