基于PostGIS拓扑进行shape面要素抽稀

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

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

图1原始行政区划shape数据

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

图2 抽稀后要素之间有空隙和叠加

2、接下来咱们通过建立要素严格拓扑关系的情况下，进行抽稀方法可达到这种效果。（原理：首先将shape行政区划建立拓扑，找到相邻要素的公共边；其次在不改变拓扑关系下对公共边进行抽稀；最后将抽稀后的公共边反写到shape图层几何字段上。）

(1)创建POSTGRESQL数据库

可在已有的数据中建立一个新的架构模式，如果数据已经入库，则可以跳过该操作。

（2）创建POSTGIS和TOPOLOGY扩展

1 2	CREATE EXTENSION postgis; CREATE EXTENSION postgis_topology;

（3）导入数据

如果处理的PostgreSQL数据库中已有数据，则跳过该操作。

1）这里使用GDAL提供的ogr2ogr命令行工具来做数据导入Shapefile格式的数据：

1	ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5432 dbname=test user=postgres password=1234" /test/原始数据/xzqh.shp -lco GEOMETRY_NAME=geom -lco FID=pk_uid -nlt PROMOTE_TO_MULTI -progress

2）使用postgis自带的shape数据导入工具

3）利用Qgis的数据库管理器导入

（4）复制数据

由于处理过程中将会给数据添加字段，为了不破坏数据，将会复制一份数据进行拓扑抽稀操作。

1	create table test_bak as select * from test --test为shape数据

（5）创建拓扑（只是添加拓扑模式，并未真实空间数据的拓扑）

创建拓扑需要传入拓扑模式名称和数据的SRID，拓扑模式名称使用的是数据名称后面加上“_topo”标识，数据的SRID通过geometry_columns视图看到是4326。

1	SELECT topology.CreateTopology('test_topo', 4490); --4490为数据的srid，当然也可以用Find_SRID()函数

创建拓扑之后，在数据库中生成一个新的模式，用来存储拓扑信息：

geometry_columns中添加了相应的拓扑表信息：

 理论拓扑关系

 topo模式表存储的关系

1）dege_data

这个存储弧段信息，弧段id，弧段是分方向的，弧度方向的左弧段、右弧段，弧段左边多边形、右边多边形，弧段坐标信息等。

start_node	end_node	edge_id	next_left_edge	abs_next_left_edge	next_right_edge	abs_next_right_edge	left_face	right_face	geom
开始节点	结束节点	边id	下一个左边	绝对下一个左边	下一个右边	绝对下一个右边	左面	右面	几何信息

2）face

face_id	bmr
多边形id	多边形外包矩形坐标信息

3）node

node_id	contiaining_face	geom
节点id	包含多边形	节点坐标信息

4）relation

topogeo_id	layer_id	element_id	element_type
拓扑id，建立的拓扑的id，存储在topology表中	图层id，建立的拓扑图层的id，存储在layer表中	要素id，对应空间数据表的gid	空间数据类型 （1:[multi]point, 2:[multi]line, 3:[multi]poly, 4:collection）

注：一般安装过postgis_topology，在数据库里都会一个topology模式，里面表有layer、topology表，此表里的信息不会自动更新，如果多次生成，记得清除之前存在的信息。

layer

topology_id	layer_id	schema_name	table_name	feature_column	level	child_id
拓扑id，与topology一致	图层id	模式名称	表名	要素列明	级别	子id

topology

id	name	srid	precision	hasz
拓扑id	拓扑名称	边id	精度

（6）为空间数据添加拓扑几何字段（只是添加拓扑几何字段，未生成信息）

添加拓扑几何字段需要传入的参数是：拓扑模式名称、数据所在模式名称、数据名称、拓扑几何字段名称和几何类型。这里将拓扑几何字段添加到复制的表(test_bak)上。

1	SELECT topology.AddTopoGeometryColumn('test_topo', 'public', 'test_bak', 'topo_geom', 'MULTIPOLYGON');

select * from test_bak;

（7）对空间数据拓扑几何字段生成拓扑信息（为每个要素建立拓扑关系）

构建拓扑传入的参数是：普通几何字段名称、拓扑模式名称、layer_id(为topolopy模式layer表图层id)、数据容差(经纬度数据使用默认容差0.000000008983153，平面数据使用0.001)：

1	UPDATE chn_adm1_bak SET topo_geom = topology.toTopoGeom(geom, 'test_topo', 1, 0.000000008983153);

这个处理过程和图形有关，由于空间数据节点非常密集，时间也比较长，处理完成之后：

 （8）拓扑弧段抽稀

接下来使用Sandro Santilli提供的SimplifyEdgeGeom函数进行拓扑弧段抽稀。

首先将SimplifyEdgeGeom函数添加到数据库中：

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

然后调用SimplifyEdgeGeom函数进行拓扑弧段抽稀，传入的参数是：拓扑模式名称、edge_id（为topo表edge里的edge字段，循环对每条记录抽稀）、抽稀容差（大地坐标系为经纬度距离：赤道上1度大约为111.3214km,1秒约为30.92米，如果想抽稀到200米一个点，约为0.001797度；投影坐标系：直接写200米）：

1	SELECT SimplifyEdgeGeom('test_topo', edge_id, 0.001797) FROM test_topo.edge;

（9）将抽稀结果更新普通几何字段中

1	UPDATE test_bak SET geom = topo_geom::geometry;

（10）抽稀结果检查

通过SQL查询抽稀前后朝阳区界的节点数：

1	select ST_NumPoints(ST_ExteriorRing(ST_GeometryN(geom,1))) from bjbj where name = '朝阳'

原始数据朝阳行政边界节点数9533个

抽稀后朝阳行政边界节点数287个

1	select ST_NumPoints(ST_ExteriorRing(ST_GeometryN(geom,1))) from bjbj_bak where name = '朝阳'

通过工具查看抽稀后省界间已不存在重叠和缝隙：

基于拓扑关系抽稀后数据

（11）后续工作

如果需要将抽稀后的结果导出，可以继续使用ogr2ogr工具进行导出。

整合后语句：

create table bj_bak as select * from bjbj; --bjbj为原始空间数据，bj_为空间数据备份

SELECT topology.CreateTopology('bj_topo', 4490); --4490为空间数据坐标系

SELECT topology.AddTopoGeometryColumn('bj_topo', 'public', 'bj_bak', 'topo_geom', 'MULTIPOLYGON'); --拓扑模式、数据所在模式、空间数据、拓扑几何列名，多面

UPDATE bj_bak SET topo_geom = topology.toTopoGeom(geom, 'bj_topo', 1, 0.000000008983153); --bj_bak空间数据，空间数据几何信息列，拓扑模式、图层id，容差

SELECT SimplifyEdgeGeom('bj_topo', edge_id, 0.001797) FROM bj_topo.edge; --topo模式，edge表edge_id列名，容差

UPDATE bj_bak SET geom = topo_geom::geometry; --bj_bak空间数据

下一次再抽稀其他空间数据时应用下列语句。

DROP SCHEMA bj_topo CASCADE; --级联清除topo模式

TRUNCATE layer CASCADE; --清空建立拓扑的图层记录

TRUNCATE topology CASCADE; --清空建立的拓扑记录