Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am using the GEOS/OGR library in Python. I have files full of geometries that overlap, and I want to convert these to files of geometries that don't overlap.

The way I have been approaching this is to just read in all of the geometries from a file, and then looping over each geometry and taking the Union() of that with all of the previous geometries. So something like:

def dedup(geometries):
   """Return a geometry that is the union of all geometries."""
   if not geometries:  return None
   current_union = geometries[0].Clone()
   for g in geometries:
      current_union = current_union.Union(g).Clone()
   return current_union

Most of the time this appears to work, but occasionally Union() will return None. My debugging reveals that under these cases:

  • IsValid() returns true for both g and current_union
  • Intersects() return true between g and current_union
  • GetGeometryCount() for current_union typically returns a number over 100.
  • I suspect the error might be some kind of memory issue because it seems like the error "moves" depending on context. Ie if I say run this on files 1-100, it might crash on file 20, but if I say run it on files 20-100, it might crash on file 40. I never use Destroy() on any geometries. Task Manager does not show me running out of memory though.

Anyway, maybe there is an easier way to do this? Or else does anyone have any idea why it is crashing here?

Thanks much in advance!

share|improve this question
add comment

1 Answer

No need to clone any geometries. Also, iterate over geometries[1:] to save one Union().

If osgeo.ogr has something like Shapely's cascaded_union use that instead of accumulating unions in a loop.

share|improve this answer
    
There is a UnionCascaded function apparently but I don't think it works. There's no documentation for it in the SWIG wrapper and it throws an exception if you pass any arguments into it. I am starting to think there might be a bug in the GEOS Union() function as well. –  Ben Jun 3 '12 at 23:01
    
I've seen reports of union bugs in the GEOS tracker. Find what your version is and look in trac.osgeo.org/geos/… to see whether any issue has been closed since or if it looks like one of the open ones. –  sgillies Jun 4 '12 at 14:38
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.