一 前言

  前几天，某科研工作者咨询我，是否可以利用PostGIS实现线数据的汇总分析，目前他使用的Arcpy处理实在太慢，难以满足应用需求。他的需求如下图：

输入数据：3条线路，权重值分别为1300，1200，1000。

输出数据：消除重叠部分线路，重叠部分聚合为单一线路段且权重值累加，结果输入数据被切分成了绝不重叠的5段。

本着好奇心的态度，利用周末时间对此问题进行了研究。

二 解决方案一（逐条截取法）

首先映入脑海的思路是，一条条的线路单独处理，如先求取1300路段和1200路段的不重叠部分a，和重叠部分b。

再用a,b这种被切分后的结果，与1000路段处理，a部分与1000路段无重叠，直接作为结果，b部分与1000路段重叠，被切割出不同部分b1和重叠部分c。

最终的结果，1300路段被切分的结果是a,b1,c。

同理处理1200,1000路段，最后对所有切割的图形进行geom的聚合运算，得到结果，重要算法如下：

#求取两个图形的不重叠部分select ST_Difference(geom1,geom2);#求取两个图形重叠部分select ST_Intersection(geom1,geom2);#根据图形聚合，权重累加，去除重叠图形select sum(weight) weight,geom from result group by geom;

整体代码逻辑如下：

do language plpgsql $$   DECLARE        rec record;        rec1 record;        rec2 record;        rec3 record;        rec4 record;        publicgeom geometry;        privategeom geometry;        geom_array geometry[];        geom_array_split geometry[];        temp_geom geometry;        _temp_geom geometry;        geomcount int;BEGIN    for rec in select * from test_road loop --遍历路段表        geom_array:=array[rec.geom];--切割成果集合初始化                --从表中找出与处理的路段存在压盖关系且非自身的数据        for rec1 in select gid,geom from test_road where ST_Overlaps(geom,rec.geom) and gid!=rec.gid loop            geom_array_split:='{}';                        --对每次的切割成果迭代分析，求取公共部分，单独部分            FOREACH temp_geom IN ARRAY geom_array loop                publicgeom:=ST_LineMerge(ST_Intersection(temp_geom,rec1.geom));--提取公共部分                privategeom:=ST_Difference(temp_geom,rec1.geom);--提取单独部分                if(ST_IsEmpty(publicgeom)=false) then                    for rec4 in select * from geomsplit(publicgeom) loop                        geom_array_split:=array_append(geom_array_split,rec4.geom);                    end loop;                end if;                if(ST_IsEmpty(privategeom)=false) then                    for rec4 in select * from geomsplit(privategeom) loop                        geom_array_split:=array_append(geom_array_split,rec4.geom);                    end loop;                end if;            end loop;            geom_array:=geom_array_split;--每次处理完重叠和不重叠部分，结果形成新的待切割结果        end loop;        --一条记录处理完毕，切割完毕的结果存入表        FOREACH temp_geom IN ARRAY geom_array loop            execute format('INSERT INTO result_test_road(            id, weight, geom) VALUES ($1,$2,$3)') using rec.gid,rec.weight,temp_geom;        end loop;    end loop;end;  $$; --对切割成果进行聚合分析select sum(weight),geom from result_test_road group by geom;

该算法在处理北京与各省会城市的路径时比较有效，样本数据30几条，很快处理完了。但是，处理北京与各县的路径时，数据量大概2500条，每切割完一条线路，消耗10分钟。（处理完大概333小时，oh my god）

总结：该算法思路比较简单，就是对每条线路逐个的切割，每次切割完的结果再和相交的线路切割，直到相交线路处理完毕。该算法保证了结果的正确性，逻辑的简单性，但是明显循环实在太多，计算代切割集合会像滚雪球一样几何级别的增长，尤其在数据表记录较多的情况下，简直惨不忍睹。

三 解决方案二（原子聚合法）

该方法是将每条线上的点拆解下来，重叠线路部分的点位也是重合的，然后根据点位进行权重聚合，相同位置的话，权重累加，如下图：

id为5的线路被拆分成了权重是6,13的两部分，同理，另外一个线路也根据权重实现了拆分。

代码如下：

--线拆分成点，该表存储所有节点CREATE TABLE road_pt(  gid serial primary key,  id integer,  weight numeric,  ptindex int,  geom geometry(Point,3857));--将每个线路拆分点存入road_pt表do language plpgsql $$   DECLARE        rec record;    ptcount int;BEGIN    for rec in select * from test_road loop        raise notice '正在处理gid:%',rec.gid;        ptcount:=ST_NPoints(rec.geom);         insert into road_pt(id,weight,ptindex,geom) values(rec.id,rec.weight,generate_series(1,ptcount),ST_PointN(rec.geom,generate_series(1,ptcount)));    end loop;end;  $$; --消耗50s--创建空间索引CREATE INDEX road_pt_geom_idx  ON road_pt  USING gist  (geom); --Query returned successfully with no result in 01:14 minutes.--相同的图形点，权重累加，权重聚合后的点数据存入表road_pt_groupselect sum(weight) weight,geom into road_pt_group from road_pt group by geom;--Query returned successfully: 76035 rows affected, 17.1 secs execution time.--创建空间索引CREATE INDEX road_pt_group_geom_idx  ON road_pt_group  USING gist  (geom);--3s--将点表中，与聚合图形相同的点位，统一更新新的点位权重update road_pt t2 set weight=t1.weight from road_pt_group t1 where ST_Equals(t2.geom,t1.geom);--Query returned successfully: 2237228 rows affected, 02:26 minutes execution time.--重新点拼接成线--点合成线，存储进入road_splitCREATE TABLE road_split(  gid serial primary key,  id integer,  weight numeric,  geom geometry(LineString,3857));--线的拆分点，合并成路段do language plpgsql $$   DECLARE        rec record;    i int;    ptcount int;    geomarr geometry[];    beforerec record;    isstart boolean:=true;BEGIN    for rec in select * from road_pt order by id,ptindex loop        if(isstart=true) then            beforerec:=rec;            isstart:=false;            geomarr:=array[rec.geom];            continue;        end if;        if(rec.id!=beforerec.id) then            insert into road_split(id,weight,geom) values (beforerec.id,beforerec.weight,st_makeline(geomarr));            geomarr:=array[rec.geom];        else            if(rec.weight>beforerec.weight) then                geomarr:=array_append(geomarr,rec.geom);                insert into road_split(id,weight,geom) values (rec.id,beforerec.weight,st_makeline(geomarr));                geomarr:=array[rec.geom];--开启新篇章            elsif(rec.weight

该算法的灵感认为：点的图形计算比线极大简化+空间索引的查询效率。请注意程序中每次的中间步骤表都有gist索引的创建！用于聚合，相等等判断快速完成。

最终，将每条线处理耗时10分钟（处理完333小时），优化到全部处理完成5min，且仍有可优化空间。

四 思考

• shp数据不应做项目应用，正如这位科研工作者的arcpy处理shp很慢的结果，参考。
• 批量空间数据处理，尤其大批量计算和分析，肯定应在空间数据引擎中实现而非后台，充分利用索引。
• 复杂分析，可能的话，可将图形抽象成点去分析，极大提升分析速度。
• 事情完成不代表做好，在某个数据量下可行不代表在任何情况下可行，思路要开阔，技术要多比较。
• 坚信GIS未来是空间数据的时代，PostGIS+大数据数据平台大有可为。