无题

In [1]:

import networkx as nx
import osmnx as ox
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
ox.__version__

Out[1]:

'1.1.1'

基本展示¶

In [5]:

city = ox.geocode_to_gdf('Qinhuai,Nanjing, China')
ax = ox.project_gdf(city).plot()
#_ = ax.axis('off')

In [3]:

# G = ox.graph_from_place('Qinhuai, Nanjing, China', network_type='all')
streets=gpd.read_file("basic/road21Final.shp")
streets.head()

Out[3]:

	FID_road21	type	name	lengthQGIS	ID	geometry
0	13	lblue	干线公路	746.941	0	LINESTRING Z (683092.608 3454982.000 0.000, 68...
1	12	lblue	干线公路	684.967	1	LINESTRING Z (683649.578 3455380.702 0.000, 68...
2	33	red	高速公路	2212.137	2	LINESTRING Z (651690.838 3456004.237 0.000, 65...
3	0	red	高速公路	2494.064	3	LINESTRING Z (687071.291 3456735.217 0.000, 68...
4	11	lblue	干线公路	1826.529	4	LINESTRING Z (684644.129 3456912.720 0.000, 68...

In [4]:

for i in range(len(streets)):
    streets.geometry[i]=shapely.ops.transform(lambda x, y, z=None: (x, y), streets.geometry[i])  # 'LINESTRING (2 3, 4 5)'
    
graph = momepy.gdf_to_nx(streets)

In [5]:

nodes, edges = momepy.nx_to_gdf(graph)
nodes.head()

Out[5]:

	nodeID	geometry
0	0	POINT (683092.608 3454982.000)
1	1	POINT (682886.751 3454263.987)
2	2	POINT (683649.578 3455380.702)
3	3	POINT (651690.838 3456004.237)
4	4	POINT (652171.081 3453844.858)

In [6]:

edges.head()

Out[6]:

	FID_road21	type	name	lengthQGIS	ID	geometry	mm_len	node_start	node_end
0	13	lblue	干线公路	746.941	0	LINESTRING (683092.608 3454982.000, 682886.751...	746.940674	0	1
1	12	lblue	干线公路	684.967	1	LINESTRING (683649.578 3455380.702, 683092.608...	684.966555	0	2
2	11	lblue	干线公路	1826.529	4	LINESTRING (684644.129 3456912.720, 683649.578...	1826.528798	2	7
3	33	red	高速公路	2212.137	2	LINESTRING (651690.838 3456004.237, 652171.081...	2212.136803	3	4
4	32	red	高速公路	5939.964	35	LINESTRING (648142.369 3460767.801, 651690.838...	5939.963801	3	37

In [7]:

nodes["osmid"]=nodes["nodeID"]
nodes["x"]=nodes.geometry.x
nodes["y"]=nodes.geometry.y
nodes=nodes[['osmid','y','x','nodeID','geometry']]

In [19]:

edges["u"]=edges["node_start"]
edges["v"]=edges["node_end"]
edges["osmid"]=edges["ID"]
edges["key"]=0
edges['length']=edges['mm_len']
edges=edges[['u','v','key','osmid','FID_road21','type','name','lengthQGIS','ID','geometry','mm_len','length','node_start','node_end']]

Out[19]:

			u	v	key	osmid	FID_road21	type	name	lengthQGIS	ID	geometry	mm_len	length	node_start	node_end
u	v	key
0	1	0	0	1	0	0	13	lblue	干线公路	746.941	0	LINESTRING (683092.608 3454982.000, 682886.751...	746.940674	746.940674	0	1
	2	0	0	2	0	1	12	lblue	干线公路	684.967	1	LINESTRING (683649.578 3455380.702, 683092.608...	684.966555	684.966555	0	2
2	7	0	2	7	0	4	11	lblue	干线公路	1826.529	4	LINESTRING (684644.129 3456912.720, 683649.578...	1826.528798	1826.528798	2	7
3	4	0	3	4	0	2	33	red	高速公路	2212.137	2	LINESTRING (651690.838 3456004.237, 652171.081...	2212.136803	2212.136803	3	4
	37	0	3	37	0	35	32	red	高速公路	5939.964	35	LINESTRING (648142.369 3460767.801, 651690.838...	5939.963801	5939.963801	3	37

In [20]:

edges.set_index(['u','v','key'],drop=True,inplace=True)
edges.head()

Out[20]:

			osmid	FID_road21	type	name	lengthQGIS	ID	geometry	mm_len	length	node_start	node_end
u	v	key
0	1	0	0	13	lblue	干线公路	746.941	0	LINESTRING (683092.608 3454982.000, 682886.751...	746.940674	746.940674	0	1
	2	0	1	12	lblue	干线公路	684.967	1	LINESTRING (683649.578 3455380.702, 683092.608...	684.966555	684.966555	0	2
2	7	0	4	11	lblue	干线公路	1826.529	4	LINESTRING (684644.129 3456912.720, 683649.578...	1826.528798	1826.528798	2	7
3	4	0	2	33	red	高速公路	2212.137	2	LINESTRING (651690.838 3456004.237, 652171.081...	2212.136803	2212.136803	3	4
	37	0	35	32	red	高速公路	5939.964	35	LINESTRING (648142.369 3460767.801, 651690.838...	5939.963801	5939.963801	3	37

In [26]:

# 展示某一列大于另外一列的所有行
tmp=edges.loc[edges['node_start']>edges['node_end']]
tmp

Out[26]:

			osmid	FID_road21	type	name	lengthQGIS	ID	geometry	mm_len	length	node_start	node_end
u	v	key

In [54]:

# 构建节点的street_count列
t1=pd.DataFrame(edges['node_start'].value_counts())
t2=pd.DataFrame(edges['node_end'].value_counts())
t1.columns=['street_count']
t2.columns=['street_count']

In [57]:

tt=t1.add(t2,fill_value=0).astype(int)
tt

Out[57]:

	street_count
0	2
1	1
2	2
3	2
4	1
...	...
2565	1
2566	1
2567	1
2568	1
2569	1

2570 rows × 1 columns

In [69]:

len(tt)==len(nodes)

Out[69]:

True

In [14]:

nodes.set_index(['osmid'],drop=True,inplace=True)
nodes.head()

Out[14]:

	y	x	nodeID	geometry
osmid
0	3.454982e+06	683092.608372	0	POINT (683092.608 3454982.000)
1	3.454264e+06	682886.751173	1	POINT (682886.751 3454263.987)
2	3.455381e+06	683649.578442	2	POINT (683649.578 3455380.702)
3	3.456004e+06	651690.837729	3	POINT (651690.838 3456004.237)
4	3.453845e+06	652171.081126	4	POINT (652171.081 3453844.858)

In [72]:

nodesnew=pd.merge(nodes,tt,how='left',left_index=True,right_index=True)
nodesnew.head(10)

Out[72]:

	y	x	nodeID	geometry	street_count
osmid
0	3.454982e+06	683092.608372	0	POINT (683092.608 3454982.000)	2
1	3.454264e+06	682886.751173	1	POINT (682886.751 3454263.987)	1
2	3.455381e+06	683649.578442	2	POINT (683649.578 3455380.702)	2
3	3.456004e+06	651690.837729	3	POINT (651690.838 3456004.237)	2
4	3.453845e+06	652171.081126	4	POINT (652171.081 3453844.858)	1
5	3.456735e+06	687071.291070	5	POINT (687071.291 3456735.217)	2
6	3.454341e+06	686371.578992	6	POINT (686371.579 3454341.316)	1
7	3.456913e+06	684644.128657	7	POINT (684644.129 3456912.720)	2
8	3.457126e+06	699989.079878	8	POINT (699989.080 3457125.993)	2
9	3.454558e+06	703019.672811	9	POINT (703019.673 3454558.207)	1

In [73]:

tt.head(10)

Out[73]:

	street_count
0	2
1	1
2	2
3	2
4	1
5	2
6	1
7	2
8	2
9	1

In [74]:

G = ox.graph_from_gdfs(nodesnew, edges)
fig, ax = ox.plot_graph(G)

In [16]:

G.graph

Out[16]:

{'crs': <Projected CRS: EPSG:32650>
 Name: WGS 84 / UTM zone 50N
 Axis Info [cartesian]:
 - E[east]: Easting (metre)
 - N[north]: Northing (metre)
 Area of Use:
 - name: Between 114°E and 120°E, northern hemisphere between equator and 84°N, onshore and offshore. Brunei. China. Hong Kong. Indonesia. Malaysia - East Malaysia - Sarawak. Mongolia. Philippines. Russian Federation. Taiwan.
 - bounds: (114.0, 0.0, 120.0, 84.0)
 Coordinate Operation:
 - name: UTM zone 50N
 - method: Transverse Mercator
 Datum: World Geodetic System 1984
 - Ellipsoid: WGS 84
 - Prime Meridian: Greenwich}

In [75]:

area_m=nodes.unary_union.convex_hull.area
ox.basic_stats(G,area=area_m,clean_int_tol=15)

Out[75]:

{'n': 2570,
 'm': 3568,
 'k_avg': 2.7766536964980544,
 'edge_length_total': 6067344.658434196,
 'edge_length_avg': 1700.4889737764004,
 'streets_per_node_avg': 2.7766536964980544,
 'streets_per_node_counts': {0: 0,
  1: 78,
  2: 1294,
  3: 360,
  4: 804,
  5: 30,
  6: 4},
 'streets_per_node_proportions': {0: 0.0,
  1: 0.030350194552529183,
  2: 0.5035019455252918,
  3: 0.14007782101167315,
  4: 0.31284046692607004,
  5: 0.011673151750972763,
  6: 0.0015564202334630351},
 'intersection_count': 2492,
 'street_length_total': 6067344.658434196,
 'street_segment_count': 3568,
 'street_length_avg': 1700.4889737764004,
 'circuity_avg': 0.9999999999999983,
 'self_loop_proportion': 0.0,
 'clean_intersection_count': 2492,
 'node_density_km': 0.13292351481358414,
 'intersection_density_km': 0.12888926027838588,
 'edge_density_km': 313.8104192935792,
 'street_density_km': 313.8104192935792,
 'clean_intersection_density_km': 0.12888926027838588}

In [4]:

# 下载数据的时候的代码，因为南京太大了，熟悉OSMnx的时候，就用了秦淮区的一个范围 fig, ax = ox.plot_graph(G)

下载数据的时候找的保存为shapefile数据的代码

In [2]:

# 下载数据的时候的代码
ox.settings.all_oneway = True
ox.settings.log_console = True
G2 = ox.graph_from_place('Nanjing, China', network_type='all',simplify=True)

In [3]:

Gfile = "./nanjing2022"
ox.save_graphml(G2, Gfile)
# 相对应的加载就是 load_graphml
#G = ox.load_graphml(Gfile)

In [12]:

import os
path = "./"
filep_nodes = os.path.join(path, "nodes.shp")
filep_edges = os.path.join(path, "edges.shp")


def stringify_nonnumeric_cols(gdf):  # 转换数据类型
    for col in (c for c in gdf.columns if not c == "geometry"):  # 除 geometry 外的所有列
        if not pd.api.types.is_numeric_dtype(gdf[col]):  # 如果不是数字类型
            gdf[col] = gdf[col].fillna("").astype(str)  # 全部转换为字符串
    return gdf


gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
gdf_nodes = stringify_nonnumeric_cols(gdf_nodes)
gdf_edges = stringify_nonnumeric_cols(gdf_edges)

gdf_edges["fid"] = np.arange(0, gdf_edges.shape[0], dtype='int')  # 增加 fid 列
gdf_nodes.to_file(filep_nodes, encoding="utf-8")
gdf_edges.to_file(filep_edges, encoding="utf-8")

D:\LenovoSoftstore\Install\anaconda\envs\osmnx_env\lib\site-packages\ipykernel_launcher.py:19: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.

获取到的数据可以转换为networkx里的multigraph,digraph这些，和networkx交互； 也可以转换为geopandas里的格式

In [6]:

# you can convert your graph to node and edge GeoPandas GeoDataFrames
gdf_nodes, gdf_edges = ox.graph_to_gdfs(graph)
gdf_nodes.head()

In [10]:

gdf_edges.head()

Out[10]:

			osmid	oneway	name	highway	length	geometry	bridge	lanes	maxspeed	ref	tunnel	access	width	junction	service
u	v	key
304815373	7232707226	0	105242650	True	内环东线	trunk	7.907	LINESTRING (118.79551 32.02452, 118.79553 32.0...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	3999186863	0	61822454	True	象房村路	residential	160.252	LINESTRING (118.79551 32.02452, 118.79688 32.0...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
304815377	770095460	0	[254471398, 254471399]	True	龙蟠南路	primary	267.154	LINESTRING (118.79258 32.01546, 118.79208 32.0...	yes	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	2602977729	0	835800514	True	NaN	trunk_link	361.196	LINESTRING (118.79258 32.01546, 118.79234 32.0...	yes	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
305194317	3302177260	0	142491926	True	中华路	secondary	59.662	LINESTRING (118.77634 32.01394, 118.77694 32.0...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

⚡ 可以采用的一个措施：
将节点和边那种的shp导入进来，通过geopandas类似的进行处理，形成geodataframes类似的格式，然后在转换成multidigraph； 也可以往上面两个结构靠，就是保证gdf_nodes是用omsid编码的，gdf_edges是用u,v,key来多重索引的；

You can create a graph from node/edge GeoDataFrames, as long as gdf_nodes is indexed by osmid and gdf_edges is multi-indexed by u, v, key (following normal MultiDiGraph structure). This allows you to load graph node/edge shapefiles or GeoPackage layers as GeoDataFrames then convert to a MultiDiGraph for graph analytics

In [11]:

G[304815377][770095460]  # 其实就是nexworkx里数据访问的方式

Out[11]:

AtlasView({0: {'osmid': [254471398, 254471399], 'bridge': 'yes', 'oneway': True, 'name': '龙蟠南路', 'highway': 'primary', 'length': 267.154, 'geometry': <shapely.geometry.linestring.LineString object at 0x000002EF6DA94248>}})

In [5]:

# convert node/edge GeoPandas GeoDataFrames to a NetworkX MultiDiGraph
G2 = ox.graph_from_gdfs(gdf_nodes, gdf_edges, graph_attrs=G.graph)

In [13]:

G.graph

Out[13]:

{'created_date': '2023-01-20 09:20:55',
 'created_with': 'OSMnx 1.0.1',
 'crs': 'epsg:4326',
 'simplified': True}

In [6]:

# If to_crs is None, project the graph to the UTM CRS for the UTM zone in which the graph’s centroid lies. Otherwise, project the graph to the CRS defined by to_crs
G_proj = ox.project_graph(G)
G_proj

Out[6]:

<networkx.classes.multidigraph.MultiDiGraph at 0x1b5fbe35908>

In [7]:

nodes_proj = ox.graph_to_gdfs(G_proj, edges=False)
nodes_proj  # 已经实现了投影，可以看x和y

Out[7]:

	y	x	street_count	lon	lat	highway	ref	geometry
osmid
304815373	3.544563e+06	669563.070442	3	118.795515	32.024524	NaN	NaN	POINT (669563.070 3544563.165)
7232707226	3.544571e+06	669564.566721	3	118.795532	32.024594	NaN	NaN	POINT (669564.567 3544570.909)
3999186863	3.544563e+06	669723.638818	3	118.797215	32.024502	NaN	NaN	POINT (669723.639 3544563.429)
304815377	3.543553e+06	669302.813527	3	118.792583	32.015456	NaN	NaN	POINT (669302.814 3543553.281)
770095460	3.543299e+06	669225.275707	3	118.791718	32.013173	NaN	NaN	POINT (669225.276 3543298.761)
...	...	...	...	...	...	...	...	...
10303904934	3.543303e+06	668114.891707	3	118.779967	32.013375	NaN	NaN	POINT (668114.892 3543302.843)
10303904966	3.543465e+06	668468.251299	1	118.783735	32.014786	NaN	NaN	POINT (668468.251 3543465.101)
10303904989	3.543423e+06	668458.820824	3	118.783628	32.014406	NaN	NaN	POINT (668458.821 3543422.794)
10303904974	3.543432e+06	668171.975837	1	118.780593	32.014533	NaN	NaN	POINT (668171.976 3543432.157)
10303904990	3.543433e+06	668418.461737	1	118.783202	32.014502	NaN	NaN	POINT (668418.462 3543432.762)

1929 rows × 8 columns

In [6]:

G_proj.graph['crs']

Out[6]:

<Projected CRS: +proj=utm +zone=50 +ellps=WGS84 +datum=WGS84 +unit ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: UTM zone 50N
- method: Transverse Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [8]:

graph_area_m = nodes_proj.unary_union.convex_hull.area
graph_area_m

Out[8]:

54649971.6096648

In [23]:

nodes_proj.unary_union  # 包含所有地理实体的集合体，是geopandas包里的

Out[23]:

In [24]:

nodes_proj.unary_union.convex_hull  # 对应的凸包

Out[24]:

In [10]:

# show some basic stats about the network
ox.basic_stats(G_proj, area=graph_area_m, clean_int_tol=15)  # 如果提供面积，可以计算density指标，如果提供clean_int_tol可以计算clean intersection

Out[10]:

{'n': 1929,
 'm': 4945,
 'k_avg': 5.127008812856403,
 'edge_length_total': 774358.6850000017,
 'edge_length_avg': 156.59427401415604,
 'streets_per_node_avg': 3.0031104199066876,
 'streets_per_node_counts': {0: 0, 1: 223, 2: 12, 3: 1246, 4: 432, 5: 16},
 'streets_per_node_proportions': {0: 0.0,
  1: 0.11560393986521514,
  2: 0.006220839813374806,
  3: 0.6459305339554173,
  4: 0.223950233281493,
  5: 0.00829445308449974},
 'intersection_count': 1706,
 'street_length_total': 486759.9220000005,
 'street_segment_count': 2836,
 'street_length_avg': 171.636079689704,
 'circuity_avg': 1.0441095736027672,
 'self_loop_proportion': 0.0007052186177715092,
 'clean_intersection_count': 1273,
 'node_density_km': 35.29736508882025,
 'intersection_density_km': 31.21685061769173,
 'edge_density_km': 14169.425201733446,
 'street_density_km': 8906.86504792104,
 'clean_intersection_density_km': 23.29369920065743}

'n': 1929,节点的数量
'm': 4945, 边的数量
'k_avg': 5.127008812856403, 平均节点度（包括了入度和出度） 'intersection_count': 1706, 交叉点数量
'streets_per_node_avg': 3.0031104199066876,无向图中，每个节点连接的街道数目（边）的平均值
'streets_per_node_counts': {0: 0, 1: 223, 2: 12, 3: 1246, 4: 432, 5: 16},无向图中，连接到一个节点的街道数目：这个数字下的街道数目一共有多少
'streets_per_node_proportion': {0: 0.0, 1: 0.11560393986521514, 2: 0.006220839813374806, 3: 0.6459305339554173, 4: 0.223950233281493, 5: 0.00829445308449974},无向图中，对应的改成了占比
'edge_length_total': 774358.6850000017,米制的所有边长度的总和
'edge_length_avg': 156.59427401415604,上面的均值
'street_length_total': 486759.9220000005,在无向图中，米制的所有边长度的总和
'street_length_avg': 171.636079689704,上面的均值
'street_segments_count': 2836,无向图中，边的数量
'node_density_km': 35.29736508882025,平方千米下的，面积除以n
'intersection_density_km': 31.21685061769173,平方千米下的，面积除以intersection_count
'edge_density_km': 14169.425201733446,km下的，edge_length_total divided by area
'street_density_km': 8906.86504792104,km下的，street_length_total divided by area
'circuity_avg': 2.0995288875627176e-05,每个边的节点之间最大圆环距离，除以，边的总长度
'self_loop_proportion': 0.0008088978766430738,表示proportion of edges that have a single node as its endpoints (ie, the edge links nodes u and v, and u==v)
'clean_intersection_count': 8906.86504792104,表示number of intersections in street network, merging complex ones into single points
'clean_intersection_density_km': 23.29369920065743,km下的，clean_intersection_count divided by area

下面配的图针对的就是实验的2021年的¶

In [78]:

# unpack dicts into individiual keys:values
stats = ox.basic_stats(G, area=area_m,clean_int_tol=15)
for k, count in stats["streets_per_node_counts"].items():
    stats["{}way_int_count".format(k)] = count
for k, proportion in stats["streets_per_node_proportions"].items():
    stats["{}way_int_prop".format(k)] = proportion

# delete the no longer needed dict elements
del stats["streets_per_node_counts"]
del stats["streets_per_node_proportions"]

# load as a pandas dataframe
stats2021=pd.DataFrame(pd.Series(stats, name="value")).round(3)

In [79]:

stats2021

Out[79]:

	value
n	2570.000
m	3568.000
k_avg	2.777
edge_length_total	6067344.658
edge_length_avg	1700.489
streets_per_node_avg	2.777
intersection_count	2492.000
street_length_total	6067344.658
street_segment_count	3568.000
street_length_avg	1700.489
circuity_avg	1.000
self_loop_proportion	0.000
clean_intersection_count	2492.000
node_density_km	0.133
intersection_density_km	0.129
edge_density_km	313.810
street_density_km	313.810
clean_intersection_density_km	0.129
0way_int_count	0.000
1way_int_count	78.000
2way_int_count	1294.000
3way_int_count	360.000
4way_int_count	804.000
5way_int_count	30.000
6way_int_count	4.000
0way_int_prop	0.000
1way_int_prop	0.030
2way_int_prop	0.504
3way_int_prop	0.140
4way_int_prop	0.313
5way_int_prop	0.012
6way_int_prop	0.002

In [11]:

# 作为geopackage或者graphml把graph存入磁盘
ox.save_graph_geopackage(G,filepath="mynetwork.gpkg")
ox.save_graphml(G,filepath="mynetwork.graphml")

⚡可视化街道中心性，就是利用relative closeness centrality来进行可视化,这里能计算closeness_centrality,那也能计算其他参数，那么这些参数的计算公式，和空间句法里计算参数的区别，这些东西都可以进一步学习

In [80]:

edge_centrality=nx.closeness_centrality(nx.line_graph(G))
nx.set_edge_attributes(G,edge_centrality,"edge_centrality")
ec=ox.plot.get_edge_colors_by_attr(G,"edge_centrality",cmap="inferno")
fig,ax=ox.plot_graph(G,edge_color=ec,edge_linewidth=2,node_size=0)

In [18]:

# 计算双向图（没有平行边）的betweenness
bc=nx.betweenness_centrality(ox.get_digraph(G),weight="length")
max_node,max_bc=max(bc.items(),key=lambda x:x[1]) # 是按照第1个（从0开始）的大小来选的最大的item
max_node,max_bc #比例是44%，说明有这个节点有~44%的最短路径是通过它的

Out[18]:

(2262867596, 0.44300260478332937)

In [19]:

bc.items()

Out[19]:

dict_items([(41542771, 0.01941747572815534), (41542772, 0.0013260715131423158), (41542778, 0.0), (41544344, 0.16123608808903622), (41544346, 0.06331991475254559), (41577408, 0.04525219038598153), (41597889, 0.007080274686242008), (41597891, 0.0), (41612013, 0.020577788302154867), (41612044, 0.053563817191569975), (41612048, 0.05441629173573289), (41612050, 0.00970873786407767), (41615126, 0.4426710869050438), (41615128, 0.03840871418422922), (41659040, 0.020672507695950745), (41672313, 0.050532796590101826), (41672314, 0.10348093772199858), (41672318, 0.23720104191333174), (41672321, 0.03805351645749467), (41717386, 0.0026521430262846316), (41756544, 0.15027231825716317), (41756546, 0.03840871418422922), (41794465, 0.011674165285342174), (41794468, 0.0), (41849849, 0.010253374378403977), (41849851, 0.006180440445181151), (41849858, 0.01948851527350225), (41864246, 0.004735969689793985), (41941829, 0.0), (41979079, 0.04593890599100166), (42238826, 0.028913094956192282), (42397551, 0.047643855079327495), (42449925, 0.0), (42452350, 0.11200568316362776), (42500019, 0.027918541321335543), (42666429, 0.0), (42708110, 0.0), (371833008, 0.14956192280369407), (371839652, 0.10982713710632253), (371839813, 0.21498934406819797), (371848315, 0.23594600994553636), (371852313, 0.054534690977977744), (2262867586, 0.0), (2262867591, 0.08555529244612835), (2262867596, 0.44300260478332937), (2262867612, 0.00961401847028179), (2262867613, 0.07624911200568317), (2262867619, 0.00961401847028179), (2262867623, 4.735969689793985e-05), (2262867639, 4.735969689793985e-05), (2262867640, 0.34004262372720817), (2262867643, 0.00970873786407767), (2262867644, 0.0), (3545571404, 0.0), (3640359571, 0.0), (3640359581, 0.08190859578498698), (3640359590, 0.02903149419843713), (3851123454, 0.004617570447549135), (3851123455, 0.005138527113426474), (3851123456, 4.735969689793985e-05), (4848535844, 0.0191333175467677), (4848535849, 0.00966137816717973), (4848535860, 0.0), (4848535864, 0.00970873786407767), (4848557987, 0.25884442339569025), (4848557990, 0.0), (4848557991, 0.03443049964480227), (4848557992, 0.08756807956429079), (4848557994, 0.04840161022969453), (4848558000, 0.07705422685294815), (4848558134, 0.08974662562159602), (4848558135, 0.0), (4848558138, 0.010466493014444708), (4848558153, 0.09372484016102298), (4848558155, 0.02178546057305233), (4848558156, 0.0852000947193938), (4848558222, 0.01932275633435946), (4848558238, 0.02884205541084537), (4848558243, 0.0030310206014681506), (4848558245, 0.005588444233956903), (4848558330, 0.01690741179256453), (4848558331, 0.006819796353303339), (4848558347, 0.004546530902202226), (4848558353, 0.06526166232536111), (4848558355, 0.04584418659720578), (4848558356, 0.01922803694056358), (4848558360, 0.027752782382192753), (4848558363, 0.08827847501775989), (4848558379, 0.01847028179019654), (4848558388, 0.027184466019417475), (4848558405, 0.024153445417949324), (4848558407, 0.05574236324887521), (4848558409, 0.052569263556713236), (4850315857, 0.11835188254795169), (4850315858, 0.023679848448969927), (4850315859, 0.0), (4850315861, 0.00970873786407767), (4850315873, 0.17830925882074355), (4850315874, 0.1939379587970637), (4850315875, 0.3254321572341937), (5052915276, 0.04825953113900071), (5052915279, 0.00966137816717973), (5052915280, 0.0), (5052915281, 0.0), (5701320606, 0.057518351882547954), (5701320609, 0.062396400663035756), (5701320612, 0.39985792090930616), (6058659222, 0.12848685768411083), (6058659224, 0.00028415818138763913), (6058659225, 0.050343357802510064), (6058659226, 0.00970873786407767), (6058659227, 0.00966137816717973), (6058659228, 0.1499171205304286), (6902845091, 0.0), (6902845092, 0.01932275633435946), (6902879136, 0.0191333175467677), (6902879139, 0.05607388112716079), (6902879140, 0.1775515036703765), (6902879143, 0.14326308311626806), (6902879145, 0.4123608808903623), (6902879146, 0.0096850580156287), (6902879148, 0.004854368932038835), (6902879151, 0.0096850580156287), (6902879152, 0.38233483305706845), (7045852575, 0.08420554108453705), (7045852579, 0.00970873786407767), (7045852581, 0.0), (7045856798, 0.01932275633435946), (7046571485, 0.00014207909069381957), (7046571497, 0.047643855079327495), (7048949334, 0.11565237982476913), (7048949336, 0.0382666350935354), (7048949340, 0.07501775988633673), (7048949342, 0.00970873786407767), (7048949343, 0.06644565474780961), (7130067150, 0.17210513852711343), (7130067152, 0.0), (7655326298, 0.004854368932038835), (7655335686, 0.23059436419606916), (7655335691, 0.33817191569973953), (7655352447, 0.037130002367984846), (7655352462, 0.0), (7737996633, 0.01922803694056358), (7737996634, 0.0), (7737996635, 0.0), (7737996668, 4.735969689793985e-05), (7737996669, 0.0), (7738026875, 0.0), (7738026884, 0.0), (7738026886, 0.0193701160312574), (7738026888, 0.0), (7738026893, 0.11702581103480937), (7738026908, 0.3151314231588918), (7738026909, 0.0), (7738026910, 0.15216670613308075), (7738026912, 0.3024390243902439), (7738026923, 0.06649301444470755), (7738026926, 0.0), (7738026928, 0.12071986739284869), (7738026929, 0.0), (7738026931, 0.0193701160312574), (7738026932, 0.0), (7738026933, 0.03840871418422922), (7738026934, 0.0), (7738026957, 0.0), (7738026958, 0.00970873786407767), (7738026983, 0.046317783566185176), (7738026984, 0.0), (7738026989, 0.0), (7738026990, 0.13634856736916884), (7738026994, 0.0), (7738027004, 0.01929907648591049), (7738027005, 0.04458915462941037), (7738027006, 0.01929907648591049), (7738027009, 0.05332701870708027), (7738027028, 0.0), (7738027032, 0.15964953824295525), (7738027067, 0.09793985318493961), (7738027068, 0.06852948141131897), (7738027074, 0.07103954534690977), (7738027095, 0.0), (7738027116, 0.0193701160312574), (7738027120, 0.0), (7738027121, 0.2259057542031731), (7738027122, 0.0), (7738027128, 0.0), (7738027130, 0.09727681742836845), (7738027131, 0.0), (7738027135, 0.11915699739521667), (7738027137, 0.15415581340279422), (7738027141, 0.0), (7738027164, 0.00970873786407767), (7738027165, 0.00966137816717973), (7738027166, 0.0), (7738027167, 4.735969689793985e-05), (7738027184, 0.00956665877338385), (7738027187, 9.47193937958797e-05), (9527190857, 0.20170494908832584), (9939284239, 0.06410134975136159), (9939284242, 0.047004499171205305), (9939284243, 0.0), (9939284244, 0.0004262372720814587), (10263950598, 0.024437603599336964), (10263950599, 0.019725313757991948), (10263950600, 0.01982003315178783), (10288549305, 0.11139000710395454), (10559494634, 0.01932275633435946)])

In [20]:

nc = ["r" if node == max_node else "w" for node in G.nodes]
ns = [80 if node == max_node else 15 for node in G.nodes]
fig, ax = ox.plot_graph(G, node_size=ns, node_color=nc, node_zorder=2)

In [21]:

# add the betweenness centraliy values as new node attributes, then plot
nx.set_node_attributes(G, bc, "bc")
nc = ox.plot.get_node_colors_by_attr(G, "bc", cmap="plasma")
fig, ax = ox.plot_graph(
    G,
    node_color=nc,
    node_size=30,
    node_zorder=2,
    edge_linewidth=0.2,
    edge_color="w",
)

routing

In [13]:

G=ox.speed.add_edge_speeds(G)
G=ox.speed.add_edge_travel_times(G)

In [16]:

orig=ox.distance.nearest_nodes(G,Y=32.0306,X=118.7854)
dest=ox.distance.nearest_nodes(G,Y=32.0367,X=118.8082)
route=ox.shortest_path(G,orig,dest,weight="travel_time")
fig,ax=ox.plot_graph_route(G,route,node_size=0)

In [17]:

# 求这条最短路径的长度
edge_lengths=ox.utils_graph.get_route_edge_attributes(G,route,"length")# 因为这个route的标注，使得只有在route上的才有length值，所以下面可以累加再round
round(sum(edge_lengths))

Out[17]:

2926

In [18]:

# 环形距离
orig_x=G.nodes[orig]["x"]
orig_y=G.nodes[orig]["y"]
dest_x=G.nodes[dest]["x"]
dest_y=G.nodes[dest]["y"]

In [20]:

dest

Out[20]:

1038950909

In [21]:

round(ox.distance.great_circle_vec(orig_y,orig_x,dest_y,dest_y))

Out[21]:

7918205

利用google API或者是遥感影像文件，获取得到路网节点对应的高程数据，再赋予色彩进行可视化
try:

except ImportError:

获取网络的其他方式

In [ ]:

place={"city":"Nanjing","state":"Qixia","country":"China"}
G=ox.graph_from_place(place,network_type="drive",truncate_by_edge=True)
fig,ax=ox.plot_graph(G,figsize=(10,10),node_size=0,edge_color="y",edge_linewidth=0.2)

In [ ]:

G=ox.graph_from_place("Qixia, Nanjing, China", network_type="all")
fig,ax=ox.plot_graph(G,node_size=0,edge_linewidth=0.5)

In [ ]:

wurster_hall=(37,112)
one_mile=1609 # menter
G=ox.graph_from_point(wurster_hall,dist=one_mile,network_type="drive")
fig,ax=ox.plot_graph(G,node_size=0)

In [14]:

# 获取网络上的其他基础设施数据
G=ox.graph_from_place("Qinhuai, Nanjing, China",retain_all=False,truncate_by_edge=True,simplify=True,custom_filter='["railway"~"subway"]')
fig,ax=ox.plot_graph(G,node_size=0,edge_color="w",edge_linewidth=0.2)

In [ ]:

# 获取其他数据，包括建筑，是通过设置ox.geometries_from_place里的参数来实现的

查询，简化，可视化，保存¶

In [ ]:

# 获取边界
ox.geocode_to_gdf
ox.project_gdf
().plot()

# 获取多个边界
place_names={"","",""}
ox.geocode_to_gdf
().to_file("",driver="GPKG")
ox.project_gdf
().plot()

# 如果知道地点对应的OSM ID，可以用这个来查询
ox.geocode_to_gdf(["",""],by_osmid=True)

In [28]:

# 对街道网络进行下载和建模
north, south, east, west=37.1,37,112,112.1
ox.graph_from_bbox

ox.graph_from_point(center_point,dist=500,) # 500meter, 此外这五百米以内的道路，是不管单向还是双向，都包括的

ox.graph_from_address

ox.graph_from_place #包括了一个和多个地址，从这些地址获取到网络数据

ox.save_graph_geopackage

# 输入一个矩形，从这个里面下载
ox.graph_from_polygon

# 从.osm xml数据中获取
ox.graph_from_xml

⚡：简化的效果挺好，思考下能不能把2011或者2021的数据套进来

In [33]:

# 简化
location_point = (33.299896, -111.831638)
G = ox.graph_from_point(location_point, network_type="drive_service", dist=500, simplify=False)

In [34]:

nc=["r" if ox.simplification._is_endpoint(G,node) else "y" for node in G.nodes()]
fig,ax=ox.plot_graph(G,node_color=nc)

In [35]:

G=ox.simplify_graph(G)
fig,ax=ox.plot_graph(G,node_color="r")

In [36]:

# 高亮所有multiline
ec=["gray" if k==0 or u==v else "r" for u,v,k in G.edges(keys=True)]
fig,ax=ox.plot_graph(G,node_color="w",node_edgecolor="k",node_size=50,edge_color=ec,edge_linewidth=3)

In [39]:

# 高亮所有单向边
ec = ["r" if data["oneway"] else "w" for u, v, key, data in G2.edges(keys=True, data=True)]
fig, ax = ox.plot_graph(G2, node_size=0, edge_color=ec, edge_linewidth=1.5, edge_alpha=0.7)

In [41]:

stats = ox.basic_stats(G)
stats["circuity_avg"]

Out[41]:

1.1358685392173113

routing, speed imputation and travel times¶

In [42]:

np.random.seed(0)

In [43]:

Gp=ox.project_graph(G)

In [44]:

# 获取最近的节点或者边
# 随机采样n个点
points=ox.utils_geo.sample_points(ox.get_undirected(Gp),n=100)
X=points.x.values
Y=points.y.values

In [46]:

X0=X.mean()
Y0=Y.mean()

In [47]:

# 获取距离一组点最近的点，以及距离
nodes, dists=ox.nearest_nodes(Gp,X,Y,return_dist=True)

In [48]:

dists

Out[48]:

[13.02141728799831,
 134.19332912976063,
 73.72831276320169,
 1.6346128584053425,
 31.051578351764224,
 3.5092931187607337,
 87.64878994034457,
 45.24829443936839,
 21.675178504056582,
 1.9023468458797488,
 99.26566791577433,
 48.04498166737192,
 31.751374377487615,
 2.959880615277566,
 20.51153786015448,
 73.97941904897311,
 39.17673517978925,
 51.19340120660089,
 26.292049431254576,
 48.73309618738953,
 19.68170187060116,
 110.61499421077518,
 2.07086297197591,
 79.14645854179912,
 103.63040736078571,
 18.294585058685378,
 3.1758686606072075,
 74.5396610113333,
 89.77840652642139,
 35.60307440223833,
 64.69173210012778,
 51.769136898939884,
 14.959001661483766,
 31.514877405422276,
 80.07542594380753,
 113.91072204216687,
 118.51794337137711,
 96.5571174517652,
 85.98895825266781,
 74.43497754359082,
 38.135585234599276,
 79.2778038957343,
 16.545247470930708,
 11.00946274600901,
 51.210083147405726,
 78.29977045162616,
 18.127620581777,
 3.7319266397329214,
 56.175853888514595,
 0.42502303011911563,
 15.789889605670393,
 29.724675648960957,
 30.74656359807148,
 72.58172722111378,
 22.394805722575228,
 18.393164301227127,
 14.209634163942912,
 61.163073088422514,
 48.11301482510856,
 15.228321340803829,
 36.6686502821069,
 99.101673081213,
 73.4806591657184,
 11.305126802576373,
 2.6208123714609033,
 53.15459623552805,
 1.1616710308975458,
 67.26829825223594,
 33.01123945419404,
 32.40340313009194,
 34.31105211554059,
 6.804884421846613,
 14.326963767036661,
 5.156922180681355,
 102.53175736335805,
 30.47406260073624,
 39.79955632535689,
 18.026335309299743,
 5.634745347636121,
 7.980388514743713,
 13.493124386913825,
 47.25059828852815,
 20.60957777585811,
 41.04494671981354,
 2.617842381746908,
 65.06686540652039,
 102.23000946268759,
 77.42432758636345,
 28.513828742105808,
 77.10283794168458,
 32.625603970107015,
 55.89264138909181,
 15.753421513281532,
 27.64679165312897,
 14.877395834709823,
 54.76616717587069,
 50.14486761670674,
 55.70394320916568,
 7.094620558333998,
 23.985772284941927]

In [50]:

# 也可以是针对一个点，发现最近的点
node=ox.nearest_nodes(Gp,X0,Y0)
node

Out[50]:

371839652

In [ ]:

# 对边同理
ox.nearest_edges

In [52]:

# 用距离作为权重，找两个点之间的最短路径
orig = list(G)[0]
dest = list(G)[12]
route = ox.shortest_path(G, orig, dest, weight="length")
fig, ax = ox.plot_graph_route(G, route, route_color="y", route_linewidth=6, node_size=0)

In [54]:

# 或者是k条最短路径
routes = ox.k_shortest_paths(G, orig, dest, k=3, weight="length")
fig, ax = ox.plot_graph_routes(G, list(routes), route_colors="y", route_linewidth=4, node_size=0)

⚡这一部分可以看看要不要用，可不可以用
Imputes free-flow travel speeds for all edges based on mean maxspeed value of edges, per highway type. This mean-imputation can obviously be imprecise, and the caller can override it by passing in hwy_speeds and/or fallback arguments that correspond to local speed limit standards.

In [57]:

G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)
edges = ox.graph_to_gdfs(G, nodes=False)
edges["highway"] = edges["highway"].astype(str)
edges.groupby("highway")[["length", "speed_kph", "travel_time"]].mean().round(1)

Out[57]:

	length	speed_kph	travel_time
highway
['residential', 'service']	219.8	40.2	19.7
['residential', 'tertiary']	238.2	40.2	21.3
['residential', 'unclassified']	265.9	40.2	23.8
['unclassified', 'residential']	265.9	40.2	23.8
residential	115.0	40.2	10.3
service	67.6	40.2	6.1
tertiary	72.4	40.2	6.5
unclassified	239.2	40.2	21.4

In [58]:

G.nodes(data=True)

Out[58]:

NodeDataView({41542771: {'y': 33.2963878, 'x': -111.8302054, 'street_count': 3}, 41542772: {'y': 33.2963871, 'x': -111.8296124, 'street_count': 3}, 41542778: {'y': 33.2963861, 'x': -111.8284497, 'street_count': 3}, 41544344: {'y': 33.303419, 'x': -111.8328293, 'street_count': 4}, 41577408: {'y': 33.2978321, 'x': -111.829308, 'street_count': 3}, 41597889: {'y': 33.2995713, 'x': -111.8267331, 'street_count': 3}, 41597891: {'y': 33.2995642, 'x': -111.8288445, 'street_count': 3}, 41612013: {'y': 33.2970891, 'x': -111.8284497, 'street_count': 3}, 41612044: {'y': 33.2989257, 'x': -111.8288359, 'street_count': 3}, 41612048: {'y': 33.301034, 'x': -111.8287095, 'street_count': 3}, 41612050: {'y': 33.302161, 'x': -111.8287501, 'street_count': 3}, 41615126: {'y': 33.3013315, 'x': -111.8328592, 'street_count': 3}, 41615128: {'y': 33.3013334, 'x': -111.833994, 'highway': 'turning_circle', 'street_count': 1}, 41659040: {'y': 33.3038903, 'x': -111.8299606, 'street_count': 3}, 41672313: {'y': 33.2989108, 'x': -111.8365753, 'street_count': 3}, 41672314: {'y': 33.2989147, 'x': -111.835408, 'street_count': 3}, 41672318: {'y': 33.2989291, 'x': -111.8303403, 'street_count': 3}, 41672321: {'y': 33.2989318, 'x': -111.8267393, 'street_count': 4}, 41717386: {'y': 33.2970891, 'x': -111.8279519, 'street_count': 3}, 41756544: {'y': 33.302656, 'x': -111.8328442, 'street_count': 3}, 41756546: {'y': 33.3025514, 'x': -111.8340428, 'highway': 'turning_circle', 'street_count': 1}, 41794465: {'y': 33.2984307, 'x': -111.827145, 'street_count': 3}, 41794468: {'y': 33.2979385, 'x': -111.82706, 'highway': 'turning_circle', 'street_count': 1}, 41849849: {'y': 33.3010551, 'x': -111.8267189, 'street_count': 3}, 41849851: {'y': 33.3003102, 'x': -111.8267245, 'street_count': 3}, 41864246: {'y': 33.3039416, 'x': -111.8280492, 'street_count': 3}, 41941829: {'y': 33.29861, 'x': -111.8297457, 'highway': 'turning_circle', 'street_count': 1}, 41979079: {'y': 33.3039357, 'x': -111.8318553, 'street_count': 3}, 42238826: {'y': 33.2970963, 'x': -111.8289732, 'street_count': 3}, 42397551: {'y': 33.3033933, 'x': -111.836549, 'street_count': 3}, 42449935: {'y': 33.2957046, 'x': -111.8266043, 'street_count': 2}, 42452350: {'y': 33.2978352, 'x': -111.830192, 'street_count': 3}, 42500019: {'y': 33.3010504, 'x': -111.8277324, 'street_count': 3}, 42666429: {'y': 33.3021037, 'x': -111.8289389, 'highway': 'turning_circle', 'street_count': 1}, 42708110: {'y': 33.3013576, 'x': -111.8277287, 'highway': 'turning_circle', 'street_count': 1}, 371834252: {'y': 33.296917, 'x': -111.8355479, 'street_count': 2}, 371839652: {'y': 33.2992775, 'x': -111.8304364, 'street_count': 3}, 371839813: {'y': 33.2989311, 'x': -111.8301766, 'street_count': 3}, 371845086: {'y': 33.2965224, 'x': -111.8264498, 'street_count': 2}, 371845424: {'y': 33.2971896, 'x': -111.8267846, 'street_count': 2}, 371848315: {'y': 33.2989196, 'x': -111.8339511, 'street_count': 3}, 2262867640: {'y': 33.3016904, 'x': -111.8328566, 'street_count': 3}, 2262867644: {'y': 33.3017162, 'x': -111.8298663, 'street_count': 1}, 3372780078: {'y': 33.2989113, 'x': -111.836983, 'street_count': 2}, 3372780204: {'y': 33.3033917, 'x': -111.837007, 'street_count': 2}, 3545571404: {'y': 33.3032242, 'x': -111.8346051, 'street_count': 1}, 3640359590: {'y': 33.3034511, 'x': -111.829964, 'street_count': 3}, 3640364793: {'y': 33.3035253, 'x': -111.8271125, 'street_count': 2}, 3823450766: {'y': 33.2989347, 'x': -111.8263069, 'street_count': 2}, 3851123455: {'y': 33.3043135, 'x': -111.8299519, 'street_count': 4}, 4848535844: {'y': 33.2956364, 'x': -111.8340035, 'street_count': 3}, 4848535845: {'y': 33.295633, 'x': -111.8338967, 'street_count': 1}, 4848535848: {'y': 33.2954493, 'x': -111.8338983, 'street_count': 1}, 4848535849: {'y': 33.2954492, 'x': -111.8340067, 'street_count': 3}, 4848557987: {'y': 33.2989237, 'x': -111.832886, 'highway': 'stop', 'street_count': 3}, 4848557992: {'y': 33.2978655, 'x': -111.8339717, 'street_count': 3}, 4848557993: {'y': 33.2978666, 'x': -111.8338489, 'street_count': 1}, 4848557999: {'y': 33.2986504, 'x': -111.8338515, 'street_count': 1}, 4848558000: {'y': 33.2986537, 'x': -111.8339561, 'street_count': 3}, 4850315862: {'y': 33.2985951, 'x': -111.8320715, 'street_count': 1}, 4850315869: {'y': 33.2980148, 'x': -111.8320688, 'street_count': 1}, 4850315874: {'y': 33.2986516, 'x': -111.8318473, 'street_count': 3}, 4850315875: {'y': 33.2989193, 'x': -111.831842, 'street_count': 3}, 5701320609: {'y': 33.2993325, 'x': -111.831753, 'street_count': 3}, 5701320612: {'y': 33.2993428, 'x': -111.8328754, 'street_count': 3}, 6058659222: {'y': 33.2980071, 'x': -111.8301879, 'street_count': 3}, 6058659224: {'y': 33.2973234, 'x': -111.8303825, 'street_count': 1}, 6058659225: {'y': 33.29732, 'x': -111.830192, 'street_count': 3}, 6058659228: {'y': 33.2985058, 'x': -111.8301761, 'street_count': 3}, 7046571485: {'y': 33.2954046, 'x': -111.834007, 'street_count': 2}, 7655352447: {'y': 33.2964013, 'x': -111.8302054, 'street_count': 3}, 7655352448: {'y': 33.2963997, 'x': -111.8303308, 'street_count': 1}, 7738027128: {'y': 33.2998071, 'x': -111.8349029, 'street_count': 1}, 7738027130: {'y': 33.2989146, 'x': -111.835433, 'street_count': 3}, 9527190857: {'y': 33.2989248, 'x': -111.8305289, 'street_count': 3}, 9939284239: {'y': 33.3034439, 'x': -111.8315241, 'street_count': 3}, 9939284242: {'y': 33.3034474, 'x': -111.8307561, 'street_count': 3}, 9939284243: {'y': 33.3032476, 'x': -111.8315291, 'street_count': 3}, 9939284244: {'y': 33.303254, 'x': -111.8307597, 'street_count': 3}, 10263950590: {'y': 33.3039436, 'x': -111.8267341, 'street_count': 2}, 10263950598: {'y': 33.3039351, 'x': -111.8300117, 'street_count': 3}, 10263950599: {'y': 33.3039376, 'x': -111.8299024, 'street_count': 3}, 10263950600: {'y': 33.3039816, 'x': -111.8299574, 'street_count': 3}, 10288549305: {'y': 33.3034148, 'x': -111.8346065, 'street_count': 3}})

In [59]:

# 这里自己填充相应的数值 (km/hour)
# to fill in edges with missing `maxspeed` from OSM
hwy_speeds = {"residential": 35, "secondary": 50, "tertiary": 60}
G = ox.add_edge_speeds(G, hwy_speeds)
G = ox.add_edge_travel_times(G)

In [61]:

# calculate two routes by minimizing travel distance vs travel time
orig = list(G)[1]
dest = list(G)[12]
route1 = ox.shortest_path(G, orig, dest, weight="length")
route2 = ox.shortest_path(G, orig, dest, weight="travel_time")

In [68]:

# plot the routes
fig, ax = ox.plot_graph_routes(
    G, routes=[route1, route2], route_colors=["r", "y"], route_linewidth=6, node_size=0
)

In [69]:

# compare the two routes
route1_length = int(sum(ox.utils_graph.get_route_edge_attributes(G, route1, "length")))
route2_length = int(sum(ox.utils_graph.get_route_edge_attributes(G, route2, "length")))
route1_time = int(sum(ox.utils_graph.get_route_edge_attributes(G, route1, "travel_time")))
route2_time = int(sum(ox.utils_graph.get_route_edge_attributes(G, route2, "travel_time")))
print("Route 1 is", route1_length, "meters and takes", route1_time, "seconds.")
print("Route 2 is", route2_length, "meters and takes", route2_time, "seconds.")

Route 1 is 962 meters and takes 73 seconds.
Route 2 is 962 meters and takes 73 seconds.

In [ ]:

# multiprocessing

# calculate 100,000 shortest-path routes using random origin-destination pairs
n = 100000
origs = np.random.choice(G.nodes, size=n, replace=True)
dests = np.random.choice(G.nodes, size=n, replace=True)

%%time
# it takes 3.2 seconds to solve all the routes using all the cores on my computer
# I have a 24-thread AMD 5900x: performance will depend on your specific CPU
routes = ox.shortest_path(G, origs, dests, weight="travel_time", cpus=None)

In [ ]:

%%time
# it takes 43 seconds to solve all the routes using just 1 core on my computer
routes = ox.shortest_path(G, origs, dests, weight="travel_time", cpus=1)

In [ ]:

# how many total results did we get
print(len(routes))

# and how many were solvable paths
# some will be unsolvable due to directed graph perimeter effects
routes_valid = [r for r in routes if r is not None]
print(len(routes_valid))

在导航的时候，正确处理了单项的街道；此外，在节点之间如果有平行的边，OSMNX会选择最abs的边来绘制

graph place queries¶

In [ ]:

# 在单次查询中获取到多个地点
gdf = ox.geocode_to_gdf(["United Kingdom", "Ireland"])

# or optionally buffer them
places = [
    "Berkeley, California, USA",
    "Oakland, California, USA",
    "Piedmont, California, USA",
    "Emeryville, California, USA",
    "Alameda, Alameda County, CA, USA",
]
gdf = ox.geocode_to_gdf(places, buffer_dist=500)

交叉点处理、简化并构建拓扑¶

In [5]:

# get a street network and plot it with all edge intersections
point = 37.858495, -122.267468
G = ox.graph_from_point(point, network_type="drive", dist=500)

In [6]:

fig,ax=ox.plot_graph(G,node_color="r")

In [7]:

# get a GeoSeries of consolidated intersections
G_proj = ox.project_graph(G)
ints = ox.consolidate_intersections(G_proj, rebuild_graph=False, tolerance=15, dead_ends=False)
len(ints)

Out[7]:

66

In [8]:

len(G) # 比较节点数

Out[8]:

109

上面的这种cleaned up intersections得到了更准确的交点数量和密度，但不改变或者整合网络中的拓扑，而是需要我们重建图

In [9]:

# consolidate intersections and rebuild graph topology
# this reconnects edge geometries to the new consolidated nodes
G2 = ox.consolidate_intersections(G_proj, rebuild_graph=True, tolerance=15, dead_ends=False)
len(G2)

Out[9]:

66

In [10]:

fig, ax = ox.plot_graph(G2, node_color="r")

通过对比发现，有的边没了，不过节点确实得到了简化；而且那种四个点的也没了；
设置rebuild_graph=False知识把节点缓冲区内的点合并；设置为True在合并之前还检查了下拓扑

In [ ]:

# 简化图_节点
ox.simplification._is_endpoint
ox.simplify_graph

In [12]:

# 自环边
loops = [edge[0] for edge in nx.selfloop_edges(G2)]
nc = ["r" if node in loops else "y" for node in G2.nodes()]
fig, ax = ox.plot_graph(G2, node_color=nc)

In [15]:

# turn off strict mode and see what nodes we'd remove
nc = ["r" if ox.simplification._is_endpoint(G, node, strict=False) else "y" for node in G.nodes()]
fig, ax = ox.plot_graph(G, node_color=nc)

In [16]:

# simplify network with strict mode turned off
G3 = ox.simplify_graph(G.copy(), strict=False)
fig, ax = ox.plot_graph(G3, node_color="r")

获取网络相关方法中的clean_periphery参数默认是True，表示执行清除外围；具体操作包括：

1. 对于请求的区域做一个0.5km的缓冲区，获取包含缓冲区的大区域内的路网；
2. 简化拓扑，仅存留表示街道交叉口的节点
3. 计算大的网络下的每个节点的无向的度；
4. 裁剪网络到请求的区域大小
5. 保存节点度的数值作为图的属性

In [4]:

# get some bbox
bbox = ox.utils_geo.bbox_from_point((45.518698, -122.679964), dist=300)
north, south, east, west = bbox
G = ox.graph_from_bbox(north, south, east, west, network_type="drive", clean_periphery=False)
fig, ax = ox.plot_graph(G, node_color="r")

In [5]:

# 节点度的分布有很多假的死胡同
k = dict(G.degree())
{n: list(k.values()).count(n) for n in range(max(k.values()) + 1)}

Out[5]:

{0: 0, 1: 30, 2: 2, 3: 6, 4: 52}

In [8]:

list(k.values()).count(4)

Out[8]:

52

In [9]:

G = ox.graph_from_bbox(north, south, east, west, network_type="drive")
fig, ax = ox.plot_graph(G, node_color="r")

相对上面的图，这里清除了stray peripheral edges；在左下角的有两个interstitial nodes被保留了，因为在大的网络中，这两个其实是intersections，而OSMnx也正确的保留了这些节点，在intersection density的求解中也被考虑了进来

In [10]:

# the streets per node distribution for this cleaned up graph is more accurate
# dict keys = count of streets emanating from the node (ie, intersections and dead-ends)
# dict vals = number of nodes with that count
k = nx.get_node_attributes(G, "street_count")
{n: list(k.values()).count(n) for n in range(max(k.values()) + 1)}

Out[10]:

{0: 0, 1: 0, 2: 0, 3: 3, 4: 57}

In [11]:

location_point = (33.299896, -111.831638)
G = ox.graph_from_point(location_point, dist=500, simplify=True)
fig, ax = ox.plot_graph(G, node_color="r")

保存和加载¶

In [13]:

# ox.save_graph_geopackage()
# ox.save_graphml()
# ox.load_graphml()

# 保存到svg
fig,ax=ox.plot_graph(G,show=False,save=True,close=True,filepath="piedont.svg")

In [ ]:

# 保存兴趣点和建筑物足迹数据
# get all "amenities" and save as a geopackage via geopandas
gdf = ox.geometries_from_place(place, tags={"amenity": True})
gdf = gdf.apply(lambda c: c.astype(str) if c.name != "geometry" else c, axis=0)
gdf.to_file("pois.gpkg", driver="GPKG")

In [ ]:

# get all building footprints and save as a geopackage via geopandas
gdf = ox.geometries_from_place(place, tags={"building": True})
gdf = gdf.apply(lambda c: c.astype(str) if c.name != "geometry" else c, axis=0)
gdf.to_file(".building_footprints.gpkg", driver="GPKG")

In [ ]:

# 保存和加载.osm xml文件
ox.settings.all_oneway = True
ox.settings.log_console = True
G = ox.graph_from_place("Piedmont, California, USA", network_type="drive")
ox.save_graph_xml(G, filepath="piedmont.osm")

street network over place shape¶

In [22]:

from descartes import PolygonPatch
from shapely.geometry import MultiPolygon
from shapely.geometry import Polygon

In [27]:

G=ox.load_graphml("mynetwork.graphml")

In [28]:

# plot the network, but do not show it or close it yet
fig, ax = ox.plot_graph(
    G,
    show=False,
    close=False,
    bgcolor="#333333",
    edge_color="w",
    edge_linewidth=0.3,
    node_size=0,
)

In [29]:

gdf=ox.geocode_to_gdf("Qinhuai, Nanjing, China")
gdf

In [32]:

# to this matplotlib axis, add the place shape as descartes polygon patches
for geometry in gdf["geometry"].tolist():
    if isinstance(geometry, (Polygon, MultiPolygon)):
        if isinstance(geometry, Polygon):
            geometry = MultiPolygon([geometry])
        for polygon in geometry.geoms:
            patch = PolygonPatch(polygon, fc="k", ec="#666666", lw=1, alpha=1, zorder=-1)
            ax.add_patch(patch)

# optionally set up the axes extents
margin = 0.02
west, south, east, north = gdf.unary_union.bounds
margin_ns = (north - south) * margin
margin_ew = (east - west) * margin
ax.set_ylim((south - margin_ns, north + margin_ns))
ax.set_xlim((west - margin_ew, east + margin_ew))
plt.show()

这个市政边界是一个行政边界，而不是一个物理边界，所以它代表的是管辖边界，而不是个别的物理特征，如岛屿。要获得单个岛屿的几何图形，使用“geometries”模块:

In [34]:

islands = ox.geometries_from_place("Qinhuai, Nanjing, China", tags={"place": ["island", "islet"]})
islands.shape

Out[34]:

(0, 1)

自定义过滤条件¶

In [ ]:

place = {"city": "Berkeley", "state": "California"}

# only get motorway ways
cf = '["highway"~"motorway"]'
G = ox.graph_from_place(place, network_type="drive", custom_filter=cf)
print(len(G), "motorway")

# only get primary ways
cf = '["highway"~"primary"]'
G = ox.graph_from_place(place, network_type="drive", custom_filter=cf)
print(len(G), "primary")

# use the pipe (|) as 'or' operator
cf = '["highway"~"motorway|primary"]'
G = ox.graph_from_place(place, network_type="drive", custom_filter=cf)
print(len(G), "motorway + primary")

G=ox.graph_from_place(place,custom_filter='["waterway"~"canal"]')

如果下载的数据量太大，OSMnx会把请求分割成多个服务器请求来下载所有数据，然后分配到网络图上

In [ ]:

%%time
# get only motorways, trunks, and their links in all of Belgium
# takes a couple minutes to do all the downloading and processing
# OSMnx automatically divides up the query into multiple requests to not overload server
cf = '["highway"~"motorway|motorway_link|trunk|trunk_link"]'
G = ox.graph_from_place("Belgium", network_type="drive", custom_filter=cf)
fig, ax = ox.plot_graph(G, node_size=0)

In [ ]:

# get NY subway rail network
# note this is rail *infrastructure* and thus includes crossovers, sidings, spurs, yards, etc
# for station-based rail network, you should download a station adjacency matrix elsewhere
ox.settings.useful_tags_way += ["railway"]
G = ox.graph_from_place(
    "New York, New York, USA",
    retain_all=False,
    truncate_by_edge=True,
    simplify=True,
    custom_filter='["railway"~"subway"]',
)

fig, ax = ox.plot_graph(G, node_size=0, edge_color="w", edge_linewidth=0.2)

绘图¶

In [35]:

from IPython.display import Image

In [36]:

# configure the inline image display
img_folder = "images"
extension = "png"
size = 240
dpi = 40

In [38]:

place = "sf"
point = (37.793897, -122.402189)
fp = f"./{img_folder}/{place}.{extension}"
fig, ax = ox.plot_figure_ground(
    point=point, filepath=fp, dpi=dpi, save=True, show=False, close=True
)
Image(fp, height=size, width=size)

Out[38]:

In [ ]:

place = "Atlanta, Georgia, USA"
fp = f"./{img_folder}/{place}.{extension}"
fig, ax = ox.plot_figure_ground(
    address=place,
    network_type="drive",
    filepath=fp,
    dpi=dpi,
    save=True,
    show=False,
    close=True,
)
Image(fp, height=size, width=size)

街道宽度：默认是根据highway属性来进行可视化，也可以进行自定义，没有被自定义的就还是使用默认宽度

In [ ]:

street_widths = {
    "footway": 0.5,
    "steps": 0.5,
    "pedestrian": 0.5,
    "path": 0.5,
    "track": 0.5,
    "service": 2,
    "residential": 3,
    "primary": 5,
    "motorway": 6,
}
place = "sf-custom"
point = (37.793897, -122.402189)
fp = f"./{img_folder}/{place}.{extension}"
fig, ax = ox.plot_figure_ground(
    point=point,
    filepath=fp,
    network_type="all",
    street_widths=street_widths,
    dpi=dpi,
    save=True,
    show=False,
    close=True,
)
Image(fp, height=size, width=size)

下载建筑足迹数据，并作为图底进行可视化

In [ ]:

# configure the inline image display
img_folder = "images"
extension = "png"
size = 240

# specify that we're retrieving building footprint geometries
tags = {"building": True}

gdf = ox.geometries_from_place("Piedmont, California, USA", tags)
gdf_proj = ox.project_gdf(gdf)
fp = f"./{img_folder}/piedmont_bldgs.{extension}"
fig, ax = ox.plot_footprints(gdf_proj, filepath=fp, dpi=400, save=True, show=False, close=True)
Image(fp, height=size, width=size)

In [ ]:

# 保存为shapefile
gdf_save = gdf.applymap(lambda x: str(x) if isinstance(x, list) else x)
gdf_save.drop(labels="nodes", axis=1).to_file("./data/piedmont_bldgs.gpkg", driver="GPKG")

In [ ]:

# 分析建筑足迹数据
# calculate the area in projected units (meters) of each building footprint, then display first five
areas = gdf_proj.area
areas.head()

In [ ]:

# total area (sq m) covered by building footprints
sum(areas)

# get the total area within Piedmont's admin boundary in sq meters
place = ox.geocode_to_gdf("Piedmont, California, USA")
place_proj = ox.project_gdf(place)
place_proj.area.iloc[0]

# what proportion of piedmont is covered by building footprints?
sum(areas) / place_proj.area.iloc[0]

In [ ]:

# 构造一个函数来简化制图
# helper funcion to get one-square-mile street networks, building footprints, and plot them
def make_plot(
    place,
    point,
    network_type="drive",
    dpi=40,
    dist=805,
    default_width=4,
    street_widths=None,
):
    fp = f"./{img_folder}/{place}.{extension}"
    gdf = ox.geometries_from_point(point, tags, dist=dist)
    fig, ax = ox.plot_figure_ground(
        point=point,
        dist=dist,
        network_type=network_type,
        default_width=default_width,
        street_widths=street_widths,
        save=False,
        show=False,
        close=True,
    )
    fig, ax = ox.plot_footprints(
        gdf, ax=ax, filepath=fp, dpi=dpi, save=True, show=False, close=True
    )

place = "portland_buildings"
point = (45.517309, -122.682138)
make_plot(place, point, network_type="all", default_width=1, street_widths={"secondary": 3})

In [39]:

# 用folium创造一个leaflet网络地图
from IPython.display import IFrame

In [41]:

# plot the street network with folium
m1 = ox.plot_graph_folium(G, popup_attribute="name", weight=2, color="#8b0000")n 

In [43]:

# save as html file then display map as an iframe
m1.save(filepath)
IFrame(filepath, width=600, height=500)

In [44]:

# use networkx to calculate the shortest path between two nodes
origin_node = list(G.nodes())[0]
destination_node = list(G.nodes())[-1]
route = nx.shortest_path(G, origin_node, destination_node)

# plot the route with folium
# like above, you can pass keyword args along to folium PolyLine to style the lines
m2 = ox.plot_route_folium(G, route, weight=10)

In [45]:

# save as html file then display map as an iframe
filepath = "route.html"
m2.save(filepath)
IFrame(filepath, width=600, height=500)

Out[45]:

In [46]:

# 把两个组合
# plot the route with folium on top of the previously created graph_map
m3 = ox.plot_route_folium(G, route, route_map=m1, popup_attribute="length", weight=7)

In [48]:

# save as html file then display map as an iframe
m3.save(filepath)
IFrame(filepath, width=600, height=500)

节点高程和边等级¶

osmnx可以使用elevation模块来自行给节点添加高程属性，或者加载本地栅格数据或谷歌高程api作为高程数据源，使用谷歌api的话需要密钥。

有了高程数据后，osmnx可以自动计算边的等级

In [4]:

import sys

In [5]:

# osmnx可以加载一个或一组栅格数据来使用
address = "600 Montgomery St, San Francisco, California, USA"
G = ox.graph_from_address(address=address, dist=500, dist_type="bbox", network_type="bike")
# add node elevations from a single raster file
# some nodes will be null because the single file does not cover the graph's extents
raster_path = "./input_data/elevation1.tif"
G = ox.elevation.add_node_elevations_raster(G, raster_path, cpus=1)

In [ ]:

# add node elevations from multiple raster files
# no nulls should remain
raster_paths = ["./input_data/elevation1.tif", "./input_data/elevation2.tif"]
G = ox.elevation.add_node_elevations_raster(G, raster_paths)
assert not np.isnan(np.array(G.nodes(data="elevation"))[:, 1]).any()

In [6]:

# add edge grades and their absolute values
G = ox.elevation.add_edge_grades(G, add_absolute=True)

使用谷歌api的案例省略

一些指标的统计：
使用网络的无向图方式，这样就不会重复计算双向的街道。对于街道等级，使用绝对值，这样就关注于了街道的阶数而不是它的方向

In [1]:

import osmnx as ox 
import pandas as pd 
import geopandas as gpd 
import networkx as nx 

In [7]:

# calculate the edges' absolute grades (and drop any infinite/null values)
grades = pd.Series([d["grade_abs"] for _, _, d in ox.get_undirected(G).edges(data=True)])
grades = grades.replace([np.inf, -np.inf], np.nan).dropna()

In [9]:

avg_grade = np.mean(grades)
print("Average street grade in {} is {:.1f}%".format(address, avg_grade * 100))

med_grade = np.median(grades)
print("Median street grade in {} is {:.1f}%".format(address, med_grade * 100))

Average street grade in 600 Montgomery St, San Francisco, California, USA is 13.9%
Median street grade in 600 Montgomery St, San Francisco, California, USA is 7.0%

In [10]:

# 利用高程差异绘制节点
# get one color for each node, by elevation, then plot the network
nc = ox.plot.get_node_colors_by_attr(G, "elevation", cmap="plasma")
fig, ax = ox.plot_graph(G, node_color=nc, node_size=5, edge_color="#333333", bgcolor="k")

In [11]:

# 利用级别差异绘制边
# get a color for each edge, by grade, then plot the network
ec = ox.plot.get_edge_colors_by_attr(G, "grade_abs", cmap="plasma", num_bins=5, equal_size=True)
fig, ax = ox.plot_graph(G, edge_color=ec, edge_linewidth=0.5, node_size=0, bgcolor="k")

In [12]:

# 考虑等级差异绘制最短边
# select an origin and destination node and a bounding box around them
origin = ox.distance.nearest_nodes(G, -122.426, 37.77)
destination = ox.distance.nearest_nodes(G, -122.441, 37.773)
bbox = ox.utils_geo.bbox_from_point((37.772, -122.434), dist=1500)

In [13]:

# define some edge impedance function here
def impedance(length, grade):
    penalty = grade**2
    return length * penalty


# add impedance and elevation rise values to each edge in the projected graph
# use absolute value of grade in impedance function if you want to avoid uphill and downhill
for _, _, _, data in G.edges(keys=True, data=True):
    data["impedance"] = impedance(data["length"], data["grade_abs"])
    data["rise"] = data["length"] * data["grade"]

In [15]:

# 发现最小化旅行距离的最短路径
route_by_length = ox.shortest_path(G, origin, destination, weight="length")
fig, ax = ox.plot_graph_route(G, route_by_length, bbox=bbox, node_size=1)

In [16]:

# 发现最小化阻力impedance的最短路径
route_by_impedance = ox.shortest_path(G, origin, destination, weight="impedance")
fig, ax = ox.plot_graph_route(G, route_by_impedance, bbox=bbox, node_size=0)

In [17]:

# 两路径对比
def print_route_stats(route):
    route_grades = ox.utils_graph.get_route_edge_attributes(G, route, "grade_abs")
    msg = "The average grade is {:.1f}% and the max is {:.1f}%"
    print(msg.format(np.mean(route_grades) * 100, np.max(route_grades) * 100))

    route_rises = ox.utils_graph.get_route_edge_attributes(G, route, "rise")
    ascent = np.sum([rise for rise in route_rises if rise >= 0])
    descent = np.sum([rise for rise in route_rises if rise < 0])
    msg = "Total elevation change is {:.1f} meters: {:.0f} meter ascent and {:.0f} meter descent"
    print(msg.format(np.sum(route_rises), ascent, abs(descent)))

    route_lengths = ox.utils_graph.get_route_edge_attributes(G, route, "length")
    print("Total trip distance: {:,.0f} meters".format(np.sum(route_lengths)))

In [18]:

# stats of route minimizing length
print_route_stats(route_by_length)

The average grade is 8.2% and the max is 14.9%
Total elevation change is 11.0 meters: 11 meter ascent and 0 meter descent
Total trip distance: 132 meters

In [19]:

# stats of route minimizing impedance (function of length and grade)
print_route_stats(route_by_impedance)

The average grade is 8.2% and the max is 14.9%
Total elevation change is 11.0 meters: 11 meter ascent and 0 meter descent
Total trip distance: 132 meters

绘制等值线图isochrone map¶

In [20]:

from descartes import PolygonPatch
from shapely.geometry import LineString
from shapely.geometry import Point
from shapely.geometry import Polygon

In [22]:

# configure the place, network type, trip times, and travel speed
place = {"city": "Berkeley", "state": "California"}
network_type = "walk"
trip_times = [5, 10, 15, 20, 25]  # in minutes
travel_speed = 4.5  # walking speed in km/hour

# download the street network
G = ox.graph_from_place(place, network_type=network_type)

In [23]:

# find the centermost node and then project the graph to UTM
gdf_nodes = ox.graph_to_gdfs(G, edges=False)
x, y = gdf_nodes["geometry"].unary_union.centroid.xy
center_node = ox.distance.nearest_nodes(G, x[0], y[0])
G = ox.project_graph(G)

In [24]:

# add an edge attribute for time in minutes required to traverse each edge
meters_per_minute = travel_speed * 1000 / 60  # km per hour to m per minute
for _, _, _, data in G.edges(data=True, keys=True):
    data["time"] = data["length"] / meters_per_minute

绘制步行5，10，15，20等时间内行走的距离

In [26]:

# get one color for each isochrone
iso_colors = ox.plot.get_colors(n=len(trip_times), cmap="plasma", start=0, return_hex=True)
iso_colors

Out[26]:

['#0d0887', '#7e03a8', '#cc4778', '#f89540', '#f0f921']

In [27]:

# color the nodes according to isochrone then plot the street network
node_colors = {}
for trip_time, color in zip(sorted(trip_times, reverse=True), iso_colors):
    subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance="time")
    for node in subgraph.nodes():
        node_colors[node] = color
nc = [node_colors[node] if node in node_colors else "none" for node in G.nodes()]
ns = [15 if node in node_colors else 0 for node in G.nodes()]
fig, ax = ox.plot_graph(
    G,
    node_color=nc,
    node_size=ns,
    node_alpha=0.8,
    edge_linewidth=0.2,
    edge_color="#999999",
)

绘制凸包表示convex hull

In [28]:

# make the isochrone polygons
isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
    subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance="time")
    node_points = [Point((data["x"], data["y"])) for node, data in subgraph.nodes(data=True)]
    bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
    isochrone_polys.append(bounding_poly)

In [29]:

# plot the network then add isochrones as colored descartes polygon patches
fig, ax = ox.plot_graph(
    G, show=False, close=False, edge_color="#999999", edge_alpha=0.2, node_size=0
)
for polygon, fc in zip(isochrone_polys, iso_colors):
    patch = PolygonPatch(polygon, fc=fc, ec="none", alpha=0.6, zorder=-1)
    ax.add_patch(patch)
plt.show()

绘制等值线作为缓冲区，得到比凸包更准确的等值区域

In [30]:

def make_iso_polys(G, edge_buff=25, node_buff=50, infill=False):
    isochrone_polys = []
    for trip_time in sorted(trip_times, reverse=True):
        subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance="time")

        node_points = [Point((data["x"], data["y"])) for node, data in subgraph.nodes(data=True)]
        nodes_gdf = gpd.GeoDataFrame({"id": list(subgraph.nodes)}, geometry=node_points)
        nodes_gdf = nodes_gdf.set_index("id")

        edge_lines = []
        for n_fr, n_to in subgraph.edges():
            f = nodes_gdf.loc[n_fr].geometry
            t = nodes_gdf.loc[n_to].geometry
            edge_lookup = G.get_edge_data(n_fr, n_to)[0].get("geometry", LineString([f, t]))
            edge_lines.append(edge_lookup)

        n = nodes_gdf.buffer(node_buff).geometry
        e = gpd.GeoSeries(edge_lines).buffer(edge_buff).geometry
        all_gs = list(n) + list(e)
        new_iso = gpd.GeoSeries(all_gs).unary_union

        # try to fill in surrounded areas so shapes will appear solid and
        # blocks without white space inside them
        if infill:
            new_iso = Polygon(new_iso.exterior)
        isochrone_polys.append(new_iso)
    return isochrone_polys


isochrone_polys = make_iso_polys(G, edge_buff=25, node_buff=0, infill=True)
fig, ax = ox.plot_graph(
    G, show=False, close=False, edge_color="#999999", edge_alpha=0.2, node_size=0
)
for polygon, fc in zip(isochrone_polys, iso_colors):
    patch = PolygonPatch(polygon, fc=fc, ec="none", alpha=0.7, zorder=-1)
    ax.add_patch(patch)
plt.show()

OSMnx结合igraph¶

networkX是纯python写的，所以虽然使用简单，但面对大型网络的处理分析略显吃力，这里可以转向别的工具，像igraph, graph-tool or cugraph for analysis.

In [31]:

import operator
import igraph as ig
print(ig.__version__)
weight = "length"

0.10.2

In [34]:

G_nx=G
osmids = list(G_nx.nodes)
G_nx = nx.relabel.convert_node_labels_to_integers(G_nx)
# give each node its original osmid as attribute since we relabeled them
osmid_values = {k: v for k, v in zip(G_nx.nodes, osmids)}
nx.set_node_attributes(G_nx, osmid_values, "osmid")

In [35]:

G_nx.nodes(data=True)

In [36]:

%%time
# convert networkx graph to igraph
G_ig = ig.Graph(directed=True)
G_ig.add_vertices(G_nx.nodes)
G_ig.add_edges(G_nx.edges())
G_ig.vs["osmid"] = osmids
G_ig.es[weight] = list(nx.get_edge_attributes(G_nx, weight).values())

Wall time: 32 ms

In [37]:

assert len(G_nx.nodes()) == G_ig.vcount()
assert len(G_nx.edges()) == G_ig.ecount()

In [62]:

# 最短路径
source = list(G_nx.nodes())[26]
target = list(G_nx.nodes())[37]

In [63]:

%%time
path1 = G_ig.get_shortest_paths(v=source, to=target, weights=weight)[0]

Wall time: 0 ns

In [64]:

%%time
path2 = nx.shortest_path(G_nx, source, target, weight=weight)

Wall time: 0 ns

In [65]:

# 照理，面对大数据，igraph是更快的

In [66]:

%%time
closeness1 = G_ig.closeness(vertices=None, mode="ALL", cutoff=None, weights=weight, normalized=True)
max_closeness1 = np.argmax(closeness1)

Wall time: 7.6 s

In [67]:

%%time
closeness2 = nx.closeness_centrality(G_nx, distance=weight, wf_improved=True)
max_closeness2 = max(closeness2.items(), key=operator.itemgetter(1))[0]

Wall time: 2min 22s

In [69]:

# confirm same node has max closeness in both graphs
print(G_nx.nodes[max_closeness2]["osmid"] == G_ig.vs[max_closeness1]["osmid"])

False

高阶绘图¶

使用plot模块

In [68]:

# get n evenly-spaced colors from some matplotlib colormap
ox.plot.get_colors(n=5, cmap="plasma", return_hex=True)

Out[68]:

['#0d0887', '#7e03a8', '#cc4778', '#f89540', '#f0f921']

In [70]:

# get node colors by linearly mapping an attribute's values to a colormap
nc = ox.plot.get_node_colors_by_attr(G, attr="y", cmap="plasma")
fig, ax = ox.plot_graph(G, node_color=nc, edge_linewidth=0.3)

In [71]:

# when num_bins is not None, bin the nodes/edges then assign one color to each bin
# also set equal_size=True for equal-sized quantiles (requires unique bin edges!)
ec = ox.plot.get_edge_colors_by_attr(G, attr="length", num_bins=5)

# otherwise, when num_bins is None (default), linearly map one color to each node/edge by value
ec = ox.plot.get_edge_colors_by_attr(G, attr="length")

# plot the graph with colored edges
fig, ax = ox.plot_graph(G, node_size=5, edge_color=ec, bgcolor="k")

In [72]:

fig, ax = ox.plot_graph(
    G,
    ax=None,  # optionally draw on pre-existing axis
    figsize=(8, 8),  # figure size to create if ax is None
    bgcolor="#111111",  # background color of the plot
    node_color="w",  # color of the nodes
    node_size=15,  # size of the nodes: if 0, skip plotting them
    node_alpha=None,  # opacity of the nodes
    node_edgecolor="none",  # color of the nodes' markers' borders
    node_zorder=1,  # zorder to plot nodes: edges are always 1
    edge_color="#999999",  # color of the edges
    edge_linewidth=1,  # width of the edges: if 0, skip plotting them
    edge_alpha=None,  # opacity of the edges
    show=True,  # if True, call pyplot.show() to show the figure
    close=False,  # if True, call pyplot.close() to close the figure
    save=False,  # if True, save figure to disk at filepath
    filepath=None,  # if save is True, the path to the file
    dpi=300,  # if save is True, the resolution of saved file
    bbox=None,  # bounding box to constrain plot
)

可以用bbox把绘图限制在图的某些精确的部分。例如，想要查看，是否完成了consolidate nearby intersections to clean-up node clusters

In [75]:

place = "Piedmont, California, USA"
G = ox.graph_from_place(place, network_type="drive")

In [76]:

Gc = ox.consolidate_intersections(ox.project_graph(G), dead_ends=True)
c = ox.graph_to_gdfs(G, edges=False).unary_union.centroid
bbox = ox.utils_geo.bbox_from_point(point=(c.y, c.x), dist=200, project_utm=True)
fig, ax = ox.plot_graph(
    Gc,
    figsize=(5, 5),
    bbox=bbox,
    edge_linewidth=2,
    edge_color="c",
    node_size=80,
    node_color="#222222",
    node_edgecolor="c",
)

# or save a figure to disk instead of showing it
fig, ax = ox.plot_graph(G, filepath="./images/image.png", save=True, show=False, close=True)

In [102]:

# 绘制路线
# impute missing edge speeds and calculate free-flow travel times
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# calculate 3 shortest paths, minimizing travel time
w = "travel_time"
orig, dest = list(G)[10], list(G)[-1]
route1 = ox.shortest_path(G, orig, dest, weight=w)
orig, dest = list(G)[0], list(G)[-1]
route2 = ox.shortest_path(G, orig, dest, weight=w)
orig, dest = list(G)[2], list(G)[9]
route3 = ox.shortest_path(G, orig, dest, weight=w)

In [103]:

route3

Out[103]:

[53018399,
 53078914,
 53085009,
 682821166,
 53119160,
 53076382,
 53064331,
 53119162,
 53081990,
 53119164,
 53041443,
 53042707,
 53119166,
 53119168,
 53085388,
 53085389,
 53073689,
 53025060,
 53082627,
 53073691,
 53033660,
 53040495,
 53085375,
 53085377,
 53027454,
 53085380,
 53021750,
 53021742,
 53021743]

In [104]:

fig, ax = ox.plot_graph_route(G, route1, orig_dest_size=0, node_size=0)

In [105]:

# you can also pass any ox.plot_graph parameters as additional keyword args
fig, ax = ox.plot_graph_route(G, route1, save=True, show=False, close=True)

In [106]:

# 绘制多重路线
routes = [route1, route2, route3]
rc = ["r", "y", "c"]
fig, ax = ox.plot_graph_routes(G, routes, route_colors=rc, route_linewidth=6, node_size=0)

标注街道段名字

In [109]:

G2=ox.get_undirected(Gc)

In [110]:

fig, ax = ox.plot_graph(G2, edge_linewidth=3, node_size=0, show=False, close=False)
for _, edge in ox.graph_to_gdfs(G2, nodes=False).fillna("").iterrows():
    text = edge["name"]
    c = edge["geometry"].centroid
    ax.annotate(text, (c.x, c.y), c="y")
plt.show()

In [ ]:

# 绘制地理实体，plot_gootprints
# get all the building footprints in a city
gdf = ox.geometries_from_place("Piedmont, California, USA", {"building": True})
gdf.shape
fig, ax = ox.plot_footprints(gdf)

# or plot street network and the entities' footprints together
fig, ax = ox.plot_footprints(gdf, alpha=0.4, show=False)
fig, ax = ox.plot_graph(G, ax=ax, node_size=0, edge_color="w", edge_linewidth=0.7)

下载地理实体¶

可以使用geometries模块来下载实体，例如grocery stores, transit stops, POI or building footprints，并转变为GeoDataFrame

In [ ]:

# get all building footprints in some neighborhood
# `True` means retrieve any object with this tag, regardless of value
place = "Bunker Hill, Los Angeles, California"
tags = {"building": True}
gdf = ox.geometries_from_place(place, tags)
gdf.shape

In [ ]:

fig, ax = ox.plot_footprints(gdf, figsize=(3, 3))

In [ ]:

# get all the parks in some neighborhood
# constrain acceptable `leisure` tag values to `park`
tags = {"leisure": "park"}
gdf = ox.geometries_from_place(place, tags)
gdf.shape

In [ ]:

# get everything tagged amenity,
# and everything tagged landuse = retail or commercial,
# and everything tagged highway = bus_stop
tags = {"amenity": True, "landuse": ["retail", "commercial"], "highway": "bus_stop"}
gdf = ox.geometries_from_place("Piedmont, California, USA", tags)
gdf.shape

In [ ]:

# view just the banks
gdf[gdf["amenity"] == "bank"].dropna(axis=1, how="any")

In [ ]:

# view just the bus stops
gdf[gdf["highway"] == "bus_stop"].dropna(axis=1, how="any").head()

街道网络的方向¶

In [111]:

# define the study sites as label : query
places = {
    # 'Atlanta'       : 'Atlanta, Georgia, USA',
    # 'Boston'        : 'Boston, MA, USA',
    "Buffalo": "Buffalo, NY, USA",
    # 'Charlotte'     : 'Charlotte, NC, USA',
    # 'Chicago'       : 'Chicago, IL, USA',
    "Cleveland": "Cleveland, OH, USA",
    # 'Dallas'        : 'Dallas, TX, USA',
    # 'Houston'       : 'Houston, TX, USA',
    # 'Denver'        : 'Denver, CO, USA',
    # 'Detroit'       : 'Detroit, MI, USA',
    # 'Las Vegas'     : 'Las Vegas, NV, USA',
    # 'Los Angeles'   : {'city':'Los Angeles', 'state':'CA', 'country':'USA'},
    # 'Manhattan'     : 'Manhattan, NYC, NY, USA',
    "Miami": "Miami, FL, USA",
    "Minneapolis": "Minneapolis, MN, USA",
    # 'Orlando'       : 'Orlando, FL, USA',
    # 'Philadelphia'  : 'Philadelphia, PA, USA',
    # 'Phoenix'       : 'Phoenix, AZ, USA',
    # 'Portland'      : 'Portland, OR, USA',
    # 'Sacramento'    : 'Sacramento, CA, USA',
    "San Francisco": {"city": "San Francisco", "state": "CA", "country": "USA"},
    # 'Seattle'       : 'Seattle, WA, USA',
    # 'St Louis'      : 'St. Louis, MO, USA',
    # 'Tampa'         : 'Tampa, FL, USA',
    "Washington": "District of Columbia, USA",
}

In [112]:

# verify OSMnx geocodes each query to what you expect (i.e., a [multi]polygon geometry)
gdf = ox.geocode_to_gdf(list(places.values()))
gdf

Out[112]:

	geometry	bbox_north	bbox_south	bbox_east	bbox_west	place_id	osm_type	osm_id	lat	lon	display_name	class	type	importance
0	POLYGON ((-78.92236 42.94914, -78.91945 42.947...	42.966449	42.826006	-78.795121	-78.922357	307947630	relation	175031	42.886717	-78.878392	Buffalo, Erie County, New York, United States	boundary	administrative	0.732931
1	POLYGON ((-81.87909 41.39641, -81.87906 41.395...	41.604436	41.390628	-81.532744	-81.879094	307598444	relation	182130	41.499657	-81.693677	Cleveland, Cuyahoga County, Ohio, United States	boundary	administrative	0.856081
2	POLYGON ((-80.31976 25.76249, -80.31968 25.762...	25.855783	25.709052	-80.139157	-80.319760	307919581	relation	1216769	25.774173	-80.193620	Miami, Miami-Dade County, Florida, United States	boundary	administrative	0.885449
3	POLYGON ((-93.32913 44.92030, -93.32911 44.920...	45.051250	44.890150	-93.193859	-93.329127	307870876	relation	136712	44.977300	-93.265469	Minneapolis, Hennepin County, Minnesota, Unite...	boundary	administrative	0.747169
4	MULTIPOLYGON (((-123.17382 37.77573, -123.1737...	37.929811	37.640314	-122.281479	-123.173825	307775847	relation	111968	37.779026	-122.419906	San Francisco, California, United States	boundary	administrative	1.025131
5	POLYGON ((-77.11979 38.93435, -77.11977 38.934...	38.995968	38.791630	-76.909366	-77.119795	307700086	relation	162069	38.893847	-76.988043	District of Columbia, United States	boundary	administrative	0.712066

In [113]:

# create figure and axes
n = len(places)
ncols = int(np.ceil(np.sqrt(n)))
nrows = int(np.ceil(n / ncols))
figsize = (ncols * 5, nrows * 5)
fig, axes = plt.subplots(nrows, ncols, figsize=figsize, subplot_kw={"projection": "polar"})

# plot each city's polar histogram
for ax, place in zip(axes.flat, sorted(places.keys())):
    print(ox.utils.ts(), place)

    # get undirected graphs with edge bearing attributes
    G = ox.graph_from_place(place, network_type="drive")
    Gu = ox.add_edge_bearings(ox.get_undirected(G))
    fig, ax = ox.bearing.plot_orientation(Gu, ax=ax, title=place, area=True)

# add figure title and save image
suptitle_font = {
    "family": "DejaVu Sans",
    "fontsize": 60,
    "fontweight": "normal",
    "y": 1,
}
fig.suptitle("City Street Network Orientation", **suptitle_font)
fig.tight_layout()
fig.subplots_adjust(hspace=0.35)
fig.savefig("images/street-orientations.png", facecolor="w", dpi=100, bbox_inches="tight")
plt.close()

2023-01-28 09:36:11 Buffalo
2023-01-28 09:37:29 Cleveland
2023-01-28 09:38:59 Miami

也可以通过ox.bearing.orientation_entropy函数来计算空间图的orientation entropy