from libpysal import weights, examples
from libpysal.cg import voronoi_frames
from contextily import add_basemap
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import geopandas
import itertools
# 输入数据
cases = geopandas.read_file(r"E:\code_practice\gitee\Python地理空间分析指南\nx_study\networkx\examples\geospatial\cholera_cases.gpkg")
cases
.dataframe tbody tr th:only-of-type { vertical-align: middle }
\3cpre>\3ccode>.dataframe tbody tr th { vertical-align: top }
.dataframe thead th { text-align: right }
Id
Count
geometry
0
0
1
POINT (-15539.921 6712903.938)
1
0
3
POINT (-15356.721 6712858.073)
2
0
2
POINT (-15336.090 6712866.374)
3
0
1
POINT (-15353.877 6713017.559)
4
0
0
POINT (-15282.597 6712895.684)
...
...
...
...
319
0
0
POINT (-14951.181 6712211.606)
320
0
0
POINT (-14979.678 6712225.369)
321
0
3
POINT (-15591.770 6712227.697)
322
0
0
POINT (-15525.087 6712152.522)
323
0
0
POINT (-14952.725 6712116.692)
324 rows × 3 columns
为了让nx正确绘制图的节点，我们需要为数据集中的每个点构造坐标数组
从geometry列中提取x,y坐标
# x坐标
x = cases.geometry.x
x
0 -15539.921064
1 -15356.720566
2 -15336.089980
3 -15353.877335
4 -15282.596809
...
319 -14951.180762
320 -14979.678278
321 -15591.770010
322 -15525.086590
323 -14952.725267
Length: 324, dtype: float64
# y
y= cases.geometry.y
y
0 6.712904e+06
1 6.712858e+06
2 6.712866e+06
3 6.713018e+06
4 6.712896e+06
...
319 6.712212e+06
320 6.712225e+06
321 6.712228e+06
322 6.712153e+06
323 6.712117e+06
Length: 324, dtype: float64
nodes_coords = np.column_stack((x,y))
#切片展示
nodes_coords[:5]
array([[ -15539.92106438, 6712903.93781305],
[ -15356.72056559, 6712858.07280407],
[ -15336.08997978, 6712866.3743291 ],
[ -15353.877335 , 6713017.55898743],
[ -15282.5968086 , 6712895.68378701]])
泰森多边形的建立 (具体看libpysal库)
并将
region_df, point_df = voronoi_frames(nodes_coords, clip="convex hull")
fig, ax = plt.subplots()
region_df.plot(ax=ax, color='blue',edgecolor='black', alpha=0.3)
point_df.plot(ax=ax, color='red',markersize =3)
利用泰森多边形，可以利用“rook”临接性来构造他们之间的邻接图
如果泰森多边形共享边/面，则它们是相邻的
delaunay = weights.Rook.from_dataframe(region_df)
delaunay_graph = delaunay.to_networkx()
nx.draw(delaunay_graph,node_size =10)
# 给节点匹配coords
positions = dict(zip(delaunay_graph.nodes, nodes_coords))
#切片展示
positions_head = dict(itertools.islice(positions.items(), 5))
positions_head
{0: array([ -15539.92106438, 6712903.93781305]),
1: array([ -15356.72056559, 6712858.07280407]),
2: array([ -15336.08997978, 6712866.3743291 ]),
3: array([ -15353.877335 , 6713017.55898743]),
4: array([ -15282.5968086 , 6712895.68378701])}
add_basemap :获取地理底图
positions ：dict，节点-坐标，用于可视化
# 可视化
# 用nx可视化delaunay_graph,并展示节点正确的地理位置
ax = region_df.plot(facecolor="lightblue", alpha=0.50, edgecolor="cornsilk", linewidth=2)
add_basemap(ax)
ax.axis("off")
nx.draw(
delaunay_graph,
positions, #以节点为键，坐标为值
ax=ax,
node_size=2,
node_color="k",
edge_color="k",
alpha=0.8,
)
plt.show()
#下面是直接用点、多边形数据可视化的结果，可以与上面对比一下
region_df, point_df = voronoi_frames(nodes_coords, clip="convex hull")
fig, ax = plt.subplots()
region_df.plot(ax=ax, color='blue',edgecolor='black', alpha=0.3)
point_df.plot(ax=ax, color='red',markersize =3)
发表评论