Mejorar el rendimiento de la diferencia entre dos capas con PostGIS

Hace poco uno de nuestros clientes nos solicitó ayuda porqué un geoprocesamiento que estaban desarrollando en PostGIS tardaba demasiado en finalizar.

No podemos compartir los datos ni el código real, pero si unas queries simplificadas que valgan como ejemplo genérico. Lo haremos paso a paso y mediante CREATE TABLE se puede visualizar el proceso. El algoritmo final es algo distinto, pero hemos tratado de mantener el ejemplo simple y la idea base es la misma.

Identificamos el problema central en una query que intentaba generar la diferencia entre dos capas de la base de datos del siguiente modo:

WITH
    l1 AS (
        SELECT ST_Union(geom) AS geom FROM boundary_l1
    ),  
    l2 AS (
        SELECT ST_Union(geom) AS geom FROM boundary_l2
    )                  
SELECT ST_Difference(l1.geom, l2.geom) FROM l1, l2;

Es decir agrupamos en una única geometría cada una de las capas y luego las «restamos». La aproximación es simple y clara, pero operaciones como ST_Difference o ST_Intersects son computacionalmente costosas. En general serán más rápidas al trabajar sobre polígonos simples sobre todo si conseguimos filtrar adecuadamente sobre que pares de geometrías realizar la operación.

Conocer los datos en estos casos es útil, ambas capas son multipolígonos de regiones administrativas, con l2 siendo de menor nivel (pongamos municipios) que l1 (pongamos provincias). En este caso se podían agrupar las geometrías de cada capa por «islas» independientes entre sí sobre las que realizar la diferencia.

Así que lo primero es generar polígonos simples para cada capa, que luego agruparemos disolviendo las fronteras y volvemos a romper a continuación en polígonos simples.

-- l1_simple already contains simple dissolved polygon
-- cast to geometry(polygon, 4326) simplifies visualize results in GIS, and
-- checks the types of the generated geometries
CREATE TABLE l1_simple AS
    SELECT (ST_Dump(geom)).geom::geometry(polygon, 4326) geom FROM boundary_l1;

CREATE TABLE l2 AS
    SELECT (ST_Dump(geom)).geom::geometry(polygon, 4326) geom FROM boundary_l2;
CREATE TABLE l2_union AS
    SELECT ST_Union(geom)::geometry(multipolygon, 4326) geom FROM l2;
CREATE TABLE l2_simple AS
    SELECT (ST_Dump(geom)).geom::geometry(polygon, 4326) geom FROM l2_union;

Examining the results shows that l2_simple has many more geometries tha

Al examinar los resultados se observa como l2_simple tiene muchas más geometrías de las esperadas. Podría ser por geometrías mal digitalizadas que nos se tocaran, pero

SELECT ST_Area(geom) FROM l2_simple ORDER BY 1;

Nos permite enseguida detectar una enorme cantidad de slivers que pueden ser eliminados.

DELETE FROM l2_simple WHERE ST_Area(geom) < 0.0001;

Ahora que tenemos polígonos simples disueltos para cada zona de interés sólo queda hacer la diferencia entre cada par de geometrías que se toquen. Pero primero nos quedaremos con las geometrías de l1 para las que no haya equivalente en l2.

CREATE TABLE RESULT AS
    SELECT ST_Multi(l1.geom)::geometry(multipolygon, 4326) AS geom
    FROM l1_simple l1
    WHERE l1.ctid NOT IN (
        SELECT l1.ctid
        FROM a JOIN l2_simple l2
        ON ST_Intersects (l1.geom, l2.geom)
    );
INSERT INTO RESULT
    SELECT ST_Difference(l1_simple.geom, l2_simple.geom)
    FROM l1_simple JOIN l2_simple
    ON ST_Intersects(l1_simple.geom, l2_simple.geom);

La query inicial la paramos en la primera prueba cuando llevaba más de 20 minutos ejecutándose. El algoritmo modificado se ejecuta en menos de un minuto.

Recuerda, simple suele ser más rápido que multi. Y si necesitas ayuda con tus geoprocesos de PostGIS escríbenos.