Obtención del complemento geoespacial de un conjunto de polígonos a un cuadro delimitador en PostGIS

Question

Quiero obtener la diferencia entre un cuadro delimitador (un solo rectángulo) y una colección de polígonos (un montón de geometrías complejas) dentro de ese cuadro delimitador.

Es decir, digamos que mi cuadro delimitador abarca los Estados Unidos, y mis polígonos son los polígonos de los centros urbanos, quiero obtener la geometría correspondiente a las regiones dentro de mi cuadro delimitador que NO son un centro urbano.

Puedo obtener fácilmente los centros urbanos dentro de mi cuadro delimitador:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT bb_uc.geometry
FROM bb_uc;

Esta consulta se ejecuta en unos 10 segundos en mi caso, y devuelve ~100k filas con geometrías de polígonos.

Ahora, cuando intento obtener el "complemento" de este conjunto de datos espaciales con relación a mi cuadro delimitador:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT ST_DIFFERENCE(ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326), ST_UNION(bb_uc.geometry)) as geometry
FROM bb_uc

Esta consulta tarda una eternidad en ejecutarse (he parado manualmente después de 3h). No había imaginado que esta operación fuera tan compleja (básicamente, aplicar una máscara a mi cuadro delimitador)

Como nota aparte, nótese que estoy utilizando una sentencia WITH "innecesaria" porque a veces añado una cláusula WHERE en mis resultados.

¿Qué puedo hacer para aumentar drásticamente mi rendimiento para llegar a la geometría "complementaria" de mi intersección? La página de manual de ST_Difference menciona no usar geometryCollection, de ahí que use ST_Union, que es la operación que tarda una eternidad. Sin embargo, no puedo pensar en una forma de evitar esto.

postgis_version: 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

Answer 1

SO,

Si alguien tiene problemas para resolver esta cuestión, proceda como sigue:

mis geodatos de origen son un polígono llamado poly_extent y polígonos administrativos sin agujeros llamados adm_polygons ver Figura 1.

Ejecuta el script

WITH 
    tbla AS (SELECT (ST_DumpPoints(geom)).geom FROM adm_polygons),
    tblb AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tbla),
    tblc AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tblb a JOIN poly_extent b ON ST_Intersects(a.geom, b.geom)),
    tbld AS (SELECT ST_Difference(a.geom, b.geom) geom FROM tblc a JOIN adm_polygons b ON ST_Intersects(a.geom, b.geom) AND ST_Overlaps(a.geom, b.geom))
    SELECT ST_Union(geom) geom FROM tbld;

Véase el resultado en la figura 2.

Buena suerte a todos :-),

Soluciones originales ...

Este script se llama - ST_CarvesPolygons...