在PostGIS中计算一个面要素表中的缝隙（Find gaps among polygons in PostGIS）

场景

在PostGIS中有一张面要素表，需要检查该表中的哪些地方有缝隙。
其中缝隙定义为这些多边形的并集中的环。
There is a surface feature table in PostGIS, and it is necessary to check which areas in the table have gaps.
The gaps are defined as the rings in the union of these polygons.
数据分布如下：

我们可以看出其中的缝隙，用红色框选一下：

第一步：计算全部要素的外包矩形

			WITH all_box AS (
			    SELECT ST_Envelope(ST_Union(geom)) AS box
			    FROM TARGET_TABLE
			)
			, unionB as (SELECT ST_Simplify(ST_Union(st_buffer(GEOM_COLUMN,TOLERANCE_VALUE)),TOLERANCE_VALUE) AS BGEOM FROM TARGET_TABLE) 

第二步：从外包矩形中依次去掉面要素

removed_box AS (
			    SELECT ST_Difference(a.box, B.BGEOM) AS x
			    FROM all_box a
			    JOIN unionB b ON TRUE
			)

第三步：去掉与边界接触的子要素

			, filtered_x AS (
			    SELECT ST_Difference(r.x, ST_Boundary(r.x)) AS y
			    FROM removed_box r
			)

第四步：如果是几何集合则分解为简单对象

			, otherpolygons as (SELECT (ST_Dump(f.y)).* FROM filtered_x f)
			, otherPolygons_sub as (select geom from otherpolygons where ST_NumInteriorRings(geom) > 0)
			, interior_rings AS (
			    SELECT ST_InteriorRingN(geom, generate_series(1, ST_NumInteriorRings(geom))) AS geom
			    FROM otherPolygons_sub t
			)
			select op.geom
			from otherpolygons op, all_box ab where 
			(st_touches(ST_ExteriorRing(ab.box),op.geom) = false)
			or
			(st_geometrytype(st_intersection(ST_ExteriorRing(ab.box),op.geom)) = ''ST_Point'')
			';

完整函数

create or replace function public.sp_check_gap(targetTb varchar,resultTb varchar,
	geomColumnName varchar default 'geom', curUserName varchar default '')
RETURNS character varying
LANGUAGE 'plpgsql'
AS $BODY$
	--DECLARE _RESULT VARCHAR;
	--DECLARE row record;
    declare _TOLERANCE FLOAT8 default 0.0000000000001;
   	declare _SQLTMP text;
    declare _ISPOINTTYPE BOOLEAN;
begin
	
	execute 'SELECT EXISTS (SELECT POSITION(''Polygon'' in ST_GEOMETRYTYPE(GEOM)) FROM ' ||targetTb|| ' LIMIT 1)' into _ISPOINTTYPE;
	-- 如果类型不是面，则不检查
	if _ISPOINTTYPE = false then 
		return -1;
	end if;

	--ID、要素ID、是否已修正
	EXECUTE 'CREATE TABLE IF NOT EXISTS ' || resultTb ||' (
		id varchar DEFAULT (((''ID_''::text || to_char(now(), ''YYYYMMDDHH24MISS''::text)) || ''_''::text) || "substring"(md5(random()::text), 1, 10)) NOT NULL,
		geom Geometry,
		corrected Boolean
	)';
	execute 'DELETE FROM ' || resultTb;

	_SQLTMP :='
			-- 第一步：计算全部要素的外包矩形
			WITH all_box AS (
			    SELECT ST_Envelope(ST_Union(geom)) AS box
			    FROM TARGET_TABLE
			)
			, unionB as (SELECT ST_Simplify(ST_Union(st_buffer(GEOM_COLUMN,TOLERANCE_VALUE)),TOLERANCE_VALUE) AS BGEOM FROM TARGET_TABLE) 
			-- 第二步：从外包矩形中依次去掉面要素
			, removed_box AS (
			    SELECT ST_Difference(a.box, B.BGEOM) AS x
			    FROM all_box a
			    JOIN unionB b ON TRUE
			)
			-- 第三步：去掉与边界接触的子要素
			, filtered_x AS (
			    SELECT ST_Difference(r.x, ST_Boundary(r.x)) AS y
			    FROM removed_box r
			)
			-- 第四步：如果是几何集合则分解为简单对象
			, otherpolygons as (SELECT (ST_Dump(f.y)).* FROM filtered_x f)
			, otherPolygons_sub as (select geom from otherpolygons where ST_NumInteriorRings(geom) > 0)
			, interior_rings AS (
			    SELECT ST_InteriorRingN(geom, generate_series(1, ST_NumInteriorRings(geom))) AS geom
			    FROM otherPolygons_sub t
			)
			insert into RESULT_TABLE (geom)
			select op.geom
			from otherpolygons op, all_box ab where 
			(st_touches(ST_ExteriorRing(ab.box),op.geom) = false)
			or
			(st_geometrytype(st_intersection(ST_ExteriorRing(ab.box),op.geom)) = ''ST_Point'')
			';

	_SQLTMP:= REPLACE(_SQLTMP,'RESULT_TABLE',resultTb);
	_SQLTMP:= REPLACE(_SQLTMP,'TARGET_TABLE',targetTb);
	_SQLTMP:= REPLACE(_SQLTMP,'GEOM_COLUMN',geomColumnName);
	_SQLTMP:= REPLACE(_SQLTMP,'TOLERANCE_VALUE',_TOLERANCE::VARCHAR);
	raise notice '_SQLTMP : %',_SQLTMP;
	execute _SQLTMP;
		
	RETURN resultTb;
end;
$BODY$;

测试调用

select public.sp_check_gap('TARGET_TABLE','RESULT_TABLE');

查询结果

select * from YOUR_RESULT_TABLENAME;

分析结果