1 簡介

「OD數(shù)據(jù)」是交通、城市規(guī)劃以及GIS等領(lǐng)域常見的一類數(shù)據(jù)，特點是每一條數(shù)據(jù)都記錄了一次OD（O即Origin，D即Destination）行為的起點與終點坐標(biāo)信息。

而針對「OD數(shù)據(jù)」常見的可視化表達方式為弧線圖，譬如圖1所示的例子，就針對紐約曼哈頓等區(qū)域的某時間段「Uber」打車記錄上下車點數(shù)據(jù)進行展示：

圖1

但這種傳統(tǒng)的表達方式局限很明顯：當(dāng)OD記錄數(shù)量眾多時，因為不同線之間的彼此堆疊，導(dǎo)致很多區(qū)域之間的OD模式被遮蓋而難以被讀出。

而前一段時間我在觀看一場學(xué)術(shù)直播的過程中，注意到一種特別的表達區(qū)域間OD數(shù)據(jù)的方式，原始文獻比較老（ https://openaccess.city.ac.uk/id/eprint/537/1/wood_visualization_2010.pdf ）發(fā)表于2010年，其思想是通過對研究區(qū)域進行網(wǎng)格化劃分，再將整個區(qū)域的原始網(wǎng)格映射到每個單一網(wǎng)格中：

圖2

譬如圖2左圖中從坐標(biāo)記為的網(wǎng)格出發(fā)，到達記為的網(wǎng)格的所有OD數(shù)據(jù)記錄，可以在右圖中對應(yīng)左圖位置的大網(wǎng)格中，劃分出的對應(yīng)相對位置的小網(wǎng)格中進行記錄。

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

圖3

圖4

使得我們可以非常清楚地觀察到每個網(wǎng)格區(qū)域?qū)ζ渌W(wǎng)格區(qū)域的OD模式，而本文就將利用Python，在圖1對應(yīng)的「Uber」上下車點分布數(shù)據(jù)的基礎(chǔ)上，實踐這種表達OD數(shù)據(jù)的特別方式。

2 模仿過程

2.1 過程分解

首先我們需要梳理一下整體的邏輯，先來看看原始的數(shù)據(jù)：

圖5

可以看到，原始數(shù)據(jù)中我們在本文真正用得到字段為上車點經(jīng)緯度pickup_longitude與pickup_latitude，以及下車點經(jīng)緯度dropoff_longitude與dropoff_latitude。

我的思路是首先對所有經(jīng)緯度點進行去重，接著保存為GeoDataFrame并統(tǒng)一坐標(biāo)參考系為「Web墨卡托」也就是EPSG:3857：

from?shapely.geometry?import?Point
import?geopandas?as?gpd

od_points?=?\
(
????#?首先合并所有的經(jīng)緯度信息
????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'})])
????#?對經(jīng)緯度進行去重
????.drop_duplicates()
)

#?基于經(jīng)緯度信息為od_points添加矢量信息列
od_points['geometry']?=?(
????od_points
????.apply(lambda?row:?Point(row['lng'],?row['lat']),?axis=1)
)

#?轉(zhuǎn)換為GeoDataFrame并統(tǒng)一坐標(biāo)到Web墨卡托
od_points?=?gpd.GeoDataFrame(od_points,?crs='EPSG:4326').to_crs('EPSG:3857')

od_points.head()

圖6

接下來我們來為研究區(qū)域創(chuàng)建網(wǎng)格面矢量數(shù)據(jù)，思路是利用numpy先創(chuàng)建出x和y方向上的等間距坐標(biāo)，譬如我們這里創(chuàng)建5行5列：

from?shapely.geometry?import?MultiLineString
from?shapely.ops?import?polygonize?#?用于將交叉線轉(zhuǎn)換為網(wǎng)格面

#?提取所有上下車坐標(biāo)點范圍的左下角及右上角坐標(biāo)信息
xmin,?ymin,?xmax,?ymax?=?od_points.total_bounds

#?創(chuàng)建x方向上的所有坐標(biāo)位置
x?=?np.linspace(xmin,?
????????????????xmax,
????????????????6)

#?創(chuàng)建y方向上的所有坐標(biāo)位置
y?=?np.linspace(ymin,?
????????????????ymax,
????????????????6)

再利用雙層列表推導(dǎo)配合MultiLineString生成彼此交叉的網(wǎng)格線，并利用shapely中提供的polygonize工具直接把交叉線轉(zhuǎn)換為MultiPolygon，再拆分每個單一網(wǎng)格并添加一一對應(yīng)的id信息以方便之后的分析過程。

#?生成全部交叉線坐標(biāo)信息
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]

#?創(chuàng)建網(wǎng)格
manhattan_grids?=?gpd.GeoDataFrame({
????'geometry':?list(polygonize(MultiLineString(hlines?+?vlines)))},?
????crs='EPSG:3857')

#?添加一一對應(yīng)得id信息
manhattan_grids['id']?=?manhattan_grids.index

上面的創(chuàng)建網(wǎng)格的方法非常實用，愛學(xué)習(xí)的朋友的可以仔細看懂之后記錄下來。

我們來簡單看看創(chuàng)建出的網(wǎng)格是什么樣子的，配合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)

#?標(biāo)注每個網(wǎng)格的id
for?row?in?manhattan_grids.itertuples():
????
????centroid?=?row.geometry.centroid
????ax.text(centroid.x,?centroid.y,?row.id,?ha='center',?va='center')

#?關(guān)閉坐標(biāo)軸
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

創(chuàng)建出的網(wǎng)格效果不錯~接下來就到了最關(guān)鍵的地方，我們需要計算出在每個原始網(wǎng)格內(nèi)部上車的全部OD記錄，在整個區(qū)域中各個網(wǎng)格內(nèi)的下車點分布情況：

首先我們以某個網(wǎng)格為例，介紹如何為其關(guān)聯(lián)上車點、下車點信息，并利用簡單的仿射變換得到鑲嵌在其內(nèi)部的小網(wǎng)格。

以id=21的網(wǎng)格為例，對應(yīng)著肯尼迪國際機場的區(qū)域，首先我們利用id對應(yīng)的從manhattan_grids表中提取的網(wǎng)格面數(shù)據(jù)，基于空間連接來與od_points表進行關(guān)聯(lián)，從而匹配到目標(biāo)網(wǎng)格內(nèi)對應(yīng)原始o(jì)d信息表中的所有上車點記錄；

接著根據(jù)這些記錄對應(yīng)的下車點信息與od_points表進行匹配，從而得到所有下車點矢量信息，然后再次利用空間連接，得到所需的網(wǎng)格下車點分布結(jié)果：

i?=?21?#?對應(yīng)肯尼迪國際機場的網(wǎng)格

#?計算得到所有網(wǎng)格整體的重心坐標(biāo)
center_grid?=?(manhattan_grids.unary_union.centroid.x,?
???????????????manhattan_grids.unary_union.centroid.y)

#?提取對應(yīng)下車點坐標(biāo)
dropoff?=?(
????#?利用空間連接，提取目標(biāo)網(wǎng)格中包含到的所有坐標(biāo)點
????gpd
????.sjoin(manhattan_grids.loc[i:i,?:],
???????????right_df=od_points,?
???????????op='contains')
????[['lng',?'lat',?'geometry']]
????#?利用提取到的坐標(biāo)點信息，關(guān)聯(lián)在目標(biāo)
????#?網(wǎng)格中上車的記錄對應(yīng)的下車點坐標(biāo)
????.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']]
????#?根據(jù)匹配到的下車點坐標(biāo)
????#?與od_points表進行連接
????#?找到對應(yīng)下車點的矢量信息
????.merge(od_points,
???????????left_on=['dropoff_longitude',?'dropoff_latitude'],
???????????right_on=['lng',?'lat'])[['geometry']]
)

#?提取上一步得到的下車坐標(biāo)點在各個網(wǎng)格中的分布數(shù)據(jù)
grid_distrib?=?(
????#?利用空間連接匹配網(wǎng)格與下車坐標(biāo)點
????gpd
????.sjoin(manhattan_grids,
???????????#?轉(zhuǎn)換為同一坐標(biāo)參考系的GeoDataFrame
???????????gpd.GeoDataFrame(dropoff,?crs='EPSG:3857'),
???????????op='contains')
????#?根據(jù)網(wǎng)格id進行分組計數(shù)
????.groupby('id',?as_index=False)
????.agg({'index_right':?'count'})
????.rename(columns={'index_right':?'下車記錄數(shù)'})
)

grid_distrib.head()

圖8

接著我們將上述的統(tǒng)計結(jié)果按照id列與原始網(wǎng)格表進行關(guān)聯(lián)，并利用仿射變換得到整體網(wǎng)格向目標(biāo)網(wǎng)格內(nèi)部的縮小鑲嵌結(jié)果（思路是首先將原始網(wǎng)格整體移動到與目標(biāo)網(wǎng)格重心重合，接著按照x和y方向上的比例進行縮小），為了方便之后繪圖標(biāo)記出目標(biāo)網(wǎng)格對應(yīng)的鑲嵌小網(wǎng)格位置，最后還需添加是否為目標(biāo)網(wǎng)格列信息：

#?利用基本的仿射變換得到原始網(wǎng)格向?qū)?yīng)目標(biāo)網(wǎng)格的嵌入變換

#?獲取當(dāng)前目標(biāo)網(wǎng)格的重心坐標(biāo)
center_child_grid?=?(manhattan_grids.at[i,?'geometry'].centroid.x,?
?????????????????????manhattan_grids.at[i,?'geometry'].centroid.y)

#?利用仿射變換得到整體網(wǎng)格在目標(biāo)網(wǎng)格中的鑲嵌
draw_gdf?=?(
????manhattan_grids
????#?基于原始的網(wǎng)格矢量來更新放縮后的網(wǎng)格矢量
????.assign(geometry=manhattan_grids
????????????#?第一步：將原始網(wǎng)格的重心平移到目標(biāo)網(wǎng)格的重心上
????????????.translate(center_child_grid[0]-center_grid[0],?
???????????????????????center_child_grid[1]-center_grid[1])
????????????#?第二步：以目標(biāo)網(wǎng)格的重心為縮放中心，進行
????????????.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(是否為目標(biāo)網(wǎng)格=0)
)

draw_gdf.loc[draw_gdf.id?==?i,?'是否為目標(biāo)網(wǎng)格']?=?1
draw_gdf.head()

圖9

經(jīng)過這一系列操作，我們就得到了id為21的網(wǎng)格下車點分布結(jié)果，將上述過程利用循環(huán)推廣到每個網(wǎng)格，并將最后的計算結(jié)果合并為一張GeoDataFrame，即表draw_base。

2.2 繪制圖像

最終我們對draw_base表進行可視化，這里為了顯示更加自然，對下車記錄進行了「對數(shù)化」+「自然間斷」處理：

%matplotlib?inline
fig,?ax?=?plt.subplots(figsize=(12,?12))

#?繪制每個鑲嵌小網(wǎng)格的輪廓
ax?=?(
????draw_base
????.plot(facecolor='none',?edgecolor='lightgrey',?ax=ax,
??????????linewidth=0.3)
)

#?繪制每個鑲嵌小網(wǎng)格的下車記錄數(shù)熱力分布
ax?=?(
????draw_base
????.assign(下車記錄數(shù)=np.log(draw_base.下車記錄數(shù)))
????.plot(column='下車記錄數(shù)',?scheme='NaturalBreaks',?
??????????k=5,?cmap='YlOrRd',?ax=ax,?alpha=0.7)
)

#?繪制原始網(wǎng)格的框架
ax?=?manhattan_grids.plot(ax=ax,?facecolor='none',?edgecolor='black',
??????????????????????????linewidth=0.8)

#?在每個原始網(wǎng)格中標(biāo)記出對應(yīng)位置的鑲嵌小網(wǎng)格
ax?=?(
????draw_base
????.query('是否為目標(biāo)網(wǎng)格?==?1')
????.plot(facecolor='none',?edgecolor='black',?
??????????linestyle='--',?ax=ax)
)

#?設(shè)置繪圖區(qū)域范圍
minx,?miny,?maxx,?maxy?=?manhattan_grids.total_bounds
ax.set_xlim(minx,?maxx)
ax.set_ylim(miny,?maxy)

#?關(guān)閉坐標(biāo)軸
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)