OD数据是交通、城市规划以及GIS等领域常见的一类数据，特点是每一条数据都记录了一次OD（O即Origin，D即Destination）行为的起点与终点坐标信息。
而针对OD数据常见的可视化表达方式为弧线图，譬如图1所示的例子，就针对纽约曼哈顿等区域的某时间段Uber打车记录上下车点数据进行展示：
但这种传统的表达方式局限很明显：当OD记录数量众多时，因为不同线之间的彼此堆叠，导致很多区域之间的OD模式被遮盖而难以被读出。
而前一段时间我在观看一场学术直播的过程中，注意到一种特别的表达区域间OD数据的方式，原始文献比较老（ https://openaccess.city.ac.uk/id/eprint/537/1/wood_visualization_2010.pdf ）发表于2010年，其思想是通过对研究区域进行网格化划分，再将整个区域的原始网格映射到每个单一网格中：
譬如图2左图中从坐标记为\((E, 5)\)的网格出发，到达记为\((A, 2)\)的网格的所有OD数据记录，可以在右图中对应左图\((E, 5)\)位置的大网格中，划分出的对应\((A, 2)\)相对位置的小网格中进行记录。
通过这样的方式，原始文献将图3所示原始OD线图转换为图4：
使得我们可以非常清楚地观察到每个网格区域对其他网格区域的OD模式，而本文就将利用Python，在图1对应的Uber上下车点分布数据的基础上，实践这种表达OD数据的特别方式。
2 模仿过程 2.1 过程分解首先我们需要梳理一下整体的逻辑，先来看看原始的数据：
可以看到，原始数据中我们在本文真正用得到字段为上车点经纬度pickup_longitude与pickup_latitude，以及下车点经纬度dropoff_longitude与dropoff_latitude。
我的思路是首先对所有经纬度点进行去重，接着保存为GeoDataFrame并统一坐标参考系为Web墨卡托也就是EPSG:3857：
from shapely.geometry import Point import geopandas as gpd od_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'})]) # 对经纬度进行去重 .drop_duplicates() ) # 基于经纬度信息为od_points添加矢量信息列 od_points['geometry'] = ( od_points .apply(lambda row: Point(row['lng'], row['lat']), axis=1) ) # 转换为GeoDataFrame并统一坐标到Web墨卡托 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 MultiLineString from shapely.ops import polygonize # 用于将交叉线转换为网格面 # 提取所有上下车坐标点范围的左下角及右上角坐标信息 xmin, ymin, xmax, ymax = od_points.total_bounds # 创建x方向上的所有坐标位置 x = np.linspace(xmin, xmax, 6) # 创建y方向上的所有坐标位置 y = np.linspace(ymin, ymax, 6)再利用双层列表推导配合MultiLineString生成彼此交叉的网格线，并利用shapely中提供的polygonize工具直接把交叉线转换为MultiPolygon，再拆分每个单一网格并添加一一对应的id信息以方便之后的分析过程。
# 生成全部交叉线坐标信息 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 = gpd.GeoDataFrame({ 'geometry': list(polygonize(MultiLineString(hlines + vlines)))}, crs='EPSG:3857') # 添加一一对应得id信息 manhattan_grids['id'] = manhattan_grids.index上面的创建网格的方法非常实用，爱学习的朋友的可以仔细看懂之后记录下来。