本文主要是介绍在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 boxFROM 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 xFROM all_box aJOIN unionB b ON TRUE)
第三步:去掉与边界接触的子要素
, filtered_x AS (SELECT ST_Difference(r.x, ST_Boundary(r.x)) AS yFROM 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 geomFROM otherPolygons_sub t)select op.geomfrom 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;
beginexecute '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 boxFROM 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 xFROM all_box aJOIN unionB b ON TRUE)-- 第三步:去掉与边界接触的子要素, filtered_x AS (SELECT ST_Difference(r.x, ST_Boundary(r.x)) AS yFROM 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 geomFROM otherPolygons_sub t)insert into RESULT_TABLE (geom)select op.geomfrom 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;
分析结果
这篇关于在PostGIS中一个面要素表中的缝隙(Find gaps among polygons in PostGIS)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!