python-geopandas读取、创建shapefile⽂件、geopandas学习教程python-geopandas读取、创建shapefile⽂件、geopandas学习教程
shapefile是GIS中⾮常重要的⼀种数据类型，在ArcGIS中被称为要素类(Feature Class)，主要包括点(point)、线(polyline)和多边形(polygon)这三种。作为⼀种⼗分常见的⽮量⽂件格式，geopandas对shapefile提供了很好的读取和写出⽀持 。
geopandas库允许对⼏何类型进⾏空间操作，其DataFrame结构相当于GIS数据中的⼀张属性表，使得可以直接操作⽮量数据属性表，使得在python中操作地理数据更⽅便。本⽂给⼤家介绍下⽤Python脚本中对Shapefile⽂件(.shp，.shx，.dbf等格式)进⾏读写操作。
开发准备
由于geopandas有好⼏个依赖库，推荐⼤家使⽤ Miniconda或是 Anaconda来安装geopandas。
安装命令： conda install -c conda-forge geopandaspython怎么读取dat文件
使⽤导⼊：import geopandas
我这⾥⽤的是geopandas 0.7的版本，版本间差异是不太⼤，最新0.8版本新增了⼀些查询、⼊库⽅⾯的特性。
shapefile⽂件信息的读取
相⽐pyshp库，geopandas库的数据读取、展⽰、分析、拓展的效果要更好。它可以读取zip中的shapefile，还可以读取GeoJson、ArcGIS 中地理数据库gdb，以及QGIS中GeoPackage 存放的⽮量数据。
import geopandas as gpd
from matplotlib import pyplot as plt
data = ad_file(r'E:\gisData\⾏政区划数据2019\省.shp')#读取磁盘上的⽮量⽂件
#data = ad_file('shapefile/china.gdb', layer='province')#读取gdb中的⽮量数据
s)# 查看数据对应的投影信息
print(data.head())# 查看前5⾏数据
data.plot()
plt.show()#简单展⽰
显⽰效果：
shapefile⽂件的创建
要素类的创建效率很⾼，既能创建要素实体，也能写⼊属性信息和定义投影。下⾯先简单介绍下三种要素类的创建⽅法。
点状要素类的创建
核⼼代码：
# 对应ry中的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⽂件
<_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py','')),
driver='ESRI Shapefile',
encoding='utf-8')
线状要素类的创建
核⼼代码：
# 这⾥ry.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'])
<_file('./output/{}.shp'.format(os.path.basename(__file__).replace('.py','')),
driver='ESRI Shapefile',
encoding='utf-8')
⾯状要素类的创建
核⼼代码：
# 对应ry中的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
)
<_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 pd
import geopandas as gpd
ry 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⽂件
<_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 =[(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)
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 gpd
import pandas as pd
ry import Polygon
from matplotlib import pyplot as plt
raw = pd.read_excel('./data-use/⽢肃省地震带.xls')# 原始数据
# 转换为⾯要素
output = upby('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 = upby('name', as_index=False,)['longitude'].mean()
new_latitude = upby('name', as_index=False)['latitude'].mean()
new_df = pd.merge(pd.DataFrame(new_longitude),pd.DataFrame(new_latitude))
new_texts =[(x_ , y_ , text, fontsize=8)for x_, y_, text in
zip(new_df['longitude'], new_df['latitude'], new_df['name'])]
# 导出shapefile
<_file('output/地震带.shp')
plt.show()
创建缓冲区、多环缓冲区
实现代码：
import os
import shapely
import geopandas as gpd
import matplotlib.pyplot as plt
polygon = ry.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⽂件
<_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()
写在最后
附相关完整代码的下载，还有更多有趣的内容，感兴趣的朋友们可以⾃⾏实践。喜欢的朋友们可以点个关注，后续将持续更新，精彩⽆限^ -^
提取码：bfrs

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。