（數據科學學習手札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中的投影座標系管理，敬請期待。