> 本文完整程式碼及資料已上傳至我的`Github`倉庫[https://github.com/CNFeffery/FefferyViz](https://github.com/CNFeffery/FefferyViz) # 1 簡介 　　**OD資料**是交通、城市規劃以及GIS等領域常見的一類資料，特點是每一條資料都記錄了一次OD（*O*即*Origin*，*D*即*Destination*）行為的起點與終點座標資訊。 　　而針對**OD資料**常見的視覺化表達方式為弧線圖，譬如圖1所示的例子，就針對紐約曼哈頓等區域的某時間段**Uber**打車記錄上下車點資料進行展示： 圖1 　　但這種傳統的表達方式侷限很明顯：當OD記錄數量眾多時，因為不同線之間的彼此堆疊，導致很多區域之間的OD模式被遮蓋而難以被讀出。 　　而前一段時間我在觀看一場學術直播的過程中，注意到一種特別的表達區域間OD資料的方式，原始文獻比較老（ https://openaccess.city.ac.uk/id/eprint/537/1/wood_visualization_2010.pdf ）發表於2010年，其思想是通過對研究區域進行網格化劃分，再將整個區域的原始網格對映到每個單一網格中：

圖2 　　譬如圖2左圖中從座標記為$(E, 5)$的網格出發，到達記為$(A, 2)$的網格的所有OD資料記錄，可以在右圖中對應左圖$(E, 5)$位置的大網格中，劃分出的對應$(A, 2)$相對位置的小網格中進行記錄。 　　通過這樣的方式，原始文獻將圖3所示原始OD線圖轉換為圖4： 圖3 圖4 　　使得我們可以非常清楚地觀察到每個網格區域對其他網格區域的OD模式，而本文就將利用`Python`，在圖1對應的**Uber**上下車點分佈資料的基礎上，實踐這種表達OD資料的特別方式。 # 2 模仿過程 ## 2.1 過程分解 　　首先我們需要梳理一下整體的邏輯，先來看看原始的資料：

圖5 　　可以看到，原始資料中我們在本文真正用得到欄位為上車點經緯度`pickup_longitude`與`pickup_latitude`，以及下車點經緯度`dropoff_longitude`與`dropoff_latitude`。 　　我的思路是首先對所有經緯度點進行去重，接著儲存為`GeoDataFrame`並統一座標參考系為**Web墨卡託**也就是`EPSG:3857`： ```Python 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() ```

圖6 　　接下來我們來為研究區域建立網格面向量資料，思路是利用`numpy`先創建出x和y方向上的等間距座標，譬如我們這裡建立5行5列： ```Python 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資訊以方便之後的分析過程。 ```Python # 生成全部交叉線座標資訊 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 ``` 　　上面的建立網格的方法非常實用，愛學習的朋友的可以仔細看懂之後記錄下來。 　　我們來簡單看看創建出的網格是什麼樣子的，配合`contextily`新增上線上底圖： ```Python import matplotlib.pyplot as plt import contextily as ctx fig, ax = plt.subplots(figsize=(4, 4), dpi=200) ax = manhattan_grids.plot(facecolor='none', edgecolor='black', ax=ax) # 標註每個網格的id for row in manhattan_grids.itertuples(): centroid = row.geometry.centroid ax.text(centroid.x, centroid.y, row.id, ha='center', va='center') # 關閉座標軸 ax.axis('off') # 新增carto的素色底圖 ctx.add_basemap(ax, source='https://d.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png', zoom=12) fig.savefig('圖7.png', dpi=300, bbox_inches='tight', pad_inches=0) ``` 圖7 　　創建出的網格效果不錯~接下來就到了最關鍵的地方，我們需要計算出在每個原始網格內部上車的全部OD記錄，在整個區域中各個網格內的下車點分佈情況： 　　首先我們以某個網格為例，介紹如何為其關聯上車點、下車點以資訊，並利用簡單的仿射變換得到鑲嵌在其內部的小網格。 　　以id=21的網格為例，對應著肯尼迪國際機場的區域，首先我們利用id對應的從`manhattan_grids`表中提取的網格面數據，基於空間連線來與`od_points`表進行關聯，從而匹配到目標網格內對應原始od資訊表中的所有上車點記錄； 　　接著根據這些記錄對應的下車點資訊與`od_points`表進行匹配，從而得到所有下車點向量資訊，然後再次利用空間連線，得到所需的網格下車點分佈結果： ```Python i = 21 # 對應肯尼迪國際機場的網格 # 計算得到所有網格整體的重心座標 center_grid = (manhattan_grids.unary_union.centroid.x, manhattan_grids.unary_union.centroid.y) # 提取對應下車點座標 dropoff = ( # 利用空間連線，提取目標網格中包含到的所有座標點 gpd .sjoin(manhattan_grids.loc[i:i, :], right_df=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']] # 根據匹配到的下車點座標 # 與od_points表進行連線 # 找到對應下車點的向量資訊 .merge(od_points, left_on=['dropoff_longitude', 'dropoff_latitude'], right_on=['lng', 'lat'])[['geometry']] ) # 提取上一步得到的下車座標點在各個網格中的分佈資料 grid_distrib = ( # 利用空間連線匹配網格與下車座標點 gpd .sjoin(manhattan_grids, # 轉換為同一座標參考系的GeoDataFrame gpd.GeoDataFrame(dropoff, crs='EPSG:3857'), op='contains') # 根據網格id進行分組計數 .groupby('id', as_index=False) .agg({'index_right': 'count'}) .rename(columns={'index_right': '下車記錄數'}) ) grid_distrib.head() ``` 圖8 　　接著我們將上述的統計結果按照id列與原始網格表進行關聯，並利用仿射變換得到整體網格向目標網格內部的縮小鑲嵌結果（思路是首先將原始網格整體移動到與目標網格重心重合，接著按照x和y方向上的比例進行縮小），為了方便之後繪圖示記出目標網格對應的鑲嵌小網格位置，最後還需新增`是否為目標網格`列資訊： ```Python # 利用基本的仿射變換得到原始網格向對應目標網格的嵌入變換 # 獲取當前目標網格的重心座標 center_child_grid = (manhattan_grids.at[i, 'geometry'].centroid.x, manhattan_grids.at[i, 'geometry'].centroid.y) # 利用仿射變換得到整體網格在目標網格中的鑲嵌 draw_gdf = ( manhattan_grids # 基於原始的網格向量來更新放縮後的網格向量 .assign(geometry=manhattan_grids # 第一步：將原始網格的重心平移到目標網格的重心上 .translate(center_child_grid[0]-center_grid[0], center_child_grid[1]-center_grid[1]) # 第二步：以目標網格的重心為縮放中心，進行 .scale(xfact=1 / 5, yfact=1 / 5, origin=(manhattan_grids.at[i, 'geometry'].centroid.x, manhattan_grids.at[i, 'geometry'].centroid.y))) .merge(grid_distrib, on='id', how='left') .assign(是否為目標網格=0) ) draw_gdf.loc[draw_gdf.id == i, '是否為目標網格'] = 1 draw_gdf.head() ``` 圖9 　　經過這一系列操作，我們就得到了id為21的網格下車點分佈結果，將上述過程利用迴圈推廣到每個網格，並將最後的計算結果合併為一張`GeoDataFrame`，即表`draw_base`。 ## 2.2 繪製圖像 　　最終我們對`draw_base`表進行視覺化，這裡為了顯示更加自然，對`下車記錄`進行了**對數化**+**自然間斷**處理： ```Python %matplotlib inline fig, ax = plt.subplots(figsize=(12, 12)) # 繪製每個鑲嵌小網格的輪廓 ax = ( draw_base .plot(facecolor='none', edgecolor='lightgrey', ax=ax, linewidth=0.3) ) # 繪製每個鑲嵌小網格的下車記錄數熱力分佈 ax = ( draw_base .assign(下車記錄數=np.log(draw_base.下車記錄數)) .plot(column='下車記錄數', scheme='NaturalBreaks', k=5, cmap='YlOrRd', ax=ax, alpha=0.7) ) # 繪製原始網格的框架 ax = manhattan_grids.plot(ax=ax, facecolor='none', edgecolor='black', linewidth=0.8) # 在每個原始網格中標記出對應位置的鑲嵌小網格 ax = ( draw_base .query('是否為目標網格 == 1') .plot(facecolor='none', edgecolor='black', linestyle='--', ax=ax) ) # 設定繪圖區域範圍 minx, miny, maxx, maxy = manhattan_grids.total_bounds ax.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) ``` 圖10 　　通過這種表達方式，我們可以很明顯地看出不同區域相對其他區域出行模式的不同，你還可以根據自己的需要，對上述繪圖邏輯進行調整，譬如每個原始網格內部色彩獨立對映等。 --- 　　以上就是本文的全部內容，歡迎在評論區與我進