python geopandas读取、创建shapefile文件的方法
作者:fungis 发布时间:2022-09-23 16:57:19
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版本新增了一些查询、入库方面的特性。
shapefile文件信息的读取
相比pyshp
库,geopandas
库的数据读取、展示、分析、拓展的效果要更好。它可以读取zip中的shapefile
,还可以读取GeoJson、ArcGIS中地理数据库gdb
,以及QGIS
中GeoPackage
存放的矢量数据。
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()#简单展示
显示效果:
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')
线状要素类的创建
核心代码:
# 这里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 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()
点集转面
将一系列点的集合转为面状要素类,下面以甘肃省的地震带为例(字段对应:名称,面索引,点索引,经度,纬度)。
数据预览
效果预览
实现代码
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()
创建缓冲区、多环缓冲区
实现代码:
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/1g7G8sQ17-9XIhojyQ1M7Ww
提取码: 59vz
最后给大家强烈安利一个geopandas学习博客: https://www.cnblogs.com/feffery/tag/geopandas/
来源:https://www.cnblogs.com/fungis/p/14946013.html
猜你喜欢
- 问题:变量名是否合法: 1.变量名可以由字母
- 本文实例讲述了Python pandas RFM模型应用。分享给大家供大家参考,具体如下:什么是RFM模型根据美国数据库营销研究所Arthu
- 大家好,今天分享一个实用的办公脚本:将多个PDF合并为一个PDF,例如我手上现在有如下3个PDF分册,需要整合成一个完整的PDF如果换成你操
- 本文实例讲述了Python设计模式之工厂方法模式。分享给大家供大家参考,具体如下:工厂方法模式(Factory Method Pattern
- 关于python的ssh库操作需要引入一个远程控制的模块——paramiko,可用于对远程服务器进行
- 这篇文章主要介绍了python框架django项目部署相关知识详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价
- 我们在切换选项卡的时候,如果使用的是ajax技术,会碰到如下情况:点击tab1选项,服务器发出一个Ajax请求获取该选项tab1的内容数据。
- 前言最近在功能性测试的过程中,需要在Python环境下用OpenCV读取网络摄像头的视频流,接着用目标检测器进行视屏帧的后续处理。在测试过程
- 当服务器必须提供与两个或更多个网络或网络子网的连接时,典型的方案是使用多宿主计算机。此计算机通常位于外围网络(也称为 DMZ、外围安全区域或
- 使用指令:pyinstaller -F -w main.py生成的.exe文件执行报错看报错信息,初步怀疑是生成的.exe不能执行.ini配
- 插入视图的条件: 1.如果视图是基于一个基础表产生的,那么这就称为非连接视图,所有的非连接视图都是可以更新的,也就是说可以在该视图上进行,I
- Python复合语句复合语句是包含其它语句(语句组)的语句;它们会以某种方式影响或控制所包含其它语句的执行。通常,复合语句会跨越多行,虽然在
- 简介今天试着用ptyhon做了一个抓取网页内容,并生成word文档的功能,功能很简单,做一下记录以备以后用到。生成word用到了第三方组件p
- 在本教程中,我们将学习如何创建一个使用Django作为后端的天气应用程序。Django提供了一个基于Python Web框架的Web框架,允
- 这篇论坛文章详细介绍了完全卸载MySQL数据库5.0的具体方法,更多内容请参考下文:数据库突然出了问题,没办法只能重装,因为事先并不知道My
- 如何在Typescript中使用for...in ?本人在TS中用for...in出现了些问题,也想到了一些解决方法。那么先来看看下面报错的
- 本文实例为大家分享了Python3实现飞机大战游戏的具体代码,供大家参考,具体内容如下1、主程序:plane_main.pyimport p
- 1.普通的输出:print(str)#str是任意一个字符串,数字···2.格式化输出: print('1,2,%s,%d'
- 可扩展标记语言 (XML) 是用于描述数据集内容以及应如何将数据输出到设备上或如何在 Web 页上显示数据的语言。标记语言的创建来源于出版商
- 本文实例讲述了Python实现读取txt文件中的数据并绘制出图形操作。分享给大家供大家参考,具体如下:下面的是某一文本文件中的数据。6.11