用Python实现城市方向熵的计算

• 引言

第一次写CSDN的博客，就先给大家介绍如何安装osmnx模块（让人很头疼），以及利用osmnx和geopandas采集openstreetmap上的城市街道数据并且计算和可视化街道的方向熵，教程的主要方法来自微信公众号“数读城事”，本博客就把你可能遇到的问题以及实现的具体步骤复现一遍，废话不多说，咱们开始吧。

• 安装osmnx

笔者使用的是Python 3.7版本，必须要吐槽的是osmnx的模块真的很难安装，网上有很多方法讲解如何安装，主流的IDLE我试过了anaconda、pycharm和python自带的IDLE，其中anaconda最难安装，不仅会遇到网络问题、编译环境问题，而且安装完成后调用失败，在这里就不推荐给大家了，pycharm可以直接使用python的包，所以，简单省事就选pip吧。这里直接附上链接，window中osmnx包的详细安装过程，大家按图索骥去安装吧，国内使用记得在安装时换国内的镜像源哦。

• 行政区划数据准备

这里我们以成都主城区为例，需要说明的是，分区县的矢量数据是1：100万的全国基础地理信息数据，由数读城事整理，数据已上传，需要的小伙伴自行下载，我们使用ArcGIS对数据进行选择和导出。

首先打开GIS。

载入区县一级的行政区划数据。

选择——按属性选择——“市”=“成都市”。

导出——导出所选要素，格式为.shp。

完毕。

• OSM数据的获取

我把接下去的步骤分为数据的获取以及可视化两个部分。数据的获取部分直接附上整理后的代码。

import osmnx as ox
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import math
from plot_map import plot_map

def reverse_bearing(x):
    return x + 180 if x < 180 else x - 180

def 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 bearings
def 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)

if __name__ == '__main__' :
    city = gpd.read_file('./成都.shp')
    city_road = gpd.GeoDataFrame(columns = ['u', 'v', 'key', 'osmid', 'name', 'highway', 'oneway', 'length',
       'geometry', 'bridge', 'ref', 'lanes', 'maxspeed', 'access', 'tunnel',
       'junction','district'])  ### 创造一个空的GeoDataFrame，用来储存我们即将下载的各个区的路网数据
    city_orientation_count = []  ### 创造一个空的列表，用来存放各个区内各个街道的方向
    try:
        for i in range(0,10):  ### 写一个循环，遍历每个区，依次获取并处理路网数据
            G = ox.graph_from_polygon(city.loc[i,'geometry'],network_type='drive')  ### 利用osmnx，获取每个区的路网数据（可以通行小汽车的路网）
            city.loc[i,'三岔路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][3]
            city.loc[i,'十字路口比例'] = ox.stats.basic_stats(G)['streets_per_node_proportion'][4]
            road_gdf = ox.graph_to_gdfs(G)[1]  ### 上一步获取的路网数据格式为networkx中的graph格式，这里我们将它转换成GeoDataFrame
            road_gdf['district'] = city.loc[i,'NAME']  ### 并为路网赋上相应的行政区信息（属于哪个区）
            city_road = city_road.append(road_gdf,ignore_index=True)  ### 将每个区的路网添加至总的路网数据中
            Gu = ox.add_edge_bearings(ox.get_undirected(G))  ### 利用osmnx的add_edge_bearings函数为路网的边添加方向属性
            b = pd.Series([d['bearing'] for u, v, k, d in Gu.edges(keys=True, data=True)])  ### 将边的方向属性都提取出来，存在一个Series中
            b = pd.concat([b, b.map(reverse_bearing)]).reset_index(drop=True).dropna()  ### 为边添加另一个方向的方向属性（+-180度）（因为路都是直线，如果从a端点到b端点与正北的夹角为30度，那么b端点到a端点与正北的夹角就是210度
            city_orientation_count.append(b)  ### 将提取出来的方向属性添加到总的方向数据中
            city.loc[i,'方向熵'] = calculate_entropy(count_and_merge(36,b))  ### 计算每个区的街道的方向熵，并直接储存在GeoDataFrame中
            print('{}处理完毕'.format(city.loc[i,'NAME']))
    except Exception as e:
        print(e)

 

代码中的try和except可以去掉，笔者的电脑不知道有什么问题，总是在循环跑到50%的地方程序死掉，然后直接重启，现在依旧如此，解决办法就是把数据多添加周围的一些区县，增大循环的次数，把需要的部分留到50%之前。

然后只需要run module就可以等待着下载数据了，下载时间视网络情况来定。

接下来就准备出图吧。

• 可视化

上代码。

city_district_entropy = city.loc[:,['NAME','方向熵']]  ### 提取街道名称与方向熵字段
    city_district_entropy.columns = ['district','方向熵']  ### 为了与路网数据的GeoDataFrame统一，改一下表的名字，将街道名称字段的列名从NAME改为district，方便下面进行连接
    city_road = pd.merge(city_road,city_district_entropy,on='district',how='outer')  ### 根据district字段进行连接
    city_bounds = [102.949, 30.052, 104.847, 31.489 ]
    base = city.plot(figsize=(10,10),facecolor='none',edgecolor='black',lw=0.6)  #### 开始绘图，先设置行政区作为底图
    city_road.plot(ax=base,column='方向熵',lw=0.4,legend=True,legend_kwds={'shrink':0.7,'label':'Entropy'})  ### 然后在底图上绘制路网地图，并以每条道路所在区的方向熵大小为路网赋予颜色（cmap操作）
    for i in city.index:  ### 标记方向熵大于75%分位数的区的名称
        if city.loc[i,'方向熵'] >= city['方向熵'].quantile(0.75):
            plt.text(city.loc[i,'geometry'].centroid.x,city.loc[i,'geometry'].centroid.y,city.loc[i,'NAME'],fontdict={'family':'SimHei','size':10},horizontalalignment='center',verticalalignment='center')
    plot_map(plt,city_bounds,zoom=14,style=4)  ### 调用plot_map添加osm的底图
    base.axis('off')  ### 关闭坐标轴，更好看
    plt.savefig('city_entropy_plot.jpg',dpi=100,bbox_inches='tight')
    plt.show() 

结果。

再来一个。

fig, axes = plt.subplots(2,5, figsize=(20,20),subplot_kw={'projection':'polar'})
    fig.tight_layout(h_pad=5)
    for i in range(2):
        for x in range(5):
            polar_plot(axes[i][x],city_orientation_count[i*2+x])
            axes[i][x].set_title(city.loc[i*5+x,'NAME'],fontdict={'family':'SimHei','weight':'heavy','size':18})
    plt.savefig('city_street_orientation.jpg',dpi=100,bbox_inches='tight')
    plt.show()    

结果。

以上。

• 参考文章

​​​​​​​【城事数说】方向熵：上海与重庆的道路方向与混乱程度​​​​​​​