无后台最短路径分析实现

相关知识

​ LBS不仅需要能确定目标的地理位置，还需要能实现对地理环境的有效分析。网络分析是地理环境分析中的一个重要技术，包括最短路径分析、网络流分析等内容。在网络分析中，最短路径分析是最基本的，也是最关键的技术，一直是计算机科学、运筹学、交通工程学、地理信息学等学科的一个研究热点。

​ 路径分析是GIS中最基本的功能，其核心是对最佳路径和最短路径的求解。路径分析包含了两个基本内容：一个是路径的搜索；另一个是距离的计算。路径搜索的算法与连通分析是一致的，通过邻接关系的传递来实现路径搜索。路径的长度(距离)以积聚距离(accumulated distance)来计算。距离的计算方法为：将栅格路径视做由一系列路径段(path segments)组成，在进行路径搜索的同时计算每个路径段的长度并累计起来，表示从起点到当前栅格单元的距离。这里路径段指的是在一定的精度范围内可以以直线段模拟和计算的栅格单元集合。

技术路线

​ 本次要讲的无后台路径分析是指：仅适用空间数据库和WebGIS服务软件，再结合分析结果前端展示WebGIS客户端技术，实现路径数据的处理、存储、最短路径分析以及前端路网及最短路径分析结果展现。

​ 这里空间数据库使用PostgreSQL+PostGIS+pgRouting，WebGIS服务软件使用Geoserver，前端展示WebGIS库使用OpenLayers3。

PostgreSQL 是一个自由的对象-关系数据库服务器(数据库管理系统)。是世界上可以获得的最先进的开放源码的数据库系统，它提供了多版本并行控制，支持几乎所有 SQL 构件(包括子查询，事务和用户定 义类型和函数)，例如复杂SQL查询，SQL子选择，外键，触发器，视图，事务，多进程并发控制(MVCC)，流式复制(9.0)，热备(9.0))等。 并且可以获得非常广阔范围的(开发)语言绑定 (包括 C，C++，Java，perl，tcl，和 python)。PostgreSQL可在所有主要操作系统(即Linux，UNIX(AIX，BSD，HP-UX，SGI IRIX，Mac OS X，Solaris，Tru64)和Windows等)上运行。

尽管在PostgreSQL提供了上述几项支持空间数据的特性，但其提供的空间特性很难达到GIS的要求，主要表现在:缺乏复杂的空间类型;没有提供空间分析;没有提供投影变换功能。为了使得PostgreSQL更好的提供空间信息服务，PostGIS应运而生。构建在其上的空间对象扩展模块PostGIS则使得其成为一个真正的大型空间数据库。PostGIS在对象关系型数据库PostgreSQL上增加了存储管理空间数据的能力，相当于Oracle的spatial部分。PostGIS最大的特点是符合并且实现了几乎全部OpenGIS的规范，提供了空间对象、空间索引、空间操作函数和空间操作符等几乎所有商业数据库所具备的的空间数据库属性，PostGIS已经成为最著名的开源GIS数据库插件。

PgRouting是基于开源空间数据库PostGIS用于网络分析的扩展模块，最初它被称作pgDijkstra，因为它只是利用Dijkstra算法实现最短路径搜索，之后慢慢添加了其他的路径分析算法，如A算法，双向A算法，Dijkstra算法，双向Dijkstra算法，tsp算法等，然后被更名为pgRouting。该扩展库依托PostGIS自身的GIS索引，丰富的坐标系与图形类型，强大的几何处理能力，如空间查询，空间处理，线性参考等优势，能保障在较大数据级别下的网络分析效果更快更好。

GeoServer 是 OpenGIS Web 服务器规范的 J2EE 实现，利用 GeoServer 可以方便的发布地图数据，允许用户对特征数据进行更新、删除、插入操作，通过 GeoServer 可以比较容易的在用户之间迅速共享空间地理信息。Geoserver是著名的开源GIS软件之一。也是项目中常用的地图服务软件。兼容 WMS 和 WFS 特性；支持 PostgreSQL、 Shapefile 、 ArcSDE 、 Oracle 、 VPF 、 MySQL 、 MapInfo ；支持上百种投影；能够将网络地图输出为 jpeg 、 gif 、 png 、 SVG 、 KML 等格式；能够运行在任何基于 J2EE/Servlet 容器之上；嵌入 MapBuilder 支持 AJAX 的地图客户端OpenLayers；除此之外还包括许多其他的特性。

OpenLayers是一个用于开发WebGIS客户端的JavaScript包。OpenLayers 支持的地图来源包括Google Maps、Yahoo、 Map、微软Virtual Earth 等，用户还可以用简单的图片地图作为背景图，与其他的图层在OpenLayers 中进行叠加，在这一方面OpenLayers提供了非常多的选择。除此之外，OpenLayers实现访问地理空间数据的方法都符合行业标准。OpenLayers 支持Open GIS 协会制定的WMS（Web Mapping Service）和WFS（Web Feature Service）等网络服务规范，可以通过远程服务的方式，将以OGC 服务形式发布的地图数据加载到基于浏览器的OpenLayers 客户端中进行显示。OpenLayers采用面向对象方式开发，并使用来自Prototype.js和Rico中的一些组件。

操作步骤

软件安装

​ 本文不详述各软件的安装配置过程，只提供软件的下载链接：

GeoServer：http://geoserver.org/download/ 本文版本：2.10.1

PostgreSQL：https://www.postgresql.org/download/ 本文版本：9.6.2

PostGIS：http://postgis.net/windows_downloads/ 本文版本：2.3

pgAdmin III：https://www.pgadmin.org/download/pgadmin-3-windows/

注：pgRouting已经包含在PostGIS安装中,OpenLayers3可以使用CDN引入。另外，在数据制作和处理及导入数据库的过程中我用到的是QGIShttps://www.qgis.org/en/site/，考虑到大家对软件各有所爱，只要能达到效果，任何其他的GIS软件都可以，大家可以根据自己的情况，选择适合自己的GIS软件或工具，这里就不再做安装选项和详细说明。

数据制作处理

数据制作与检查

使用QGIS或其他GIS软件绘制测试路网数据。地理坐标系使用GCS_WGS_1984，编号为：4326。注意：所有线段连接的地方都需要断开，这样便于以后的分析。做好测试路网数据之后最好做一次拓扑检查，查看数据是否有拓扑问题，如果没有则可进行下一步处理。

安装数据库扩展

使用paAdmin3连接PostgreSQL数据库，连接或创建测试数据库，在测试数据库上打开SQL面板，执行以下语句，以在PostgreSQL数据库中添加空间扩展：
CREATE EXTENSION postgis; –创建数据库PostGIS扩展
CREATE EXTENSION pgrouting; –创建数据库pgRouting扩展
CREATE EXTENSION postgis_topology; –创建数据库PostGIS拓扑扩展
CREATE EXTENSION fuzzystrmatch; –字符串相似度计算扩展（配合下面的tiger使用）（本例可不装）
CREATE EXTENSION postgis_tiger_geocoder; –创建数据库地理编码扩展（本例可不装）
CREATE EXTENSION address_standardizer; –创建数据库地址标准化扩展（本例可不装）

数据入库

制作好数据并进行拓扑检查和处理后，再完成数据库空间扩展的安装，就可以将制作好检查好的数据导入PostgreSQL数据库了。这里我依然使用的是QGIS来进行数据入库。首先需要在QGIS中连接测试的PostgreSQL数据库作为数据源之一（左侧面板中选择PostgreSQL，输入参数即可），然后选择QGIS->DataBase->DB Manager，在DB Manager弹窗中选择PostGIS数据库选项，如果上一步PostgreSQL连接正确，即可在此处的PostGIS选项下看到已经连接的数据库名称了。选择此数据库下的Public架构，点击工具栏上面的Import按钮，选择当前加载的已经制作好的路网数据并设置好选项，点击OK，提示完成即表示已经导入到数据库。

数据处理

接下来我们需要对导入到数据库的空间表进行拓扑处理。打开pgAdmin并从测试数据库中打开SQL窗口，一次进行如下操作：

修改表结构

–添加起点id
ALTER TABLE public.road_new ADD COLUMN source integer;
–添加终点id
ALTER TABLE public.road_new ADD COLUMN target integer;
–添加道路权重值
ALTER TABLE public.road_new ADD COLUMN length double precision;

创建拓扑

–为road_new表创建拓扑布局，即为source和target字段赋值
SELECT pgr_createTopology(‘public.road_new’,0.00001, ‘geom’, ‘id’);
注：此处的gid可能实际表中为id，因此需要根据数据库中实际表字段给定

创建索引

–为source和target字段创建索引
CREATE INDEX source_idx ON road_new(“source”);
CREATE INDEX target_idx ON road_new(“target”);

权重字段赋值

–为length赋值
update road_new set length =st_length(geom);
–为road_new表添加reverse_cost字段并用length的值赋值
ALTER TABLE road_new ADD COLUMN reverse_cost double precision;
UPDATE road_new SET reverse_cost =length;

编写分析函数

通过以上操作，所有的数据准备工作已经完成，接下来我们需要编写数据库函数，来实现任意两点之间的最短路径计算。打开pgAdmin，在测试数据库中打开SQL查询窗口，开始编写路径分析函数：

DROP FUNCTION pgr_fromPointtoPoint(tbl varchar,startx float, starty float,endx float,endy float);

CREATE OR REPLACE function pgr_fromPointtoPoint(tbl varchar,startx float, starty float,endx float,endy float)  

returns geometry as 

$body$  

declare 

    v_startLine geometry;--离起点最近的线 

    v_endLine geometry;--离终点最近的线 


    v_startTarget integer;--距离起点最近线的终点

    v_startSource integer;

    v_endSource integer;--距离终点最近线的起点

    v_endTarget integer;


    v_statpoint geometry;--在v_startLine上距离起点最近的点 

    v_endpoint geometry;--在v_endLine上距离终点最近的点 


    v_res geometry;--最短路径分析结果

    v_res_a geometry;

    v_res_b geometry;

    v_res_c geometry;

    v_res_d geometry; 


    v_perStart float;--v_statpoint在v_res上的百分比 

    v_perEnd float;--v_endpoint在v_res上的百分比 
 

    v_shPath_se geometry;--开始到结束

    v_shPath_es geometry;--结束到开始

    v_shPath geometry;--最终结果

    tempnode float;      
begin

    --查询离起点最近的线 
    execute 'select geom, source, target  from ' ||tbl||
                            ' where ST_DWithin(geom,ST_Geometryfromtext(''point('||         startx ||' ' || starty||')'',4326),15)
                            order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',4326))  limit 1'
                            into v_startLine, v_startSource ,v_startTarget; 
     

    --查询离终点最近的线 
    execute 'select geom, source, target from ' ||tbl||
                            ' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')'',4326),15)
                            order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'',4326))  limit 1'
                            into v_endLine, v_endSource,v_endTarget; 
 

    --如果没找到最近的线，就返回null 
    if (v_startLine is null) or (v_endLine is null) then 
        return null; 
    end if ; 

    select  ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')',4326)) into v_statpoint; 
    select  ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')',4326)) into v_endpoint;    

   -- ST_Distance 
     

    --从开始的起点到结束的起点最短路径 
    execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
    'FROM pgr_kdijkstraPath( 
    ''SELECT id as id, source, target, length as cost FROM ' || tbl ||''',' 
    ||v_startSource || ', ' ||'array['||v_endSource||'] , false, false 
    ) a, ' 
    || tbl || ' b 
    WHERE a.id3=b.id   
    GROUP by id1   
    ORDER by id1' into v_res ;
   

    --从开始的终点到结束的起点最短路径
    execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
    'FROM pgr_kdijkstraPath( 
    ''SELECT id as id, source, target, length as cost FROM ' || tbl ||''',' 
    ||v_startTarget || ', ' ||'array['||v_endSource||'] , false, false 
    ) a, ' 
    || tbl || ' b 
    WHERE a.id3=b.id   
    GROUP by id1   
    ORDER by id1' into v_res_b ;
 

    --从开始的起点到结束的终点最短路径
    execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
    'FROM pgr_kdijkstraPath( 
    ''SELECT id as id, source, target, length as cost FROM ' || tbl ||''',' 
    ||v_startSource || ', ' ||'array['||v_endTarget||'] , false, false 
    ) a, ' 
    || tbl || ' b 
    WHERE a.id3=b.id   
    GROUP by id1   
    ORDER by id1' into v_res_c ;
 

    --从开始的终点到结束的终点最短路径
    execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
    'FROM pgr_kdijkstraPath( 
    ''SELECT id as id, source, target, length as cost FROM ' || tbl ||''',' 
    ||v_startTarget || ', ' ||'array['||v_endTarget||'] , false, false 
    ) a, ' 
    || tbl || ' b 
    WHERE a.id3=b.id   
    GROUP by id1   
    ORDER by id1' into v_res_d ;

    if(ST_Length(v_res) > ST_Length(v_res_b)) then
       v_res = v_res_b;
    end if;

    if(ST_Length(v_res) > ST_Length(v_res_c)) then
       v_res = v_res_c;
    end if;

    if(ST_Length(v_res) > ST_Length(v_res_d)) then
       v_res = v_res_d;
    end if;

    --如果找不到最短路径，就返回null 
    --if(v_res is null) then 
    --    return null; 
    --end if; 

    --将v_res,v_startLine,v_endLine进行拼接 
    select  st_linemerge(ST_Union(array[v_res,v_startLine,v_endLine])) into v_res;

    select  ST_Line_Locate_Point(v_res, v_statpoint) into v_perStart; 
    select  ST_Line_Locate_Point(v_res, v_endpoint) into v_perEnd;         

    if(v_perStart > v_perEnd) then 
        tempnode =  v_perStart;
        v_perStart = v_perEnd;
        v_perEnd = tempnode;
    end if;

    --截取v_res 
    SELECT ST_Line_SubString(v_res,v_perStart, v_perEnd) into v_shPath;

    return v_shPath; 

end; 

$body$ 

LANGUAGE plpgsql VOLATILE STRICT;

slqView服务发布

完成以上步骤后，通过Geoserver SQLView方式发布路网分析服务，过程如下：

新建工作空间

创建数据存储

新建SQLView图层

前端结果展示

前端展示代码

<!DOCTYPE html>
<html lang="en">
<head>
    <meta charset="UTF-8">
    <title>路径分析</title>
    <link rel="stylesheet" href="./css/ol.css" type="text/css">
    <script src="./js/ol.js"></script>
    <style>
        #map:focus {
            outline: #4A74A8 solid 0.15em;
        }
    </style>
</head>
<body>
    <div style="height: 50px;margin-top: 0px;margin-left: 0px">
        <div style="margin-top: 20px">
            <div>
                <input type="text" id="startPoint"/>
                <input type="text" id="endPoint"/>
            </div>
            <div style="float: right;">
                <button onclick="startNavigation()">任意路径规划</button>
                <button onclick="clearResult()">清除路径规划</button>
            </div>
        </div>
    </div>
    <div id="map" class="map" tabindex="0"></div>
    
<script>

    //初始化功能参数
    var params = {
        LAYERS: 'test:routeway',
        FORMAT: 'image/png'
    };
    var routeLayer = null;
    var navigateLayer = null;
    var startPoint;
    var destPoint;
    var vectorLayer;

    //初始化地图和底图
    var map = new ol.Map({
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ],
        target: 'map',
        controls: ol.control.defaults({
            attributionOptions: {
                collapsible: false
            }
        }),
        view: new ol.View({
            projection: 'EPSG:4326',
            center: [lontitude,latitude],
            zoom: 15
        })
    });

    //添加路网基础服务
    var roadLayer = new ol.layer.Tile({
        source: new ol.source.TileWMS({
        url: 'http://localhost:8097/geoserver/yc/wms',
        params: {'LAYERS': 'test:road', 'TILED': true},
            serverType: 'geoserver'
        })
    });
    map.addLayer(roadLayer);


    function startNavigation() {
        clearResult();

        startPoint = new ol.Feature();
        destPoint = new ol.Feature();

        vectorLayer = new ol.layer.Vector({
            source: new ol.source.Vector({
                features: [startPoint, destPoint]
            }),
            style:new ol.style.Style({
                image:new ol.style.Icon(({
                    size:[24,36],
                    anchor:[0.5,0.75],
                    anchorXUnits:'fraction',
                    anchorYUnits:'fraction',
                    src:'./img/marker.png'
                }))
            })
        });
        map.addLayer(vectorLayer);
        map.on('click', clickMap);
    }

    function initClearButtion(type) {
        var clearButton = document.getElementById('clear');
        clearButton.addEventListener('click', function(event) {
            clearResult(type);
        });
    }

    function clearResult(type) {
        clearNavigateLayer();
    }

    function clearNavigateLayer() {
        if(navigateLayer){
            startPoint.setGeometry(null);
            destPoint.setGeometry(null);
            map.removeLayer(vectorLayer);
            map.removeLayer(navigateLayer);
            navigateLayer = null;
        }
    }

    function clickMap(event) {
        if (startPoint.getGeometry() == null && destPoint.getGeometry() == null) {
            // First click.
            startPoint.setGeometry(new ol.geom.Point(event.coordinate));
        }else if (startPoint.getGeometry() != null && destPoint.getGeometry() == null) {
            // Second click.
            destPoint.setGeometry(new ol.geom.Point(event.coordinate));

            var startCoord = (startPoint.getGeometry().getCoordinates());
            var destCoord = (destPoint.getGeometry().getCoordinates());
            console.log(startCoord,destCoord);

            navigateLayer = new ol.layer.Tile({
                source: new ol.source.TileWMS({
                    url: 'http://localhost:8097/geoserver/test/wms',
                    params: {'LAYERS': 'test:routeway', 'TILED': true,'viewparams':'sx:'+startCoord[0]+';sy:'+startCoord[1]+';ex:'+destCoord[0]+';ey:'+destCoord[1]},
                    serverType: 'geoserver'
                })
            });
            console.log(navigateLayer);
            map.addLayer(navigateLayer);
            map.un('click', clickMap);
        }
    }
    
</script>
</body>
</html>

前端展示效果