pgrouting3.0 AB经纬度 查找河流最短路径

因为百度地图 高德地图导航没有河流

所有自己下点地图

Geofabrik Download Server

然后切割选择上海，注意，线一定要转线段。然后倒入pg

建立拓扑关系 

ALTER TABLE shanghairiver ADD COLUMN "source" integer;
ALTER TABLE shanghairiver ADD COLUMN "target" integer;
 
--创建拓扑
SELECT pgr_createTopology('shanghairiver', 0.00001, 'geom', 'gid');
 
--为起点号终点号加空间索引
CREATE INDEX source_idx3 ON shanghairiver("source");
CREATE INDEX target_idx3 ON shanghairiver("target");
 
--添加长度字段、并计算赋值
ALTER TABLE shanghairiver ADD COLUMN length double precision;
update shanghairiver set length =st_length(geom);
 
--将长度值赋给reverse_cost，作为路线选择标准
ALTER TABLE shanghairiver ADD COLUMN reverse_cost double precision;
UPDATE shanghairiver SET reverse_cost = length;

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坐标
)
RETURNS TABLE (
    seq integer, -- 道路序号
    gid integer,
    name text, -- 道路名
    heading double precision,
    cost double precision, -- 消耗
    geom geometry -- 道路几何集合
) 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 source;

    -- 查询距离目的地最近的道路节点
    EXECUTE 'SELECT id::integer FROM ' || quote_ident(tbl) || '_vertices_pgr
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' || x2 || ' ' || y2 || ')'', 4326) LIMIT 1'
            INTO target;

    -- 最短路径查询
    seq := 0;
    sql := 'SELECT seq, gid, name, cost, geom FROM ' ||
           'pgr_dijkstra(''SELECT gid as id, source::integer, target::integer, '
           || 'length::float AS cost FROM ' || quote_ident(tbl) || ''', '
           || source || ', ' || target || ', false) as di '
           || 'LEFT JOIN ' || quote_ident(tbl) || ' ON di.edge = gid ORDER BY seq';

    -- Remember start point
    point := source;

    FOR rec IN EXECUTE sql
    LOOP
        -- Flip geometry (if required)
        IF (point != rec.seq) THEN
            rec.geom := ST_Reverse(rec.geom);
            point := rec.seq;
        ELSE
            point := rec.seq + 1;
        END IF;

        -- Check if rec.geom is not null before using it in EXECUTE
        IF rec.geom IS NOT NULL THEN
            EXECUTE 'SELECT degrees( ST_Azimuth(
                ST_StartPoint(''' || rec.geom::text || '''),
                ST_EndPoint(''' || rec.geom::text || ''') ) )'
                INTO heading;
        ELSE
            heading := NULL; -- Set heading as NULL if geom is NULL
        END IF;

        -- 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;
上面是最短路径

CREATE OR REPLACE FUNCTION pgr_fromAtoB_all_routes(
    IN tbl varchar, -- 数据库表名
    IN x1 double precision, -- 起点x坐标
    IN y1 double precision, -- 起点y坐标
    IN x2 double precision, -- 终点x坐标
    IN y2 double precision -- 终点y坐标
)
RETURNS TABLE (
    seq integer, -- 道路序号
    gid integer,
    name text, -- 道路名
    heading double precision,
    cost double precision, -- 消耗
    geom geometry -- 道路几何集合
) 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 source;

    -- 查询距离目的地最近的道路节点
    EXECUTE 'SELECT id::integer FROM ' || quote_ident(tbl) || '_vertices_pgr
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' || x2 || ' ' || y2 || ')'', 4326) LIMIT 1'
            INTO target;

    -- 所有路径查询
    seq := 0;
    sql := 'SELECT seq, gid, name, cost, geom FROM ' ||
           'pgr_ksp(''SELECT gid as id, source::integer, target::integer, '
           || 'length::float AS cost FROM ' || quote_ident(tbl) || ''', '
           || source || ', ' || target || ', 50) as ksp '  -- 指定最多返回 10 条路径
           || 'LEFT JOIN ' || quote_ident(tbl) || ' ON ksp.edge = gid ORDER BY seq';

    -- Remember start point
    point := source;

    FOR rec IN EXECUTE sql
    LOOP
        -- Flip geometry (if required)
        IF (point != rec.seq) THEN
            rec.geom := ST_Reverse(rec.geom);
            point := rec.seq;
        ELSE
            point := rec.seq + 1;
        END IF;

        -- Check if rec.geom is not null before using it in EXECUTE
        IF rec.geom IS NOT NULL THEN
            EXECUTE 'SELECT degrees( ST_Azimuth(
                ST_StartPoint(''' || rec.geom::text || '''),
                ST_EndPoint(''' || rec.geom::text || ''') ) )'
                INTO heading;
        ELSE
            heading := NULL; -- Set heading as NULL if geom is NULL
        END IF;

        -- 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;
上面是50条路径

select * from pgr_fromAtoB('shanghairiver'::text,      121.568589,31.290387,121.51029989,31.39239914)

select * from pgr_fromAtoB_all_routes('shanghairiver'::text,      121.568589,31.290387,121.51029989,31.39239914)