使用PostGIS对数据做拓扑抽稀

抽稀，是指减少图形的节点密度。目前在做大数据展示时，基本上都需行政区划数据将shp转geojson格式在web端进行加载渲染，如果图形节点数能达到几万个，这样的图形在渲染和空间关系判断等环节中效率会很慢，所以在一些要求不是很严格的情况下，可以使用抽稀方法减少图形的节点密度，加快处理效率。

对于抽稀的要求，最重要的是要基本保持原图形的形状，这个大部分软件arcgis、qgis等抽稀算法都可以达到。但对于质量比较高的情况下，抽稀后相邻图形不能出现重叠和缝隙，利用上述软件很难实现。

我们将数据导入到PostGIS数据库中

导入结果如下：

1. 使用ST_Simplify

PostGIS本身有一个函数ST_Simplify方法对要素进行抽稀的结果，可以看到比较多的重叠和缝隙

ST_Simplify -使用Douglas-Peucker算法返回给定几何图形的“简化”版本。

1.1. 概要

geometry ST_Simplify(geometry geomA, float tolerance, boolean preserveCollapsed);
1

1.2. 描述

使用Douglas-Peucker算法返回给定几何图形的“简化”版本。 将实际上只做一些(多)线和(多)多边形，但你可以安全地调用它与任何类型的几何。 由于简化是在逐个对象的基础上进行的，所以您也可以向该函数提供一个GeometryCollection。

“保存坍塌”标志将保留那些在其他情况下由于容忍度太小的对象。 例如，1米长的线简化为10米的公差。 如果给出了保护旗，这条线将不会消失。 这个标志对渲染引擎很有用，可以避免大量非常小的物体从地图上消失，留下令人惊讶的空白

注意，拓扑图可能不会被保留，并可能导致无效的几何图形。 使用(参见ST_SimplifyPreserveTopology)保存拓扑。

1.3. 具体实现

update bj_county_simplify set  the_geom =st_simplify(the_geom, 20);

update bj_county_simplify set pointcount=  ST_NPoints(the_geom);
1
2
3

• 节点个数锐减，但图斑没有出现大范围的变化

但是

• 图形共边处瑕疵严重

2. 使用ST_SimplifyDouglas-Peucker算法简化拓扑边

2.1. 准备环境

CREATE EXTENSION postgis;

CREATE EXTENSION postgis_topology;
1
2
3

2.2. 数据备份

create table bj_county_topo_simplify_20 as select * from bj_county_ys
1

2.1. 准备SimplifyEdgeGeom函数

https://gist.github.com/leplatrem/5729022 使用ST_SimplifyDouglas-Peucker算法简化拓扑边。

CREATE OR REPLACE FUNCTION SimplifyEdgeGeom(atopo varchar, anedge int, maxtolerance float8)
RETURNS float8 AS $$
DECLARE
  tol float8;
  sql varchar;
BEGIN
  tol := maxtolerance;
  LOOP
    sql := 'SELECT topology.ST_ChangeEdgeGeom(' || quote_literal(atopo) || ', ' || anedge
      || ', ST_Simplify(geom, ' || tol || ')) FROM '
      || quote_ident(atopo) || '.edge WHERE edge_id = ' || anedge;
    BEGIN
      RAISE DEBUG 'Running %', sql;
      EXECUTE sql;
      RETURN tol;
    EXCEPTION
     WHEN OTHERS THEN
      RAISE WARNING 'Simplification of edge % with tolerance % failed: %', anedge, tol, SQLERRM;
      tol := round( (tol/2.0) * 1e8 ) / 1e8; -- round to get to zero quicker
      IF tol = 0 THEN RAISE EXCEPTION '%', SQLERRM; END IF;
    END;
  END LOOP;
END
$$ LANGUAGE 'plpgsql' STABLE STRICT;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

2.3. 拓扑处理

-- 创建拓扑
SELECT topology.CreateTopology('xzq_test_topo', 4527);

-- 创建拓扑字段
SELECT topology.AddTopoGeometryColumn('xzq_test_topo', 'public', 'bj_county_topo_simplify_20', 'f_topo_geom', 'MULTIPOLYGON');

-- 将拓扑信息写入到新增的拓扑字段中
UPDATE bj_county_topo_simplify_20 SET f_topo_geom = topology.toTopoGeom(the_geom, 'xzq_test_topo', 1, 0.001);

-- 简化拓扑中的边信息
SELECT SimplifyEdgeGeom('xzq_test_topo', edge_id,20) FROM xzq_test_topo.edge;

-- 简化后的变重新回写到图形中
UPDATE bj_county_topo_simplify_20 SET the_geom = f_topo_geom::geometry;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

2.4. 删除拓扑

SELECT topology.DropTopology('xzq_test_topo');
1

3. 使用ST_SimplifyvwVisvalingam-Whyatt算法简化拓扑边

2.1. 准备环境

CREATE EXTENSION postgis;

CREATE EXTENSION postgis_topology;
1
2
3

2.2. 数据备份

create table bj_county_topo_simplifyvw_20 as select * from bj_county_ys;
1

2.1. 准备SimplifyvwEdgeGeom函数

CREATE OR REPLACE FUNCTION SimplifyvwEdgeGeom(atopo varchar, anedge int, maxtolerance float8)
RETURNS float8 AS $$
DECLARE
  tol float8;
  sql varchar;
BEGIN
  tol := maxtolerance;
  LOOP
    sql := 'SELECT topology.ST_ChangeEdgeGeom(' || quote_literal(atopo) || ', ' || anedge
      || ', ST_Simplifyvw(geom, ' || tol || ')) FROM '
      || quote_ident(atopo) || '.edge WHERE edge_id = ' || anedge;
    BEGIN
      RAISE DEBUG 'Running %', sql;
      EXECUTE sql;
      RETURN tol;
    EXCEPTION
     WHEN OTHERS THEN
      RAISE WARNING 'Simplification of edge % with tolerance % failed: %', anedge, tol, SQLERRM;
      tol := round( (tol/2.0) * 1e8 ) / 1e8; -- round to get to zero quicker
      IF tol = 0 THEN RAISE EXCEPTION '%', SQLERRM; END IF;
    END;
  END LOOP;
END
$$ LANGUAGE 'plpgsql' STABLE STRICT;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

2.3. 拓扑处理

-- 创建拓扑
SELECT topology.CreateTopology('xzq_testvw_topo', 4527);

-- 创建拓扑字段
SELECT topology.AddTopoGeometryColumn('xzq_testvw_topo', 'public', 'bj_county_topo_simplifyvw_20', 'f_topo_geom', 'MULTIPOLYGON');

-- 将拓扑信息写入到新增的拓扑字段中
UPDATE bj_county_topo_simplifyvw_20 SET f_topo_geom = topology.toTopoGeom(the_geom, 'xzq_testvw_topo', 1, 0.001);

-- 简化拓扑中的边信息
SELECT SimplifyvwEdgeGeom('xzq_testvw_topo', edge_id,200) FROM xzq_testvw_topo.edge;

-- 简化后的变重新回写到图形中

UPDATE bj_county_topo_simplifyvw_20 SET the_geom = f_topo_geom::geometry;

update bj_county_topo_simplifyvw_20 set pointcount=  ST_NPoints(the_geom);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

SELECT topology.DropTopology('xzq_testvw_topo');
1

4. postgis中的拓扑结构

4.1. edge_data边数据信息

4.2. face面信息

4.3. node边数据信息

4.4. relation关系信息