python-geopandas讀取、建立shapefile檔案

fungis發表於2021-06-28

作者:fungis 描述:一個熱帶生活、樂於分享、努力搬磚的giser 交流郵箱:fungis@163.com

shapefile是GIS中非常重要的一種資料型別,在ArcGIS中被稱為要素類(Feature Class),主要包括點(point)、線(polyline)和多邊形(polygon)。作為一種十分常見的向量檔案格式,geopandasshapefile提供了很好的讀取和寫出支援,其DataFrame結構相當於GIS資料中的一張屬性表,使得可以直接操作向量資料屬性表,使得在python中操作地理資料更方便。本文給大家介紹下用Python指令碼中對Shapefile檔案(.shp,.shx,.dbf等格式)進行讀寫操作。

開發準備

由於geopandas有好幾個依賴庫,推薦大家使用 Miniconda或是 Anaconda來安裝geopandas。

安裝命令:

conda install -c conda-forge geopandas   

國內映象:

conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge geopandas                   

使用匯入:import geopandas

我這裡用的是geopandas 0.7的版本,版本間差異是不太大,最新0.8版本新增了一些查詢、入庫方面的特性。

shapefile檔案資訊的讀取

相比pyshp庫,geopandas庫的資料讀取、展示、分析、擴充的效果要更好。它可以讀取zip中的shapefile,還可以讀取GeoJson、ArcGIS中地理資料庫gdb,以及QGISGeoPackage 存放的向量資料。

import geopandas as gpd
from matplotlib import pyplot as plt

data = gpd.read_file(r'E:\gisData\行政區劃資料2019\省.shp')#讀取磁碟上的向量檔案
#data = gpd.read_file('shapefile/china.gdb', layer='province')#讀取gdb中的向量資料
print(data.crs)  # 檢視資料對應的投影資訊
print(data.head())  # 檢視前5行資料
data.plot()
plt.show()#簡單展示

顯示效果:

python-geopandas讀取、建立shapefile檔案

shapefile檔案的建立

要素類的建立效率很高,既能建立要素實體,也能寫入屬性資訊和定義投影。下面先簡單介紹下三種要素類的建立方法。

點狀要素類的建立

python-geopandas讀取、建立shapefile檔案
  • 核心程式碼:
# 對應shapely.geometry中的Point,用於表示單個點,下面我們建立一個由若干Point物件組成
cq = geopandas.GeoSeries([geometry.Point(110, 60),
                          geometry.Point(110.5, 50.4),
                          geometry.Point(120, 55),
                          geometry.Point(107.8, 54.6),
                          geometry.Point(114.6, 50)],
                         crs='EPSG:4326',  # 指定座標系為WGS 1984
                         index=['一號', '二號', '三號', '四號', '五號'],  # 相關的索引
                         )
# 匯出資料為shapefile檔案
cq.to_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py', '')),
           driver='ESRI Shapefile',
           encoding='utf-8')

線狀要素類的建立

python-geopandas讀取、建立shapefile檔案
  • 核心程式碼:
# 這裡shapely.geometry.LineString([(x1, y1), (x2, y2), ...])用於建立多點按順序連線而成的線段
cq = geopandas.GeoSeries([geometry.LineString([(0, 0), (1, 1), (1, 0)]),
                          geometry.LineString([(0.5, 2), (0, 1), (-1, 0)])],
                         crs='EPSG:4326',
                         index=['一號線', 'b'])
cq.to_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py', '')),
           driver='ESRI Shapefile',
           encoding='utf-8')

面狀要素類的建立

python-geopandas讀取、建立shapefile檔案
  • 核心程式碼:
# 對應shapely.geometry中的Polygon,用於表示面,下面我們建立一個由若干Polygon物件組成
cq = geopandas.GeoSeries([geometry.Polygon([(14, 14), (13, 18), (20, 11), (18, 10)]),
                          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))]),
                          geometry.Polygon([(11, 2), (11, 10), (12, 10), (12, 2)])
                          ],
                         index=['簡單面', '複雜面', 'c區'],  # 構建一個索引欄位
                         crs='EPSG:4326',  # 座標系是:WGS 1984
                         )
cq.to_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py', '')),
           driver='ESRI Shapefile',
           encoding='utf-8')

擴充應用例項

展高程點

高程點檔案儲存格式與CASS中讀取的DAT格式一致,示例:【1,ZDH ,450000.000,4100000,20002,DYG,450000.000,4100000,2000 】其中,“1”代表的是“點號”,“ZDH”代表的是“程式碼”,之後的分別是“東座標、北座標、高程值”即“Y、X、H ”或者是“X、Y、H ”

  • AutoCAD中展點效果
python-geopandas讀取、建立shapefile檔案
  • geopandas中展點效果
python-geopandas讀取、建立shapefile檔案
  • 實現程式碼

    # -*- coding: utf-8 -*-
    
    import pandas as pd
    import geopandas as gpd
    from shapely.geometry import Point
    from matplotlib import pyplot as plt
    from matplotlib.ticker import FuncFormatter
    
    # 讀取資料
    file_path = './data-use/高程資料.csv'
    rankings_colname = ['name', 'mark', 'longitude', 'latitude', 'height'];
    df = pd.read_csv(file_path, header=None, names=rankings_colname)
    # print(df.head(5))#輸出前五行資料檢視
    xy = [Point(xy) for xy in zip(df['longitude'], df['latitude'])]
    pts = gpd.GeoSeries(xy)  # 建立點要素資料集
    #儲存為SHP檔案
    pts.to_file('./output/展高程點.shp', driver='ESRI Shapefile', encoding='utf-8')
    """fig是用來設定影像大小引數,ax是行列有多少個點"""
    fig, ax = plt.subplots(figsize=(8, 6))  # 返回一個包含figure和axes物件的元組
    ax = pts.plot(ax=ax,
                  facecolor='white',
                  edgecolor='black',
                  marker='X',
                  linewidth=0.5,  # 內外符號比例係數
                  markersize=12,
                  label='高程點')
    # 地圖示註
    new_texts = [plt.text(x_ + 1, y_ + 1, text, fontsize=8) for x_, y_, text in
                 zip(df['longitude'], df['latitude'], df['name'])]
    
    
    # 設定座標軸
    def formatnum(x, pos):
        # return '$%.1f$x$10^{4}$' % (x / 10000)#科學計數法顯示
        return int(x)  # 取整顯示
    
    
    formatter = FuncFormatter(formatnum)
    ax.yaxis.set_major_formatter(formatter)
    
    # 美觀起見隱藏頂部與右側邊框線
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    plt.grid(True, alpha=0.4)  # 顯示網格,透明度為50%
    ax.legend(title="圖例", loc='lower right', ncol=1, shadow=True) # 新增圖例
    plt.title('展高程點', fontdict={'weight': 'normal', 'size': 20}) # 設定圖名&改變圖示題字型
    # 儲存圖片
    plt.savefig('images/展高程點.png', dpi=300, bbox_inches='tight', pad_inches=0)
    plt.show()
    

點集轉面

將一系列點的集合轉為面狀要素類,下面以甘肅省的地震帶為例(欄位對應:名稱,面索引,點索引,經度,緯度)。

  • 資料預覽
python-geopandas讀取、建立shapefile檔案
  • 效果預覽
python-geopandas讀取、建立shapefile檔案
python-geopandas讀取、建立shapefile檔案
  • 實現程式碼

    import geopandas as gpd
    import pandas as pd
    from shapely.geometry import Polygon
    from matplotlib import pyplot as plt
    
    raw = pd.read_excel('./data-use/甘肅省地震帶.xls')  # 原始資料
    # 轉換為面要素
    output = raw.groupby('id') \
        .apply(lambda df: Polygon([(x, y) for x, y in zip(df['longitude'], df['latitude'])])) \
        .to_frame(name='geometry')
    
    # 轉換為GeoDataFrame
    output = gpd.GeoDataFrame(output, crs='EPSG:4326')
    output.plot()
    # 地圖示註
    new_longitude = raw.groupby('name', as_index=False,)['longitude'].mean()
    new_latitude = raw.groupby('name', as_index=False)['latitude'].mean()
    new_df = pd.merge(pd.DataFrame(new_longitude),pd.DataFrame(new_latitude))
    new_texts = [plt.text(x_ , y_ , text, fontsize=8) for x_, y_, text in
                 zip(new_df['longitude'], new_df['latitude'], new_df['name'])]
    # 匯出shapefile
    output.to_file('output/地震帶.shp')  
    plt.show()
    

建立緩衝區、多環緩衝區

python-geopandas讀取、建立shapefile檔案
  • 實現程式碼:

    import os
    import shapely
    import geopandas as gpd
    import matplotlib.pyplot as plt
    
    polygon = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
    # 分別繪製多邊形、多邊形正向緩衝區,座標系是WGS1984,單位是度
    cq = gpd.GeoSeries([polygon,
                        polygon.buffer(distance=1),
                        polygon.buffer(distance=3)],
                       crs='EPSG:4326')
    # 匯出資料為shapefile檔案
    cq.to_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py', '')),
               driver='ESRI Shapefile',
               encoding='utf-8')
    ax = cq.plot(alpha=0.2)
    ax.axis('off')  # 取消座標軸的顯示
    plt.show()
    

寫在最後

附相關完整程式碼的下載,還有更多有趣的內容,感興趣的朋友們可以自行實踐。喜歡的朋友們可以點個關注,後續將持續更新,精彩無限^ - ^

連結:https://pan.baidu.com/s/1QCack83clQEJGorQLaNkKg

提取碼:bfrs

最後給大家強烈安利一個geopandas學習部落格: https://www.cnblogs.com/feffery/tag/geopandas/

相關文章