其他
【城事数说】方向熵：上海与重庆的道路方向与混乱程度
作者：YuanLiang
今天的原创投稿是由公众号“Yuan的数据分析”的作者YuanLiang带来的纯干货技术流文章。文末附上他的公众号二维码，希望大家可以关注。
城市作为人类空间活动的载体，有着丰富的空间形态与结构机理，众多研究者对这一命题展开了丰富的研究，
然而现有研究多从建筑学及设计学的角度出发，很少有研究从数据科学的角度对城市街道的形态与结构进行分析。
（好了我编不下去了
本篇教程的主要目的在于介绍给大家利用osmnx以及geopandas采集openstreetmap上的城市街道数据，计算并可视化街道的方向熵(orientation entropy)
具体的方法大家可以关注osmnx的github主页：https://github.com/gboeing/osmnx
也可以读一下Geoff Boeing 2019年发表在Applied Network Science上的文章
Urbanspatial order: street network orientation, configuration, and entropy
DOI:https://doi.org/10.1007/s41109-019-0189-1
本教程实际上是对这篇文章的一个小点的复现
PS：这里的方向熵是指一个区域内所有街道的方向的混乱程度。举个例子，大家都知道北京、西安这种皇城，道路一般都十分规整（体现了权威与秩序），属于棋盘格式的布局，路痴也不容易迷路，反正就东南西北四个方向，对应的方向熵就会很低（不混乱）；反观重庆之类的城市，城市形态十分diversity，各种路纵横交错，路痴肯定迷路，对应的方向熵就会很高（很混乱）
⾔归正传我们开始分析
⾸先调包，此处为⼩旭学⻓的开发的plot_map打call，强烈推荐！超好⽤！！！
import osmnx as oximport matplotlib.pyplot as pltimport geopandas as gpdimport pandas as pdimport numpy as npfrom plot_map import plot_mapimport math下面是四个计算街道方向熵以及绘图的函数，大家不用理解，直接拿来用即可
其中方向熵的计算公式为：
其实就是香农熵的计算公式，P（Oi）指方向为第i个方向的街道的比例（这里将街道的方向分为了36类，将正北设置为0度，以10度为间隔创建36个区间，街道的方向落在哪个区间，我们就认为它的方向就在这个区间）
def reverse_bearing(x): return x + 180 if x < 180 else x - 180def count_and_merge(n, bearings): # make twice as many bins as desired, then merge them in pairs # prevents bin-edge effects around common values like 0° and 90° n = n * 2 bins = np.arange(n + 1) * 360 / n count, _ = np.histogram(np.array(bearings), bins=bins) # move the last bin to the front, so eg 0.01° and 359.99° will be binned together count = np.roll(count, 1) return count[::2] + count[1::2]def calculate_entropy(count): count_p = count/count.sum() entropy = 0 for i in count_p: entropy -= i*math.log(i) return entropy # function to draw a polar histogram for a set of edge bearingsdef polar_plot(ax, bearings, n=36, title=''):
bins = np.arange(n + 1) * 360 / n count = count_and_merge(n, bearings) _, division = np.histogram(bearings, bins=bins) frequency = count / count.sum() division = division[0:-1] width = 2 * np.pi / n
ax.set_theta_zero_location('N') ax.set_theta_direction('clockwise')
x = division * np.pi / 180 bars = ax.bar(x, height=frequency, width=width, align='center', bottom=0, zorder=2, color='#003366', edgecolor='k', linewidth=0.5, alpha=0.7) ax.set_ylim(top=frequency.max()) title_font = {'family':'DejaVu Sans', 'size':24, 'weight':'bold'} xtick_font = {'family':'DejaVu Sans', 'size':10, 'weight':'bold', 'alpha':1.0, 'zorder':3} ytick_font = {'family':'DejaVu Sans', 'size': 9, 'weight':'bold', 'alpha':0.2, 'zorder':3} ax.set_title(title.upper(), y=1.05, fontdict=title_font) ax.set_yticks(np.linspace(0, max(ax.get_ylim()), 5)) yticklabels = ['{:.2f}'.format(y) for y in ax.get_yticks()] yticklabels[0] = '' ax.set_yticklabels(labels=yticklabels, fontdict=ytick_font) xticklabels = ['N', '', 'E', '', 'S', '', 'W', ''] ax.set_xticklabels(labels=xticklabels, fontdict=xtick_font) ax.tick_params(axis='x', which='major', pad=-2)首先导入上海市的shp地图
sh = gpd.read_file('./上海行政区/区县界_area.shp')下面开始下载上海市各个区的路网数据，并计算每个区的路网的方向熵
### 创造一个空的GeoDataFrame，用来储存我们即将下载的各个区的路网数据sh_road = gpd.GeoDataFrame(columns = ['u', 'v', 'key', 'osmid', 'name', 'highway', 'oneway', 'length', 'geometry', 'bridge', 'ref', 'lanes', 'maxspeed', 'access', 'tunnel', 'junction','district'])
### 创造一个空的列表，用来存放各个区内各个街道的方向sh_orientation_count = []
### 写一个循环，遍历上海的每个区，依次获取并处理路网数据for i in sh.index:
### 利用osmnx，获取每个区的路网数据（可以通行小汽车的路网） G = ox.graph_from_polygon(sh.loc[i,'geometry'],network_type='drive') #sh.loc[i,'三岔路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][3] #sh.loc[i,'十字路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][4]
### 上一步获取的路网数据格式为networkx中的graph格式，这里我们将它转换成GeoDataFrame road_gdf = ox.graph_to_gdfs(G)[1]
### 并为路网赋上相应的行政区信息（属于哪个区） road_gdf['district'] = sh.loc[i,'NAME'] ### 将每个区的路网添加至总的路网数据中 sh_road = sh_road.append(road_gdf,ignore_index=True)
### 利用osmnx的add_edge_bearings函数为路网的边添加方向属性 Gu = ox.add_edge_bearings(ox.get_undirected(G))
### 将边的方向属性都提取出来，存在一个Series中 b = pd.Series([d['bearing'] for u, v, k, d in Gu.edges(keys=True, data=True)])
### 为边添加另一个方向的方向属性（+-180度）（因为路都是直线，如果从a端点到b端点与正北的夹角为30度，那么b端点到a端点与正北的夹角就是210度 b = pd.concat([b, b.map(reverse_bearing)]).reset_index(drop=True).dropna()
### 将提取出来的方向属性添加到总的方向数据中 sh_orientation_count.append(b)
### 计算每个区的街道的方向熵，并直接储存在上海的GeoDataFrame中 sh.loc[i,'方向熵'] = calculate_entropy(count_and_merge(36,b))
print('{}处理完毕'.format(sh.loc[i,'NAME']))经过上一步的操作，我们可以得到三个数据：
上海的GeoDataFrame，储存了每个区的geometry和方向熵的具体数值
上海的路网数据GeoDataFrame，储存了上海每个区的路网数据，有一个district字段专门来标记每条道路属于哪个区
上海各个区的道路方向的集合（列表形式），储存了每个区内所有路网的方向
下面我们开始可视化
首先我们从上海的GeoDataFrame中提取每个区的方向熵，然后根据区的名称的字段，将其与上海路网数据的GeoDataFrame进行连接
这样我们可以为路网数据的GeoDataFrame添加一个方向熵的字段，代表每条道路所在区的路网的方向熵
### 提取街道名称与方向熵字段sh_district_entropy = sh.loc[:,['NAME','方向熵']]### 为了与路网数据的GeoDataFrame统一，改一下表的名字，将街道名称字段的列名从NAME改为district，方便下面进行连接sh_district_entropy.columns = ['district','方向熵']### 根据district字段进行连接sh_road = pd.merge(sh_road,sh_district_entropy,on='district',how='outer')sh_bounds = [120.79560279980944,30.623112899720564,122.03020320013769, 31.925845100091397]#### 开始绘图，先设置上海的行政区作为底图base = sh.plot(figsize=(10,10),facecolor='none',edgecolor='black',lw=0.6)
### 然后在底图上绘制路网地图，并以每条道路所在区的方向熵大小为路网赋予颜色（cmap操作）sh_road.plot(ax=base,column='方向熵',lw=0.4,legend=True,legend_kwds={'shrink':0.7,'label':'Entropy'})
### 标记方向熵大于75%分位数的区的名称for i in sh.index: if sh.loc[i,'方向熵'] >= sh['方向熵'].quantile(0.75): plt.text(sh.loc[i,'geometry'].centroid.x,sh.loc[i,'geometry'].centroid.y,sh.loc[i,'NAME'],fontdict={'family':'Arial Unicode MS','size':10},horizontalalignment='center',verticalalignment='center')
### 调用plot_map添加osm的底图plot_map(plt,sh_bounds,zoom=14,style=4)
### 关闭坐标轴，更好看base.axis('off')
### 保存图片plt.savefig('sh_entropy_plot.jpg',dpi=100,bbox_inches='tight')
plt.show()imgsavepath do not exist, your tile map will not save这样上海市的路网情况就可视化出来了，可以看到，虹口区的路网方向熵最大，说明路网最乱，路痴最容易迷路；普陀区其次，杨浦、长宁与金山也处在前几位。
以上只是最简单的可视化，我们还可以做的骚
fig, axes = plt.subplots(4,4, figsize=(20,20),subplot_kw={'projection':'polar'})fig.tight_layout(h_pad=5)for i in range(4): for x in range(4): polar_plot(axes[i][x],sh_orientation_count[i*4+x]) axes[i][x].set_title(sh.loc[i*4+x,'NAME'],fontdict={'family':'Arial Unicode MS','weight':'heavy','size':18})plt.savefig('sh_street_orientation.jpg',dpi=100,bbox_inches='tight')plt.show()可以看到，果不其然，虹口的道路各个方向的都有，难怪方向熵那么大，反观闵行，道路方向的分布十分集中，说明闵行的路网多为方格式布局，十分规整。
上海的可视化就到此为止，下面我们来可视化一下5D城市重庆的街道方向熵。
具体的代码跟上海一样，我就不加注释了。
cq = gpd.read_file('./重庆市/区县界_area.shp')cq = cq.iloc[[13,14,15,26,31,32,33,34,35],:]cq.reset_index(drop=True,inplace=True)cq_road = gpd.GeoDataFrame(columns = ['u', 'v', 'key', 'osmid', 'name', 'highway', 'oneway', 'length', 'geometry', 'bridge', 'ref', 'lanes', 'maxspeed', 'access', 'tunnel', 'junction','district'])cq_orientation_count = []for i in cq.index:
G = ox.graph_from_polygon(cq.loc[i,'geometry'],network_type='drive') #cq.loc[i,'三岔路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][3] #cq.loc[i,'十字路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][4]
road_gdf = ox.graph_to_gdfs(G)[1] road_gdf['district'] = cq.loc[i,'NAME'] cq_road = cq_road.append(road_gdf,ignore_index=True)
Gu = ox.add_edge_bearings(ox.get_undirected(G)) b = pd.Series([d['bearing'] for u, v, k, d in Gu.edges(keys=True, data=True)]) b = pd.concat([b, b.map(reverse_bearing)]).reset_index(drop=True).dropna() cq_orientation_count.append(b)
cq.loc[i,'方向熵'] = calculate_entropy(count_and_merge(36,b))
print('{}处理完毕'.format(cq.loc[i,'NAME']))
cq_district_entropy = cq.loc[:,['NAME','方向熵']]cq_district_entropy.columns = ['district','方向熵']cq_road = pd.merge(cq_road,cq_district_entropy,on='district',how='outer')cq_bounds = [106.20977790026993,29.079847999725487,107.03557209963279, 30.17672399987273]base = cq.plot(figsize=(10,10),facecolor='none',edgecolor='black',lw=0.6)cq_road.plot(ax=base,column='方向熵',lw=0.4,legend=True,legend_kwds={'shrink':0.7,'label':'Entropy'})for i in cq.index: if cq.loc[i,'方向熵'] >= cq['方向熵'].quantile(0.75): plt.text(cq.loc[i,'geometry'].centroid.x,cq.loc[i,'geometry'].centroid.y,cq.loc[i,'NAME'],fontdict={'family':'Arial Unicode MS','size':10},horizontalalignment='center',verticalalignment='center')plot_map(plt,cq_bounds,zoom=14,style=4)base.axis('off')plt.savefig('cq_entropy_plot.jpg',dpi=100,bbox_inches='tight')plt.show()imgsavepath do not exist, your tile map will not savefig, axes = plt.subplots(3,3, figsize=(15,15),subplot_kw={'projection':'polar'})fig.tight_layout(h_pad=5)for i in range(3): for x in range(3): polar_plot(axes[i][x],cq_orientation_count[i*3+x]) axes[i][x].set_title(cq.loc[i*3+x,'NAME'],fontdict={'family':'Arial Unicode MS','weight':'heavy','size':18})plt.savefig('cq_street_orientation.jpg',dpi=100,bbox_inches='tight')哈哈哈哈哈哈哈，可以看到，渝中区太魔鬼了，路痴的噩梦，这谁能顶得住。
本篇教程到这里就结束啦，关于更多的信息，请欢迎大家关注我的公众号【Yuan的数据分析】以及github主页：https://github.com/YuanLianggg
还有小旭学长的python教程：https://github.com/ni1o1/pygeo-tutorial （强烈推荐！！！）
【Yuan的数据分析】扫描二维码关注
文章到此结束，下面是我个人的小广告==
话说我刚刚整个了微博==
虽然还没发啥东西==
大家要不要来支持下
使用关键词搜索历史文章请点击【阅读原文】