ST_intersection devuelve polígonos erróneos

Question

Estoy intentando realizar un recorte por atributos en PostGIS entre dos shapefile que contienen 3000 polígonos. Ambos shapefiles están libres de errores (comprobados con el verificador de topología), ambos shapefile están usando EPSG:3857.

La siguiente consulta funciona bien para el 99,9% de los polígonos, pero para algunos crea polígonos incompletos o falla algunos cuando se cargan en QGIS. Los polígonos con errores tienen el mismo nombre utilizado como atributo de recorte. Realmente no entiendo por qué algunos polígonos no se generan en comparación con otros.

create table poly_20m_checked (gid serial PRIMARY KEY, polyID text, geometry geometry);
    insert into poly_20m_checked(polyID, geometry)
    select
    a.name as polyID,
    ST_Intersection(a.geom,b.geom) as geometry
    FROM donut20m a, symm20m b
    WHERE a.name=b.name
    Group by polyID, geometry

donut20m: symm20m: resultado (falta un polígono):

Archivos de ejemplo que contienen 4 polígonos, entre los cuales uno que devuelve un error (movido para hacer la localización anónima): https://drive.google.com/open?id=0B7mqg7Gq1NGVVVVsbTZqMlVhZnc

La razón por la que lo estoy haciendo en PostGIS es que no creo que haya una forma de recortar polígonos en 2 shapefiles usando atributos en QGIS. ¿Qué podría estar causando errores? Estoy muy abierto a sugerencias para evitar el problema. He pensado en dividir los polígonos por ID y luego recortar cada par de shapefiles pero no he conseguido codificarlo.

EDITAR: descargar para ver los polígonos de muestra. EDIT2: Intenté simplificar las geometrías en 1m, lo que hace que la cosa mejore un poco, pero todavía tiene errores

Answer 1

Desafortunadamente, las funciones ST_Intersection ocasionalmente dan una GeometryCollection donde no se esperaría. Este es el caso de tus datos. Cuando inspeccione la salida con ST_AsText verá que el segundo donut de la parte superior produce tanto un polígono como una cadena de líneas. Esta cadena de líneas es falsa, ya que sólo tiene 2 coordenadas idénticas, pero debido a algún proceso de redondeo, postgis (o probablemente GEOS) ve esto como una entidad separada. Esto ocurre especialmente con los polígonos que se supone que se superponen con precisión.

La solución es sencilla, utilizar siempre un ST_Dump en la salida de ST_Intersection. Ejemplo:

select
    a.name as polyID,
    (ST_Dump(ST_Intersection(a.geom,b.geom))).geom as geometry
    FROM tmp.sampledonut20m   a, tmp.samplesymm20m b
    WHERE a.name=b.name
    Group by polyID, geometry

Utiliza un filtro más adelante si quieres deshacerte de las rayas. Hay un buen ejemplo con una función de Dan Baston por aquí: PostGIS ST_Intersección de polígonos puede devolver líneas