（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式

費弗裡發表於2020-10-18
視覺化

本文完整程式碼及資料已上傳至我的Github倉庫https://github.com/CNFeffery/FefferyViz

1 簡介

　　OD資料是交通、城市規劃以及GIS等領域常見的一類資料，特點是每一條資料都記錄了一次OD（O即Origin，D即Destination）行為的起點與終點座標資訊。

　　而針對OD資料常見的視覺化表達方式為弧線圖，譬如圖1所示的例子，就針對紐約曼哈頓等區域的某時間段Uber叫車記錄上下車點資料進行展示：

（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖1

　　但這種傳統的表達方式侷限很明顯：當OD記錄數量眾多時，因為不同線之間的彼此堆疊，導致很多區域之間的OD模式被遮蓋而難以被讀出。

　　而前一段時間我在觀看一場學術直播的過程中，注意到一種特別的表達區域間OD資料的方式，原始文獻比較老（ https://openaccess.city.ac.uk/id/eprint/537/1/wood_visualization_2010.pdf ）發表於2010年，其思想是通過對研究區域進行網格化劃分，再將整個區域的原始網格對映到每個單一網格中：

（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖2

　　譬如圖2左圖中從座標記為\((E, 5)\)的網格出發，到達記為\((A, 2)\)的網格的所有OD資料記錄，可以在右圖中對應左圖\((E, 5)\)位置的大網格中，劃分出的對應\((A, 2)\)相對位置的小網格中進行記錄。

　　通過這樣的方式，原始文獻將圖3所示原始OD線圖轉換為圖4：

（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖3
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖4

　　使得我們可以非常清楚地觀察到每個網格區域對其他網格區域的OD模式，而本文就將利用Python，在圖1對應的Uber上下車點分佈資料的基礎上，實踐這種表達OD資料的特別方式。

2 模仿過程

2.1 過程分解

　　首先我們需要梳理一下整體的邏輯，先來看看原始的資料：

（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖5

　　可以看到，原始資料中我們在本文真正用得到欄位為上車點經緯度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()
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖6

　　接下來我們來為研究區域建立網格面向量資料，思路是利用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

　　上面的建立網格的方法非常實用，愛學習的朋友的可以仔細看懂之後記錄下來。

　　我們來簡單看看建立出的網格是什麼樣子的，配合contextily新增上線上底圖：

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)
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖7

　　建立出的網格效果不錯~接下來就到了最關鍵的地方，我們需要計算出在每個原始網格內部上車的全部OD記錄，在整個區域中各個網格內的下車點分佈情況：

　　首先我們以某個網格為例，介紹如何為其關聯上車點、下車點以資訊，並利用簡單的仿射變換得到鑲嵌在其內部的小網格。

　　以id=21的網格為例，對應著肯尼迪國際機場的區域，首先我們利用id對應的從manhattan_grids表中提取的網格面資料，基於空間連線來與od_points表進行關聯，從而匹配到目標網格內對應原始od資訊表中的所有上車點記錄；

　　接著根據這些記錄對應的下車點資訊與od_points表進行匹配，從而得到所有下車點向量資訊，然後再次利用空間連線，得到所需的網格下車點分佈結果：

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()
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖8

　　接著我們將上述的統計結果按照id列與原始網格表進行關聯，並利用仿射變換得到整體網格向目標網格內部的縮小鑲嵌結果（思路是首先將原始網格整體移動到與目標網格重心重合，接著按照x和y方向上的比例進行縮小），為了方便之後繪圖示記出目標網格對應的鑲嵌小網格位置，最後還需新增是否為目標網格列資訊：

# 利用基本的仿射變換得到原始網格向對應目標網格的嵌入變換

# 獲取當前目標網格的重心座標
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()
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖9

　　經過這一系列操作，我們就得到了id為21的網格下車點分佈結果，將上述過程利用迴圈推廣到每個網格，並將最後的計算結果合併為一張GeoDataFrame，即表draw_base。

2.2 繪製影像

　　最終我們對draw_base表進行視覺化，這裡為了顯示更加自然，對下車記錄進行了對數化+自然間斷處理：

%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)
（在模仿中精進資料視覺化03）OD資料的特殊視覺化方式
圖10

　　通過這種表達方式，我們可以很明顯地看出不同區域相對其他區域出行模式的不同，你還可以根據自己的需要，對上述繪圖邏輯進行調整，譬如每個原始網格內部色彩獨立對映等。

　　以上就是本文的全部內容，歡迎在評論區與我進行討論~

相關文章