（數據科學學習手札74）基於geopandas的空間數據分析——數據結構篇

本文對應代碼已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

　　geopandas是建立在GEOS、GDAL、PROJ等開源地理空間計算相關框架之上的，類似pandas語法風格的空間數據分析Python庫，其目標是盡可能地簡化Python中的地理空間數據處理，減少對Arcgis、PostGIS等工具的依賴，使得處理地理空間數據變得更加高效簡潔，打造純Python式的空間數據處理工作流。本系列文章就將圍繞geopandas及其使用過程中涉及到的其他包進行系統性的介紹說明，每一篇將盡可能全面具體地介紹geopandas對應方面的知識，計划涵蓋geopandas的數據結構、投影坐標系管理、文件IO、基礎地圖制作、集合操作、空間連接與聚合。

　　作為基於geopandas的空間數據分析系列文章的第一篇，通過本文你將會學習到geopandas中的數據結構。

　　geopandas的安裝和使用需要若干依賴包，如果不事先妥善安裝好這些依賴包而直接使用pip install geopandas或conda install geopandas可能會引發依賴包相關錯誤導致安裝失敗，官方文檔中的推薦安裝方式為：

conda install --channel conda-forge geopandas

　　conda-forge是一個社區項目，在conda的基礎上提供了更廣泛更豐富的軟件資源包，通過它我們可以自動下載安裝好所有geopandas的必要依賴包而無需手動繁瑣地去安裝它們。在完成安裝后，下面我們開始對geopandas的系統性學習之旅。

2 數據結構

　　geopandas作為pandas向地理分析計算方面的延拓，基礎的數據結構延續了Series和DataFrame的特點，創造出GeoSeries與GeoDataFrame兩種基礎數據結構：

2.1 GeoSeries

2.1.1 GeoSeries中的基礎幾何對象

　　與Series相似，GeoSeries用來表示一維向量，只不過這里的向量每個位置上的元素都表示着一個shapely中的幾何對象，有如下幾種類型：

• Points

　　對應shapely.geometry中的Point，用於表示單個點，下面我們創建一個由若干Point對象組成的GeoSeries並像Series一樣定義索引：

from shapely import geometry
import geopandas as gpd

# 創建存放Point對象的GeoSeries
# 這里shapely.geometry.Point(x, y)用於創建單個點對象
gpd.GeoSeries([geometry.Point(0, 0),
               geometry.Point(0, 1),
               geometry.Point(1, 1),
               geometry.Point(1, 0)],
              index=['a', 'b', 'c', 'd'])

圖1

　　可以看到創建出的GeoSeries數據類型為geometry，即幾何對象。

• MultiPoint

　　對應shapely中的MultiPoint，用於表示多個點的集合，下面我們創建一個由若干MultiPoint對象組成的GeoSeries：

# 創建存放MultiPoint對象的GeoSeries
# 這里shapely.geometry.MultiPoint([(x1, y1), (x2, y2), ...])用於創建多點集合
gpd.GeoSeries([geometry.MultiPoint([(0, 1), (1, 0)]),
               geometry.MultiPoint([(0, 0), (1, 1)])],
              index=['a', 'b'])

圖2

　　在jupyter notebook或jupyter lab中可以圖像的形式直接顯示GeoSeries中的單個元素：

圖3

• LineString

　　對應shapely中的LineString，用於表示由多個點按順序連接而成的線，下面我們創建一個由若干LineString對象組成的GeoSeries：

# 創建存放LineString對象的GeoSeries
# 這里shapely.geometry.LineString([(x1, y1), (x2, y2), ...])用於創建多點按順序連接而成的線段
gpd.GeoSeries([geometry.LineString([(0, 0), (1, 1), (1, 0)]),
               geometry.LineString([(0, 0), (0, 1), (-1, 0)])],
              index=['a', 'b'])

圖4

　　同樣地，直接顯示第一個元素：

圖5

• MultiLineString

　　對應shapely中的MultiLineString，用於表示多條線段的集合，下面我們創建一個由若干MultiLineString對象組成的GeoSeries：

# 創建存放MultiLineString對象的GeoSeries
# 這里shapely.geometry.MultiLineString([LineString1, LineString2])用於創建多條線段的集合
gpd.GeoSeries([geometry.MultiLineString([[(0, 0), (1, 1), (1, 0)],
                                        [(-0.5, 0), (0, 1), (-1, 0)]])],
              index=['a'])

圖6

　　同樣地，直接顯示第一個元素：

圖7

• Polygon(無孔)

　　geopandas中的Polygon對應shapely中的Polygon，用於表示面，根據內部有無孔洞可繼續細分。下面我們創建一個由無孔Polygon對象組成的GeoSeries：

# 創建存放無孔Polygon對象的GeoSeries
# 這里shapely.geometry.Polygon([(x1, y1), (x2, y2),...])用於創建無孔面
gpd.GeoSeries([geometry.Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])],
              index=['a'])

圖8

　　同樣地，直接顯示第一個元素：

圖9

• Polygon(有孔)

　　區分於上文中的無孔Polygon，下面我們創建一個由有孔Polygon對象組成的GeoSeries：

# 創建存放有孔Polygon對象的GeoSeries
# 這里shapely.geometry.Polygon(polygonExteriors, interiorCoords)用於創建有孔面
# 其中polygonExteriors用於定義整個有孔Polygon的外圍，是一個無孔的多邊形
# interiorCoords是用於定義內部每個孔洞（本質上是獨立的多邊形）的序列
gpd.GeoSeries([geometry.Polygon([(0,0),(10,0),(10,10),(0,10)], 
                                [((1,3),(5,3),(5,1),(1,1)), 
                                 ((9,9),(9,8),(8,8),(8,9))])])

　　同樣地，直接顯示第一個元素：

圖10

• MultiPolygon

　　對應shapely中的MultiPolygon，用於表示多個面的集合，下面我們創建一個由MultiPolygon對象組成的GeoSeries：

# 創建存放MultiPolygon對象的GeoSeries
# 這里shapely.geometry.MultiPolygon([Polygon1, Polygon2])用於創建多個面的集合
gpd.GeoSeries([geometry.MultiPolygon([geometry.Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
                                      geometry.Polygon([(2, 2), (2, 3), (3, 3), (3, 2), (2, 2)])])],
              index=['a'])

圖11

　　顯示第一個元素：

圖12

• LinearRing

　　LinearRing對應shapely.geometry中的LinearRing，是一種特殊的幾何對象，可以理解為閉合的線或無孔多邊形的邊框，創建時傳入數據的格式與Polygon相同，下面我們創建一個由LinearRing對象組成的GeoSeries：

# 創建存放LinearRing對象的GeoSeries
# 這里shapely.geometry.LinearRing([(x1, y1), (x2, y2),...])用於創建LinearRing
gpd.GeoSeries([geometry.LinearRing([(0, 0), (0, 1), (1, 1), (1, 0)])],
              index=['a'])

圖13

　　顯示第一個元素，可以看出LinearRing就是無孔多邊形的邊框線：

圖14

　　在同一個GeoSeries可以混合上述類型中的多種幾何對象，這意味着點線面在概念上相異的幾何對象可以共存於同一份數據中

2.1.2 GeoSeries常用屬性

　　類似pandas中的Series，GeoSeries在被創建完成之后也擁有很多實用的地理屬性，下面對其中較為常用的進行列舉：

• area

　　area屬性返回與GeoSeries中每個元素一一對應的面積值（這里的面積單位和下文涉及的長度單位取決於投影坐標系，之后關於geopandas投影坐標系管理的文章將會詳細介紹，這里僅做演示）：

# 創建混合點線面的GeoSeries，這里第5個有孔多邊形內部空洞創建時使用[::-1]顛倒順序
# 是因為GeoSeries.plot()方法繪制有孔多邊形的一個bug，即外部邊框與內部孔洞創建時坐標
# 方向同為順時針或順時針時內部孔洞會自動被填充，如果你對這個bug感興趣，可以前往
# https://github.com/geopandas/geopandas/issues/951查看細節
s = gpd.GeoSeries([geometry.Polygon([(0, 0), (0.5, 0.5), (1, 0), (0.5, -0.5)]),
                   geometry.Polygon([(1, 1), (1.5, 1.5), (2, 1), (1.5, -1.5)]),
                   geometry.Point(3, 3),
                   geometry.LineString([(2, 2), (0, 3)]),
                   geometry.Polygon([(4, 4), (8, 4), (8, 8), (4, 8)], 
                                [[(5, 5), (7, 5), (7, 7), (5, 7)][::-1]])])

# 在jupyter中開啟matplotlib交互式繪圖模式
%matplotlib widget
s.plot() # 對s進行簡單的可視化

圖15

　　可以看到，s中包含了多種幾何對象，下面直接得到s的面積：

圖16 計算GeoSeries面積

• bounds

　　bounds屬性返回每個幾何對象所在box左下角、右上角的坐標信息：

圖17

• length

　　length屬性返回每個幾何對象邊長：

圖18

• geom_type

　　geom_type返回每個幾何對象類型：

圖19

• exterior與interiors

　　對於多邊形對象，exterior返回LinearRing格式的外邊框線，對於有孔多邊形，interiors返回所有內部孔洞LinearRing格式邊框線集合：

圖20

• is_valid

　　在shapely中涉及到很多拓撲計算操作時，對幾何對象的合法性有要求，譬如定義多邊形時坐標按順序連線時穿過了之前定義的邊就屬於非法，因為geopandas對矢量對象的計算依賴於shapely，於是引進了屬性用於判斷每個幾何對象是否合法，下面我們創建兩個形狀相同的多邊形，其中一個滿足上述所說的非法情況，另一個由兩個多邊形拼接而成：

s_ = gpd.GeoSeries([geometry.Polygon([(4, 0), (6, 1), (4, 1), (6, 0)]),
               geometry.MultiPolygon([geometry.Polygon([(4, 0), (5, 0.5), (6, 0)]), 
                                      geometry.Polygon([(5, 0.5), (6, 1), (4, 1)])])])

　　從形狀上看兩者相同：

圖21

　　下面我們嘗試用shapely中的intersection方法來取得這兩個幾何對象的相交部分，出現了拓撲邏輯錯誤：

圖22

　　查看s_.is_valid，可以看出第一個自相交的多邊形非法：

圖23

• boundary

　　boundary返回每個幾何對象的低維簡化表示（點對象無具體的更低維簡化，故無返回值）：

圖24

• centroid

　　centroid返回每個幾何對象的重心（幾何中心）：

圖25

• convex_hull

　　convex_hull返回每個幾何對象的凸包，Polygon格式，即恰巧包含對應幾何對象的凸多邊形：

import numpy as np

# 利用獨立的正態分布隨機數創建兩個MultiPoint集合
s__ = gpd.GeoSeries([geometry.MultiPoint(np.random.normal(loc=0, scale=2, size=[10, 2]).tolist()),
               geometry.MultiPoint(np.random.normal(loc=5, scale=2, size=[10, 2]).tolist())])

ax = s__.plot(color='red') # 繪制s__
s__.convex_hull.plot(ax=ax, alpha=0.4) # 疊加繪制各自對應凸包，調低填充透明度以顯示更明顯

圖26

• envelope

　　envelope屬性返回對應幾何對象的box范圍，Polygon格式，即包含對應元素中所有點的最小矩形：

import numpy as np

# 創建兩團獨立的MultiPoint
s__ = gpd.GeoSeries([geometry.MultiPoint(np.random.normal(loc=0, scale=2, size=[10, 2]).tolist()),
               geometry.MultiPoint(np.random.normal(loc=5, scale=2, size=[10, 2]).tolist())])

ax = s__.plot(color='red') # 繪制s__
s__.envelope.plot(ax=ax, alpha=0.4) # 疊加繪制各自對應envelope，調低填充透明度以顯示更明顯

圖27

2.2 GeoDataFrame

2.2.1 GeoDataFrame基礎

　　顧名思義，geopandas中的GeoDataFrame是在pandas.DataFrame的基礎上，加入空間分析相關內容進行改造而成，其最大特點在於其在原有數據表格基礎上增加了一列GeoSeries使得其具有矢量性，所有對於GeoDataFrame施加的空間幾何操作也都作用在這列指定的幾何對象之上。下面我們舉個簡單的例子，基於不同均值和標准差的正態分布隨機數，創建GeoDataFrame來記錄這些信息：

contents = [(loc, 0.5) for loc in range(0, 10, 2)]
geo_df = gpd.GeoDataFrame(data=contents, 
                          geometry=[geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents],
                          columns=['均值', '標准差'])
geo_df

圖28

　　其中定義GeoDataFrame時作為每行所關聯幾何對象的GeoSeries需要通過geometry參數指定，而除了用上述的方式創建GeoDataFrame，先創建數據表，再添加矢量信息列亦可，這時幾何對象列的名稱可以自由設置，但一定要利用GeoDataFrame.set_geometry()方法將后添加的矢量列指定為矢量主列，因為每個GeoDataFrame若在定義之處沒有指定矢量列，后將無法進行與適量信息掛鈎的所有操作（GeoSeries所有屬性都可同樣作用於GeoDataFrame，因為所有空間操作實際上都直接作用於其矢量主列）：

• 添加矢量列但未定義

geo_df = gpd.GeoDataFrame(contents, columns=['均值', '標准差'])
geo_df['raw_points'] = [geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents]
# 嘗試查看矢量類型
geo_df.geom_type

圖29

　　這時所有直接針對GeoDataFrame的矢量相關操作都無法使用。

• 重新為GeoDataFrame指定矢量列

geo_df.set_geometry('raw_points').geom_type

　　這時相關操作可正常使用：

圖30

• 多個矢量列切換

　　通過前面的內容，我們知道了每個GeoDataFrame都有一個矢量主列，相關操作例如繪圖都基於此列，實際上GeoDataFrame允許表中存在多個矢量列，只要求任意時刻有且僅有1列為矢量主列即可，因此我們可以在一個GeoDataFrame中保存多列矢量，需要用到哪列時再進行切換即可，如下面的例子：

geo_df = gpd.GeoDataFrame(contents, columns=['均值', '標准差'])
geo_df['raw_points'] = [geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents]
geo_df.set_geometry('raw_points', inplace=True) # inplace=True表示對原數據進行更新

# 繪制第一圖層
ax = geo_df.plot(color='red')
geo_df['convex_hull'] = geo_df.convex_hull
# 切換矢量主列
geo_df.set_geometry('convex_hull', inplace=True)

# 繪制第二圖層
geo_df.plot(ax=ax, color='blue', alpha=0.4)

圖31

2.2.2 GeoDataFrame數據索引

　　作為pandas.DataFrame的延伸，GeoDataFrame同樣支持pandas.DataFrame中的.loc以及.iloc對數據在行、列尺度上進行索引和篩選，這里我們以geopandas自帶的世界地圖數據為例：

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot()

圖32 geopandas自帶世界地圖

　　查看其表格內容：

圖33

　　使用.loc+條件篩選選擇數據：

圖34

　　使用.iloc選擇數據：

圖35

　　而除了這些常規的數據索引方式之外，geopandas為GeoDataFrame添加了.cx索引方式，可以傳入所需的空間范圍，用於索引與傳入范圍相交的對應數據：

# 選擇與東經80度-110度，北緯0度-30度范圍相交的幾何對象
part_world = world.cx[80:110, 0:30]

# 繪制第一圖層：世界地圖
ax = world.plot(alpha=0.05)
# 繪制第二圖層：.cx所選擇的地區
ax = part_world.plot(ax=ax, alpha=0.6)
# 繪制第三圖層：.cx條件示意圖
ax = gpd.GeoSeries([geometry.box(minx=80, miny=0, maxx=110, maxy=30).boundary])\
    .plot(ax=ax, color='red')

　　示意圖如下：

圖36

　　放大到所選區域，可以看出正如前面所說，通過.cx，所有與指定空間范圍有重疊的對象都被選擇：

圖37

　　以上就是本文的全部內容，如有筆誤望指出，系列文章下一篇將詳細介紹geopandas中的投影坐標系管理，敬請期待。