1 簡(jiǎn)介

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

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

圖1

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

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

圖2

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

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

圖3

圖4

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

2 模仿過程

2.1 過程分解

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

圖5

可以看到，原始數(shù)據(jù)中我們?cè)诒疚恼嬲玫玫阶侄螢樯宪圏c(diǎn)經(jīng)緯度pickup_longitude與pickup_latitude，以及下車點(diǎn)經(jīng)緯度dropoff_longitude與dropoff_latitude。

我的思路是首先對(duì)所有經(jīng)緯度點(diǎn)進(jìn)行去重，接著保存為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'})])
????#?對(duì)經(jīng)緯度進(jìn)行去重
????.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)點(diǎn)范圍的左下角及右上角坐標(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，再拆分每個(gè)單一網(wǎng)格并添加一一對(duì)應(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')

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

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

我們來簡(jiǎn)單看看創(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)注每個(gè)網(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)格效果不錯(cuò)~接下來就到了最關(guān)鍵的地方，我們需要計(jì)算出在每個(gè)原始網(wǎng)格內(nèi)部上車的全部OD記錄，在整個(gè)區(qū)域中各個(gè)網(wǎng)格內(nèi)的下車點(diǎn)分布情況：

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

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

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

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

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

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

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

grid_distrib.head()

圖8

接著我們將上述的統(tǒng)計(jì)結(jié)果按照id列與原始網(wǎng)格表進(jìn)行關(guān)聯(lián)，并利用仿射變換得到整體網(wǎng)格向目標(biāo)網(wǎng)格內(nèi)部的縮小鑲嵌結(jié)果（思路是首先將原始網(wǎng)格整體移動(dòng)到與目標(biāo)網(wǎng)格重心重合，接著按照x和y方向上的比例進(jìn)行縮?。?，為了方便之后繪圖標(biāo)記出目標(biāo)網(wǎng)格對(duì)應(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)格的重心為縮放中心，進(jìn)行
????????????.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)格下車點(diǎn)分布結(jié)果，將上述過程利用循環(huán)推廣到每個(gè)網(wǎng)格，并將最后的計(jì)算結(jié)果合并為一張GeoDataFrame，即表draw_base。

2.2 繪制圖像

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

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

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

#?繪制每個(gè)鑲嵌小網(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)

#?在每個(gè)原始網(wǎng)格中標(biāo)記出對(duì)應(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)