在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;

在这里插入图片描述

分析结果

在这里插入图片描述

  • 9
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

丷丩

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值