# 利用python的networkx3.0 进行GIS的网络分析

## 第一步

``````# import we must need library files
# the function of 'as' is samplify quote
import geopandas
from tqdm import tqdm
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from shapely import geometry
import os
``````

ikzhttps://www.osgeo.cn/networkx/reference/algorithms/shortest_paths.htmld

## 第二步

``````G[(1,1)] # 会发生key值错误，原因在于实际上的key值是“（1，1）”，有双引号
# 因此我们还是利用整数索引当作节点标志为好``````
``````"""
this code will show how to translate shape files into networks through the python Networkx library
we will use the file of hefei province' road network to conduct the experiment
"""

# 计算欧式距离 利用np矩阵算法
def calEuclidean(x, y):
x = np.array(x, dtype='f4')
y = np.array(y, dtype='f4')
dist = np.sqrt(np.sum(np.square(x - y)))
return dist

# 生成导航文本
with open(file=r"导航.txt", mode="w", encoding="utf-8") as f:
for i in range(len(short_path) - 1):
f.write("从{0}点到{1}点，预计用时{2}s，行驶距离为{3}m n".format(short_path[i], short_path[i + 1],
network[short_path[i]][short_path[i + 1]]['time'],
network[short_path[i]][short_path[i + 1]]['weight']))

# 计算网络中离中心点距离最小的点索引，参数分别是点文件，和停靠距离
def nearest(filename, dist):
arry = []
series = geopandas.GeoSeries.from_file(filename, crs="EPSG:3395", index=[1, 2])
points = geopandas.GeoSeries.from_file("node.shp", crs="EPSG:3395")
for i in range(len(series.geometry)):
arry.append(0)
point = series.geometry[i]
distance = dist
for n in range(len(points.geometry)):
# 先确定在给定矩形内的点
if (point.x - dist < points.geometry[n].x <= point.x + dist) and (
point.y - dist < points.geometry[n].y <= point.y + dist):
# 在计算这些点当中最近的点
if calEuclidean([points.geometry[n].x, points.geometry[n].y], [point.x, point.y]) < distance:
arry.pop()
arry.append(n)
distance = calEuclidean([points.geometry[n].x, points.geometry[n].y], [point.x, point.y])
if arry[i] == 0:
print("没有找到第{}个点的停靠点,可能距离过小".format(i + 1))
return arry

# 读取线要素的shp文件，转换为networkx支持的网络图
class convert_to_binary(object):
def __init__(self, src_path):
self.coordinate = []
self.list = []
self.number = 0
self.number2 = 0
# 创建一个无向图 这个类也只基于无向图
self.network = nx.MultiGraph()
self.length = self.network_data.shape[0]
self.filename = src_path
self.temporary_coordinates = []

# SHP file converted to network data
def shp_network(self):
if os.path.exists("test.gexf"):
print("已经成功创建网络文件")
if os.path.exists("node.npy"):
print("正在读取坐标文件")
for i in zip(coordinate[:, 0], coordinate[:, 1]):
self.coordinate.append(i)
return self.network
else:
# 这种挨个搜索的方法只适合建立无向图，毕竟只是手搓的，功能有限，而且速度有待改进
for n in tqdm(range(self.length), desc="构建网络数据集中："):
def Generate_network(data):
line_len = len(data.xy[0])
for ii in zip(data.xy[0], data.xy[1]):
if ii not in self.coordinate:
self.coordinate.append(ii)
self.list.append(ii)
for m in range(line_len - 1):
# 确保不会出现一个位置上有两个端点
start = self.number + m
end = self.number + m + 1
if self.list[start] in self.coordinate:
start = self.coordinate.index(self.list[start])
if self.list[end] not in self.coordinate:
end = self.number2 + m
else:
end = self.coordinate.index(self.list[end])
else:
start = self.number2 + m
if self.list[end] in self.coordinate:
end = self.coordinate.index(self.list[end])
else:
end = self.number2 + m + 1
distance = calEuclidean(self.coordinate[start], self.coordinate[end])
# 赋予每条边的权重为距离和时间，时间等于距离除以速度，3395的单位是m，速度我设置的是m/s，因此时间单位是s
# 路网在处理之初，要对所有道路进行分类和加上速度属性，一般不同的道路对速度要求不一
tm = distance / float(self.network_data['avgspeed'][n])
self.number = self.number + line_len
self.number2 = len(self.coordinate)

lin_str = self.network_data.geometry[n]
if lin_str.type == 'LineString':
Generate_network(lin_str)
else:
for mul_line in lin_str:
Generate_network(mul_line)
nx.write_gexf(self.network, "test.gexf")
self.Calculate_coordinates()
return self.network

# 返回索引对应的xy坐标元组
def getpos(self):
return self.coordinate

# 遍历坐标
def Calculate_coordinates(self):
for n in tqdm(range(self.length)):
data = self.network_data.geometry[n]
for m in zip(data.xy[0], data.xy[1]):
if m not in self.coordinate:
self.coordinate.append(m)
# 保存所有节点
coordinate = np.array(self.coordinate)
points = geopandas.GeoSeries(geopandas.points_from_xy(coordinate[:, 0], coordinate[:, 1], crs="EPSG:3395"))
points.to_file("node.shp")
np.save("node.npy", coordinate)

# 保存路径
def save_path(self, path):
self.temporary_coordinates = []
for i in path:
self.temporary_coordinates.append(self.coordinate[int(i)])
lines = geopandas.GeoSeries([geometry.LineString(xy for xy in self.temporary_coordinates)],
crs='EPSG:3395')
lines.to_file("path.shp")
print("路径保存成功")``````

## 第三步

主函数部分：

``````if __name__ == "__main__":
# select GDK as the encoding format
# 查看是不是投影坐标系
print("The projection information of {}  is {}".format(src, network_road.crs))
# 仅仅只是我写的这个函数是针对线要素的，你可以构想对于多边形和单点的,(此外线要素分单线与多线)
print("仅支持线要素")
# an undirected graph is built to read from SHP file
print("the number of shapes are:%d" % network_road.shape[0])
print(network.number_of_edges())
print(network.number_of_nodes())
# 因为不想一次运行太多，我把代码注释了一部分
# 在下方写需要的功能
# 使用matplotlib库绘制
options = {"node_color": "blue", "node_size": 0.5, "linewidths": 0, "width": 0.1}
nx.draw(network, pos=pos, **options)
plt.show()``````

运行中：

查看plt的绘图效果：

效果不是很好（默默打开ArcGIS），把节点的size调为0，边的宽度为1

``````# source=500, target=5000，前者是起点标志，后者是终点标志
short_path = nx.shortest_path(network, source=500, target=5000)
with open(file=r"导航.txt", mode="w", encoding="utf-8") as f:
for i in range(len(short_path) - 1):
f.write("从{0}点到{1}点，预计用时{2}s，行驶距离为{3}m n".format(short_path[i], short_path[i + 1],
network[short_path[i]][short_path[i + 1]]['time'],
network[short_path[i]][short_path[i + 1]]['weight']))
# 保存节点和路径数据

运行后生成一个node.shp文件，这个是用来存储我们文件中的所有不重复节点，打开看看。

看起来确实是距离最短，接下来分析时间最短路径

``  short_path = nx.shortest_path(network, source=500, target=5000,weight="time")``

我们知道那条路其实是高架，所以也很有道理：）用时缺少最短，当然我们只分析的是开车的情况。

## 最后一步

``````# 计算网络中离中心点距离最小的点索引，参数分别是点文件，和停靠距离
def nearest(filename, dist):
arry = []
series = geopandas.GeoSeries.from_file(filename, crs="EPSG:3395", index=[1, 2])
points = geopandas.GeoSeries.from_file("node.shp", crs="EPSG:3395")
for i in range(len(series.geometry)):
arry.append(0)
point = series.geometry[i]
distance = dist
for n in range(len(points.geometry)):
# 先确定在给定矩形内的点
if (point.x - dist < points.geometry[n].x <= point.x + dist) and (
point.y - dist < points.geometry[n].y <= point.y + dist):
# 在计算这些点当中最近的点
if calEuclidean([points.geometry[n].x, points.geometry[n].y], [point.x, point.y]) < distance:
arry.pop()
arry.append(n)
distance = calEuclidean([points.geometry[n].x, points.geometry[n].y], [point.x, point.y])
if arry[i] == 0:
print("没有找到第{}个点的停靠点,可能距离过小".format(i + 1))
return arry``````

``````# 因为3395坐标系的单位是米，因此500的含义为500m，也就是停靠范围，是在500m以内的最近点
park = nearest("路网/endpoint.shp", 500)
print("数组是{}".format(park))
# source 是起点标志，而target是终点标志
short_path = nx.shortest_path(network, source=park[0], target=park[1])
# 保存路径数据