PostGIS pgrouting路徑分析

--讓數據庫支持PostGIS和pgRouting的函數和基礎表（安裝后第一次使用時執行，以后都不再執行）
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
CREATE EXTENSION address_standardizer; 

--注：使用postgis shapefile mport/export manager上傳shp時，在Option中勾選“generate simple geometries instead of multi geometries”，以生成單個geometry

-----------------------------------------
--以節點為參數進行最短路徑分析
--以road表作為實例--
-----------------------------------------
ALTER TABLE road ADD COLUMN source integer;--起點 
ALTER TABLE road ADD COLUMN target integer;--終點
ALTER TABLE road ADD COLUMN length double precision;--增加路線長度字段（根據長度設置權重）
UPDATE road SET length = ST_Length(geom);--計算路線長度


select pgr_createTopology('road', 0.0001, 'geom', 'gid');--創建拓撲


SELECT * FROM pgr_dijkstra('SELECT gid as id,
source::integer,target::integer,length::double precision as cost 
FROM road',30, 60, false, false); --路徑分析


SELECT st_astext(geom) FROM pgr_dijkstra('SELECT gid AS id,source::integer,
target::integer,length::double precision AS cost FROM road',30, 60, false, false) as di
join road pt on di.id2 = pt.gid;--查詢所經過的所有點


SELECT seq, id1 AS node, id2 AS edge, cost,geom into dijkstra_res FROM pgr_dijkstra('
SELECT gid AS id,source::integer,target::integer,length::double precision AS cost FROM road',
30, 60, false, false) as di join road pt on di.id2 = pt.gid;--查詢結果存儲到新的表格

select * from dijkstra_res;--查詢表格內容

-----------------------------------------
--以起始點坐標為參數進行最短路徑分析
--以gyroad表作為實例--
-----------------------------------------
ALTER TABLE gyroad ADD COLUMN source integer;--起點 
ALTER TABLE gyroad ADD COLUMN target integer;--終點
ALTER TABLE gyroad ADD COLUMN length double precision;--增加路線長度字段（根據長度設置權重）
UPDATE gyroad SET length = ST_Length(geom);--計算路線長度
select pgr_createTopology('gyroad', 0.0001, 'geom', 'gid');--創建拓撲

--添加起始點坐標x,y字段
ALTER TABLE gyroad ADD COLUMN x1 double precision;
ALTER TABLE gyroad ADD COLUMN y1 double precision;
ALTER TABLE gyroad ADD COLUMN x2 double precision;
ALTER TABLE gyroad ADD COLUMN y2 double precision;
--計算起始點坐標
UPDATE gyroad SET x1 =ST_x(ST_PointN(geom, 1));
UPDATE gyroad SET y1 =ST_y(ST_PointN(geom, 1));
UPDATE gyroad SET x2 =ST_x(ST_PointN(geom, ST_NumPoints(geom)));
UPDATE gyroad SET y2 =ST_y(ST_PointN(geom, ST_NumPoints(geom)));

SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_astar('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost,
x1, y1, x2, y2 FROM gyroad',30, 60, false,false);--A*算法路徑查詢

SELECT seq, id1 AS source, id2 AS target,cost FROM pgr_kdijkstraCost('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost FROM gyroad',30, 
array[60,70,100], false, false);--查詢從出發點到目的地的消耗

SELECT seq, id1 AS path, id2 AS edge, cost FROM pgr_kdijkstraPath('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost
FROM gyroad',30, array[60,100], false, false);--pgr_kdijkstraPath函數路徑分析

-----------------------------------------
--DROP FUNCTION pgr_fromAtoB(varchar, double precision, double precision, 
--                           double precision, double precision);
--基於任意兩點之間的最短路徑分析
CREATE OR REPLACE FUNCTION pgr_fromAtoB(
                IN tbl varchar,--數據庫表名
                IN x1 double precision,--起點x坐標
                IN y1 double precision,--起點y坐標
                IN x2 double precision,--終點x坐標
                IN y2 double precision,--終點y坐標
                OUT seq integer,--道路序號
                OUT gid integer,
                OUT name text,--道路名
                OUT heading double precision,
                OUT cost double precision,--消耗
                OUT geom geometry--道路幾何集合
        )
        RETURNS SETOF record AS
$BODY$
DECLARE
        sql     text;
        rec     record;
        source    integer;
        target    integer;
        point    integer;
        
BEGIN
    -- 查詢距離出發點最近的道路節點
    EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr 
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' 
            || x1 || ' ' || y1 || ')'',4326) LIMIT 1' INTO rec;
    source := rec.id;
    
    -- 查詢距離目的地最近的道路節點
    EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr 
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' 
            || x2 || ' ' || y2 || ')'',4326) LIMIT 1' INTO rec;
    target := rec.id;

    -- 最短路徑查詢 
        seq := 0;
        sql := 'SELECT gid, geom, name, cost, source, target, 
                ST_Reverse(geom) AS flip_geom FROM ' ||
                        'pgr_bdAstar(''SELECT gid as id, source::int, target::int, '
                                        || 'length::float AS cost,x1,y1,x2,y2 FROM '
                                        || quote_ident(tbl) || ''', '
                                        || source || ', ' || target 
                                        || ' ,false, false), '
                                || quote_ident(tbl) || ' WHERE id2 = gid ORDER BY seq';


    -- Remember start point
        point := source;

        FOR rec IN EXECUTE sql
        LOOP
        -- Flip geometry (if required)
        IF ( point != rec.source ) THEN
            rec.geom := rec.flip_geom;
            point := rec.source;
        ELSE
            point := rec.target;
        END IF;

        -- Calculate heading (simplified)
        EXECUTE 'SELECT degrees( ST_Azimuth( 
                ST_StartPoint(''' || rec.geom::text || '''),
                ST_EndPoint(''' || rec.geom::text || ''') ) )' 
            INTO heading;

        -- Return record
                seq     := seq + 1;
                gid     := rec.gid;
                name    := rec.name;
                cost    := rec.cost;
                geom    := rec.geom;
                RETURN NEXT;
        END LOOP;
        RETURN;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-----------------------------------------
--測試
SELECT st_astext(ST_MakeLine(route.geom)) FROM (SELECT seq,gid,name,heading,cost,geom FROM pgr_fromAtoB('gyroad', 106.535, 26.905, 106.955, 27.040)ORDER BY seq) AS route
--Openlayers測試
--http://localhost:6060/geoserver/PostGIS/wms?service=WMS&version=1.1.0&request=GetMap&layers=PostGIS:shortgyroad&styles=&bbox=104,24,108,28&width=330&height=768&srs=EPSG:4326&format=application/openlayers&viewparams=x1:106.565;y1:26.915;x2:106.925;y2:28.040