PostGIS ST_Intersección de polígonos puede devolver líneas

Question

Cuando se recortan los polígonos de una tabla con polígonos de otra, ST_Intersection puede devolver un conjunto de resultados que pueden ser manejados con ST_Dump. Las geometrías múltiples devueltas no son necesariamente ST_Polygon sino también ST_LineString (probablemente también un punto). Así que cuando se ejecuta una consulta

INSERT INTO c (geom)
  (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
  FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom));

al intentar rellenar la tabla "c" con polígonos recortados, falla con ERROR: El tipo de geometría (LineString) no coincide con el tipo de columna (Polygon)

Hice otra sentencia SELECT anidada para que sólo salieran las geometrías de los polígonos, como:

INSERT INTO c (geom)
  (SELECT geom FROM
    (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
    FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom))) AS cl
    WHERE ST_GeometryType(cl.geom)='ST_Polygon');

Como esto es un poco engorroso, me pregunto si hay una solución más elegante para eliminar las geometrías no válidas.

Answer 1

Este podría ser un buen punto para utilizar una función del lenguaje SQL. Aquí hay una rápida que debería funcionar para esta situación:

CREATE OR REPLACE FUNCTION PolygonalIntersection(a geometry, b geometry)
RETURNS geometry AS $$
SELECT ST_Collect(geom)
FROM 
(SELECT (ST_Dump(ST_Intersection(a, b))).geom 
UNION ALL
-- union in an empty polygon so we get an 
-- empty geometry instead of NULL if there
-- is are no polygons in the intersection
SELECT ST_GeomFromText('POLYGON EMPTY')) SQ
WHERE ST_GeometryType(geom) = 'ST_Polygon';
$$ LANGUAGE SQL;

Esto retendrá los componentes poligonales de una intersección, pero desechará todo lo demás. Siempre devuelve un MultiPolígono, incluso si tiene uno o ningún componente.

WITH 
      square   as (SELECT ST_GeomFromText('POLYGON ((0 0, 0  1,  1  1,  1  0, 0 0))') AS geom),
biggersquare   as (SELECT ST_GeomFromText('POLYGON ((0 0, 0 10, 10 10, 10  0, 0 0))') AS geom),
adjacentsquare as (SELECT ST_GeomFromText('POLYGON ((0 0, 1  0,  1 -1, -1 -1, 0 0))') AS geom)   

SELECT ST_AsText(PolygonalIntersection(square.geom, biggersquare.geom))
  FROM square, biggersquare;
--"MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)))"

SELECT ST_AsText(PolygonalIntersection(square.geom, adjacentsquare.geom))
  FROM square, adjacentsquare;
--"MULTIPOLYGON(EMPTY)"

Answer 2

Muy buena respuesta de @dbaston. Sin embargo, devolver una geometría vacía en lugar de null podría causar problemas porque la geometría vacía devuelta no tiene srid. St_Intersection también podría devolver MultiPolygon. Esta función actualizada fue realmente útil para mí:

CREATE OR REPLACE FUNCTION PolygonalIntersection(a geometry, b geometry)
RETURNS geometry AS $$
SELECT ST_Collect(geom)
FROM 
(SELECT (ST_Dump(ST_Intersection(a, b))).geom 
) SQ
WHERE ST_GeometryType(geom) = 'ST_Polygon' OR ST_GeometryType(geom) = 'ST_MultiPolygon';
$$ LANGUAGE SQL;