1. 程式人生 > >(資料科學學習手札77)基於geopandas的空間資料分析——檔案IO

(資料科學學習手札77)基於geopandas的空間資料分析——檔案IO

本文對應程式碼和資料已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

  在上一篇文章中我們對geopandas中的座標參考系有了較為深入的學習,而在日常空間資料分析工作中向量檔案的讀入和寫出,是至關重要的環節。

  作為基於geopandas的空間資料分析系列文章的第三篇,通過本文你將會學習到geopandas中的檔案IO。

2 檔案IO

2.1 向量檔案的讀入

  geopandasfiona作為操縱向量資料讀寫功能的後端,使用geopandas.read_file()讀取對應型別檔案,而在後端實際上是使用fiona.open

來讀入資料,即兩者引數是保持一致的,讀入的資料自動轉換為GeoDataFrame,下面是geopandas.read_file()主要引數:

filename:str型別,傳入檔案對應的路徑或url

layer:str型別,當要讀入的資料格式為地理資料庫.gdbQGIS中的.gpkg時,傳入對應圖層的名稱

  下面結合上述引數,來介紹一下使用geopandas.read_file()在不同情況下讀取常見格式向量資料的方法,使用到的示例資料為中國地圖,CRSEPSG:4326,本文使用到的所有資料都可以在文章開頭提及的Github倉庫對應本文路徑下找到:

圖1

2.1.1 shapefile

  作為非常常見的一種向量檔案格式,geopandasshapefile提供了很好的讀取和寫出支援,下面分為不同情況來介紹:

  • 完整的shapefile

  如圖2,這是一個完整的shapefile

圖2

  使用geopandas來讀取這種形式的shapefile很簡單:

import geopandas as gpd

data = gpd.read_file('geometry/china_provinces/china_provinces.shp')
print(data.crs) # 檢視資料對應的crs
data.head() # 檢視前5行
圖3
  • 缺少投影的shapefile

  當shapefile中缺失.prj檔案時,使用geopandas讀入後形成的GeoDataFrame會缺失crs屬性:
  

圖4

  如果已經知道資料對應的CRS,可以在讀入資料後補充上crs資訊以進行其他操作:

import pyproj

data.crs = pyproj.CRS.from_user_input('EPSG:4326')
data.crs
圖5
  • 直接讀取資料夾

  當資料夾下只有單個shapefile時,可以直接讀取該資料夾:

圖6
  • 讀取zip壓縮包中的檔案

  geopandas通過傳入特定語法格式的檔案路徑資訊,以支援直接讀取.zip格式壓縮包中的shapefile檔案,主要分為兩種情況。
  當檔案在壓縮包內的根目錄時,使用下面的語法規則來讀取資料:

zip://路徑/xxx.zip

  譬如我們要讀取圖7所示的壓縮包內檔案:

圖7

  按照對應的語法規則,讀取該型別資料方式如下:

圖8

  而當檔案在壓縮包內的資料夾中時,如圖9:
  

圖9

  使用下面的語法規則來讀取資料:

zip://路徑/xxx.zip!壓縮包內指定檔案路徑

  將上述語法運用到上述檔案:

圖10

2.1.2 gdb與gpkg

  對於Arcgis中的地理資料庫gdb,以及QGIS中的GeoPackage,要讀取其包含的向量資料,就要涉及到圖層的概念,對應geopandas.read_file()layer引數,只需要將gdbgpkg檔案路徑作為filename引數,再將對應的圖層名稱作為layer引數傳入:

  • gdb
data = gpd.read_file('geometry/china_provinces.gdb', 
                     layer='china_provinces')
print(data.crs) # 檢視資料對應的crs
data.head() # 檢視前5行
圖11
  • gpkg

  類似讀入gdb檔案:

data = gpd.read_file('geometry/china_provinces.gpkg', 
                     layer='china_provinces',
                     encoding='utf-8')
print(data.crs) # 檢視資料對應的crs
data.head() # 檢視前5行
圖12

2.1.3 GeoJSON

  作為web地圖中最常使用的向量資料格式,GeoJSON幾乎被所有線上地圖框架作為資料來源格式,在geopandas中讀取GeoJSON非常簡單,只需要傳入檔案路徑名稱即可,下面我們來讀入圖13所示的檔案:

圖13
圖14

2.1.4 過濾

  geopandas在0.1.0版本中新增了bbox過濾,在0.7.0版本中新增了蒙版過濾和行過濾功能,可以輔助我們根據自己的需要讀入原始資料中的子集,下面一一進行介紹:

  • bbox過濾

  bbox過濾允許我們在read_file()中傳入一個邊界框作為引數bbox,格式為(左下角x, 左下角y, 右上角x, 右上角y),這樣在讀入的過程中只會保留幾何物件與bbox有相交的資料記錄,下面我們仍然以上文中使用過的中國地圖資料為例,我們在讀入的過程中,傳入邊界框:

from shapely import geometry

data = gpd.read_file('geometry/china_provinces.json',
                    bbox=(100, 20, 110, 30))

%matplotlib widget
ax = data.plot()
# 繪製bbox框示意
ax = gpd.GeoSeries([geometry.box(minx=100, 
                                 miny=20, 
                                 maxx=110, 
                                 maxy=30).boundary]).plot(ax=ax, color='red')
圖15

  可以看到只有跟紅色框有相交的幾何物件被讀入。

  • 蒙版過濾

  蒙版過濾和bbox過濾功能相似,都是篩選與指定區域相交的資料記錄,不同的是蒙版過濾通過mask引數可以傳入任意形狀的多邊形,不再像bbox過濾那樣只接受矩形:

data = gpd.read_file('geometry/china_provinces.json',
                    mask=geometry.Polygon([(100, 20), (110, 30), (120, 20)]))

ax = data.plot()
# 繪製bbox框示意
ax = gpd.GeoSeries([geometry.Polygon([(100, 20), 
                                      (110, 30), 
                                      (120, 20)]).boundary]).plot(ax=ax, color='red')
圖16

  可以看到只有跟紅色多邊形相交的幾何物件被讀入。

  • 行過濾

  行過濾的功能就比較簡單,通過引數rows控制讀入原資料的前若干行,可以用於在讀取大型資料時先快速檢視前幾行以瞭解整個資料的格式:

圖17

2.2 向量檔案的寫出

  在geopandas中使用to_file()來將GeoDataFrameGeoSeries寫出為向量檔案,主要支援shapefileGeoJSON以及GeoPackage,不像geopandas.read_file()可以根據傳入的檔名稱資訊自動推斷型別,我們在寫出向量資料時就需要使用driver引數來宣告檔案型別:

  • ESRI Shapefile

  我們將上文最後一次讀入的GeoDataFrame寫出為ESRI Shapefile,設定driver引數為ESRI Shapefile,如果你對檔案編碼有要求,這裡可以使用encoding引數來指定,譬如這裡我們指定為utf-8

'''在工程根目錄下建立output資料夾'''
import os

try:
    os.mkdir('output')
except FileExistsError:
    pass
    
data.to_file('output/output.shp', 
             driver='ESRI Shapefile',
             encoding='utf-8')

  可以看到在output資料夾下,成功匯出了完整的shapefile

圖18

  而如果匯出的檔名不加字尾副檔名,則會生成包含在新目錄下的shapefile

data.to_file('output/output_shapefile', 
             driver='ESRI Shapefile',
             encoding='utf-8')
圖19

  也可以向指定的資料夾下追加圖層:

data.to_file('output/output_shapefile_multi_layer', 
             driver='ESRI Shapefile',
             layer='layer1',
             encoding='utf-8')

data.to_file('output/output_shapefile_multi_layer', 
             driver='ESRI Shapefile',
             layer='layer2',
             encoding='utf-8')

data.to_file('output/output_shapefile_multi_layer', 
             driver='ESRI Shapefile',
             layer='layer3',
             encoding='utf-8')
圖20
  • GeoPackage

  對於gdb檔案,由於ESRI的限制,暫時無法在開源的geopandas中匯出,但我們可以用QGIS中的GeoPackage作為替代方案(開源世界萬歲O(∩_∩)O~~),只需要將driver引數設定為GPKG即可,這裡需要注意一個bug:在使用geopandas匯出GeoPackage檔案時,可能會出現圖21所示錯誤:

圖21

  但我觀察到即使出現了上述錯誤,GeoPackage檔案也是成功儲存到路徑下的且整個程式並未被打斷,因此可以無視上述錯誤:

圖22
  • GeoJSON

  寫出為GeoJSON非常容易,只需要設定driver='GeoJSON'即可:

圖23

  以上就是本文的全部內容,如有筆誤望指出