shapefile是GIS中非常重要的一種數據類型，在ArcGIS中被稱為要素類(Feature Class)，主要包括點(point)、線(polyline)和多邊形(polygon)。作為一種十分常見的矢量文件格式，geopandas對shapefile提供了很好的讀取和寫出支持，其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版本新增了一些查詢、入庫方面的特性。

相比pyshp庫，geopandas庫的數據讀取、展示、分析、拓展的效果要更好。它可以讀取zip中的shapefile，還可以讀取GeoJson、ArcGIS中地理數據庫gdb，以及QGIS中GeoPackage 存放的矢量數據。

import geopandas as gpdfrom matplotlib import pyplot as pltdata = 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()#簡單展示

顯示效果：

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

點狀要素類的創建

核心代碼：

# 對應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’)

線狀要素類的創建

核心代碼：

# 這里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’)

面狀要素類的創建

核心代碼：

# 對應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中展點效果

geopandas中展點效果

實現代碼

# -*- coding: utf-8 -*-import pandas as pdimport geopandas as gpdfrom shapely.geometry import Pointfrom matplotlib import pyplot as pltfrom 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()

點集轉面

將一系列點的集合轉為面狀要素類，下面以甘肅省的地震帶為例（字段對應：名稱，面索引，點索引，經度，緯度）。

數據預覽

效果預覽

實現代碼

import geopandas as gpdimport pandas as pdfrom shapely.geometry import Polygonfrom matplotlib import pyplot as pltraw = 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’)# 轉換為GeoDataFrameoutput = 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’])]# 導出shapefileoutput.to_file(’output/地震帶.shp’) plt.show()

創建緩沖區、多環緩沖區

實現代碼：

import osimport shapelyimport geopandas as gpdimport matplotlib.pyplot as pltpolygon = 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/1g7G8sQ17-9XIhojyQ1M7Ww

提取碼: 59vz

最后給大家強烈安利一個geopandas學習博客： https://www.cnblogs.com/feffery/tag/geopandas/

以上就是python geopandas讀取、創建shapefile文件的方法的詳細內容，更多關于python讀取shapefile文件的資料請關注好吧啦網其它相關文章！