本文共 3940 字，大约阅读时间需要 13 分钟。
本文完整代码及数据已上传至我的GitHub仓库。
OD数据是交通、城市规划以及GIS等领域常见的一类数据，特点是每一条数据都记录了一次OD（Origin即起点，Destination即终点）行为的起点与终点坐标信息。而针对OD数据的常见可视化表达方式为弧线图，例如纽约曼哈顿等区域的某时间段Uber打车记录上下车点数据的展示方式（如图1所示）。
但这种传统的表达方式存在局限性：当OD记录数量众多时，由于不同线之间的彼此堆叠，导致很多区域之间的OD模式被遮盖而难以被读出。前一段时间，我在观看一场学术直播的过程中，注意到一种特别的表达区域间OD数据的方式，原始文献比较老（发表于2010年），其思想是通过对研究区域进行网格化划分，再将整个区域的原始网格映射到每个单一网格中（如图2所示）。
以坐标为( E,5 )的网格为例，到达坐标为( A,2 )的网格的所有OD数据记录，可以在右图中对应左图( E,5 )位置的大网格中，划分出的对应( A,2 )相对位置的小网格中进行记录。通过这种方式，原始文献将图3所示原始OD线图转换为图4所示的结果，使得我们可以非常清楚地观察到每个网格区域对其他网格区域的OD模式。本文将利用Python，在图1对应的Uber上下车点分布数据的基础上，实践这种表达OD数据的特别方式。
首先，我们需要梳理整体逻辑。原始数据中我们主要用到的字段包括上车点经纬度pickup_longitude与pickup_latitude，以及下车点经纬度dropoff_longitude与dropoff_latitude。我的思路是首先对所有经纬度点进行去重，然后保存为GeoDataFrame并统一坐标参考系为Web墨卡托（EPSG:3857）。
具体实现如下：
from shapely.geometry import Pointfrom geopandas import GeoDataFrameod_points = pd.concat([ taxi_trip_flow[['pickup_longitude', 'pickup_latitude']] .rename(columns={'pickup_longitude': 'lng', 'pickup_latitude': 'lat'}), taxi_trip_flow[['dropoff_longitude', 'dropoff_latitude']] .rename(columns={'dropoff_longitude': 'lng', 'dropoff_latitude': 'lat'})], ignore_index=True)od_points = od_points.drop_duplicates()od_points['geometry'] = od_points.apply(lambda row: Point(row['lng'], row['lat']), axis=1)od_points = gpd.GeoDataFrame(od_points, crs='EPSG:4326').to_crs('EPSG:3857')od_points.head() 接下来，我们为研究区域创建网格面矢量数据。思路是利用numpy先创建出x和y方向上的等间距坐标，例如我们这里创建5行5列：
from shapely.geometry import MultiLineStringfrom shapely.ops import polygonizexmin, ymin, xmax, ymax = od_points.total_boundsx = np.linspace(xmin, xmax, 6)y = np.linspace(ymin, ymax, 6)hlines = [((x1, yi), (x2, yi)) for x1, x2 in zip(x[:-1], x[1:]) for yi in y]vlines = [((xi, y1), (xi, y2)) for y1, y2 in zip(y[:-1], y[1:]) for xi in x]manhattan_grids = GeoDataFrame({ 'geometry': list(polygonize(MultiLineString(hlines + vlines)))}, crs='EPSG:3857')manhattan_grids['id'] = manhattan_grids.index 创建出的网格效果如图7所示。接下来是最关键的部分：计算每个原始网格内部的上车点记录，并利用仿射变换得到镶嵌在其内部的小网格。
以id=21的网格为例（对应肯尼迪国际机场的区域），我们需要完成以下步骤：
i = 21center_grid = (manhattan_grids.unary_union.centroid.x, manhattan_grids.unary_union.centroid.y)dropoff = ( gpd.sjoin(manhattan_grids.loc[i:i, :], od_points, op='contains') [['lng', 'lat', 'geometry']] .merge(taxi_trip_flow[['pickup_longitude', 'pickup_latitude', 'dropoff_longitude', 'dropoff_latitude']], left_on=['lng', 'lat'], right_on=['pickup_longitude', 'pickup_latitude']) [['dropoff_longitude', 'dropoff_latitude']] .merge(od_points, left_on=['dropoff_longitude', 'dropoff_latitude'], right_on=['lng', 'lat']) [['geometry']])grid_distrib = ( gpd.sjoin(manhattan_grids, GeoDataFrame(dropoff, crs='EPSG:3857'), op='contains') .groupby('id', as_index=False) .agg({'index_right': 'count'}) .rename(columns={'index_right': '下车记录数'}))grid_distrib.head() 通过上述步骤，我们得到了id为21的网格下车点分布结果。将上述过程推广到每个网格，并将最后的计算结果合并为一张GeoDataFrame，即表draw_base。
最终，我们对draw_base表进行可视化。为了显示更加自然，对下车记录进行了对数化+自然间断处理：
import matplotlib.pyplot as pltimport contextily as ctxfig, ax = plt.subplots(figsize=(12, 12), dpi=200)draw_base.plot(facecolor='none', edgecolor='lightgrey', ax=ax, linewidth=0.3)draw_base.assign(下车记录数=np.log(draw_base.下车记录数))draw_base.plot(column='下车记录数', scheme='NaturalBreaks', k=5, cmap='YlOrRd', ax=ax, alpha=0.7)manhattan_grids.plot(ax=ax, facecolor='none', edgecolor='black', linewidth=0.8)draw_base.query('是否为目标网格 == 1').plot(facecolor='none', edgecolor='black', linestyle='--', ax=ax)minx, miny, maxx, maxy = manhattan_grids.total_boundsax.set_xlim(minx, maxx)ax.set_ylim(miny, maxy)ax.axis('off')ctx.add_basemap(ax, source='https://d.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png', zoom=12)fig.savefig('图10.png', dpi=500, bbox_inches='tight', pad_inches=0) 通过这种表达方式，我们可以很明显地看出不同区域相对其他区域出行模式的不同。你还可以根据自己的需要，对上述绘图逻辑进行调整，例如每个原始网格内部色彩独立映射等。
以上就是本文的全部内容，欢迎在评论区与我进行讨论~
转载地址：http://sxcuz.baihongyu.com/