本文對應程式碼已上傳至我的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中的投影座標系管理，敬請期待