最短路径算法——dijkstra

先上代码

from osgeo import ogr
import numpy as np

def dijkstra(pList:dict,rList,adjacency,pStart):
    # pList 点表
    # 字典格式pList={pId:pIndex,…………}
    # rList 边表
    # 字典格式rList={rId:rLength,…………}
    # adjacency 空间邻接矩阵
    # 考虑到两节点间可能有多条边，所以采用三维的空间邻接矩阵
    # pStart 起点Id
    S=[pStart]
    T=list(pList.keys())
    T.remove(pStart)
    # 距离初始化
    n=len(T)
    dist=dict()
    track=dict()
    for i in range(n):
        pEnd=T[i]
        sindex=pList[pStart]
        eindex=pList[pEnd]
        num=adjacency[sindex,eindex,0]
        if num==0:
            dist[pEnd]=np.Infinity
            track[pEnd]=[]
            continue
        else:
            for j in range(1,int(num)+1):
                rid=adjacency[sindex,eindex,j]
                rdist=rList[rid]
                if j==1:
                    dist[pEnd]=rdist
                    track[pEnd]=[pStart,int(rid),pEnd]
                else:
                    if rdist<dist[pEnd]:
                        dist[pEnd]=rdist
                        track[pEnd]=[pStart,int(rid),pEnd]
    print(track)
    print('开始迭代')    
    # 循环迭代，每次中T中取出一个距离最小的点加入到S中,直到T为空
    sdist=dict()
    while len(T):
        mdist=min(dist.values())
        mindex=list(dist.values()).index(mdist)
        mid=list(dist.keys())[mindex]
        S.append(mid)
        T.remove(mid)
        sdist[mid]=dist[mid]
        del dist[mid]
        for i in range(len(T)):
            tid=T[i]
            mindex=pList[mid]
            tindex=pList[tid]
            num=adjacency[mindex,tindex,0]
            tmdist=np.Infinity
            tmrid=None
            if num==0:
                continue
            else:
                for j in range(1,int(num)+1):
                    rid=adjacency[mindex,tindex,j]
                    if rList[rid]<tmdist:
                        tmdist=rList[rid]
                        tmrid=rid
            if dist[tid]>(tmdist+sdist[mid]):
                # 更新最小距离 更新最短路径               
                dist[tid]=tmdist+sdist[mid]
                track[tid]=track[mid].copy()
                track[tid].extend([int(tmrid),int(tid)])
    print(track)
    print(sdist)
ogr.RegisterAll()
# 读取点表数据
pds=ogr.Open('./map/point.shp')
play=ogr.DataSource.GetLayerByIndex(pds,0)
ogr.Layer.ResetReading(play)
pList=dict()
pN=0
while True:
    fea=ogr.Layer.GetNextFeature(play)
    if not(fea):
        break
    fid=ogr.Feature.GetFieldAsInteger(fea,'id')
    pList[fid]=pN
    pN+=1
adjacency=np.zeros((pN,pN,10))
# 读取边表
pIds=pList.keys()
rList=dict()
rds=ogr.Open('./map/road.shp')
rlay=ogr.DataSource.GetLayerByIndex(rds,0)
ogr.Layer.ResetReading(rlay)
while True:
    fea=ogr.Layer.GetNextFeature(rlay)
    if not(fea):
        break
    fid=ogr.Feature.GetFieldAsInteger(fea,'id')
    start=ogr.Feature.GetFieldAsInteger(fea,'start')
    end=ogr.Feature.GetFieldAsInteger(fea,'end')
    sindex=pList[start]
    eindex=pList[end]
    adjacency[sindex,eindex,0]+=1
    rnum=adjacency[sindex,eindex,0]
    adjacency[sindex,eindex,int(rnum)]=fid
    adjacency[eindex,sindex,:]=adjacency[sindex,eindex,:]
    fgeo=ogr.Feature.geometry(fea)
    flength=ogr.Geometry.Length(fgeo)
    rList[fid]=flength

dijkstra(pList,rList,adjacency,4)

点shp文件

属性表中只有点Id

边shp文件

属性表中有边id以及每条边的两个端点

算法结果

• 输出的第一个字典是每个点到起点的最短路径，格式：点id——边id——点id——边id——……——终点id
• 输出的第二个字典是每条最短路径的长度，格式：点id：长度

与QGIS中的结果进行验证，路径完全一致，路径长度基本一致

粉红色表示的是4到5的最短路径：4——18——8——5——0——2——5,长度为：173400.070，也可以逐项进行其它点的验证