GEOS Python: Tomando la unión de varias geometrías

Question

Estoy usando la biblioteca GEOS/OGR en Python. Tengo archivos llenos de geometrías que se superponen, y quiero convertirlos en archivos de geometrías que no se superponen.

La forma en que he estado abordando esto es simplemente leer todas las geometrías de un archivo, y luego hacer un bucle sobre cada geometría y tomar la Unión() de esa con todas las geometrías anteriores. Así que algo como:

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

La mayoría de las veces parece que esto funciona, pero ocasionalmente Union() no devuelve ninguna. Mi depuración revela que en estos casos:

• IsValid() devuelve verdadero tanto para g como para current_union
• Intersects() return true between g and current_union
• GetGeometryCount() para current_union típicamente devuelve un número sobre 100.
• Sospecho que el error puede ser algún tipo de problema de memoria porque parece que el error "se mueve" dependiendo del contexto. Es decir, si digo que lo ejecute en los archivos 1-100, podría fallar en el archivo 20, pero si digo que lo ejecute en los archivos 20-100, podría fallar en el archivo 40. Nunca uso Destroy() en ninguna geometría. Sin embargo, el Administrador de Tareas no muestra que me esté quedando sin memoria.

De todos modos, tal vez haya una forma más fácil de hacer esto. ¿O alguien tiene alguna idea de por qué se está estrellando aquí?

¡Muchas gracias de antemano!

Answer 1

No hay necesidad de clonar ninguna geometría. También, iterar sobre las geometrías[1:] para salvar una Unión().

Si osgeo.ogr tiene algo como Shapely's cascaded_union usar eso en lugar de acumular uniones en un bucle.

Answer 2

(Esta es una vieja pregunta, pero sigue siendo el principal resultado del GIS.SE para una búsqueda de "Python Union", así que estoy ofreciendo una respuesta. Espero que sea de alguna utilidad).

Python osgeo incluye un UnionCascaded método que implementa el GEOS subyacente CascadedUnion método. Es un poco extraño que mientras Union toma dos geometrías (en forma geom1.Union(geom2) ), UnionCascaded opera en una colección de geometría. Según la documentación de GEOS, este método es más eficiente que realizar una unión de piezas en geometrías individuales.

No sé exactamente lo que estás pasando a la función en geometries . Si es una colección de geometría, posiblemente geometries.UnionCascaded() funcionará. Los errores discutidos en los comentarios (por el OP Ben y los sgillies usuarios) podrían haber sido arreglados en los últimos tres años. Union trabajará para combinar geometrías de diferentes dimensiones (por ejemplo, polígonos y líneas), y posiblemente UnionCascaded también, pero me cuesta entender cómo sería útil o qué harías con el objeto resultante. Así que, como ejemplo, voy a asumir que geometries es una capa con características de la misma dimensión. Creo que la forma más fácil de operar lo que quieres hacer es recoger tus geometrías en un multi y luego llamar a UnionCascaded en ese objeto. Por ahora, asumamos que estamos trabajando con polígonos.

def dedup(geometries):
    """Return a geometry that is the union of all geometries."""
    if not geometries:  return None
    multi  = ogr.Geometry(ogr.wkbMultiPolygon)
    for g in geometries:
        multi.AddGeometry(g.geometry())
    return multi.UnionCascaded()

Posiblemente (no lo he probado) si reemplazas ogr.wkbMultiPolygon con ogr.wkbUnknown esto funcionaría para las líneas o puntos. Pero posiblemente el tipo de geometría de la salida también sería desconocido, y eso puede hacer difícil trabajar con él.