python |
您所在的位置:网站首页 › python处理mdb文件 › python |
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 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()#简单展示显示效果: ![]() 要素类的创建效率很高,既能创建要素实体,也能写入属性信息和定义投影。下面先简单介绍下三种要素类的创建方法。 点状要素类的创建![]() ![]() ![]() 高程点文件存储格式与CASS中读取的DAT格式一致,示例:【1,ZDH ,450000.000,4100000,20002,DYG,450000.000,4100000,2000 】其中,“1”代表的是“点号”,“ZDH”代表的是“代码”,之后的分别是“东坐标、北坐标、高程值”即“Y、X、H ”或者是“X、Y、H ” AutoCAD中展点效果![]() ![]() 实现代码 # -*- 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/1QCack83clQEJGorQLaNkKg 提取码:bfrs 最后给大家强烈安利一个geopandas学习博客: https://www.cnblogs.com/feffery/tag/geopandas/ |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |