批量裁剪矢量地理数据的几种途径
近些天,有一个项目需要对大量的地理数据进行统计,然后在前端展示,作为组里为数不多拥有GIS专业背景的人,毫无意外地任务又落在我身上。
该任务统计过程其实简单,难点在于数据量实在是太大,还要经过空间运算、采样后才能着手统计,而且需求端(前端)的目的一开始不是很明确,扔过来一个截图文件,先是统计这个、后又统计那个、再这个,统计后交换数据结构也是时常变化,本来数据量就大,前后折腾耗了很多时间,又一个面向口头约定编程的经典案例。
如果是小量数据,直接无脑上GIS软件就行,可是它不行,所以我想着写份代码一劳永逸地解决诸如此类的问题,此时眉头一皱,GDAL 涌上心头。
我打算用GDAL 去解决,发现当数据格式为gdb时,它会将整个数据全部读取,但是数据量太大导致异常卡顿,每次读取都要十几分钟,无奈只能将gdb里的每一个图层提取出来,形成一个个的shapefile文件,GDAL 具备空间数据格式转换的能力,同时它还提供了一个超级方便的命令行工具,ogr2ogr(OGR包含自GDAL ,一个用来处理矢量格式数据的库,而GDAL 主要是处理栅格)。
初识ogr2ogr
从pypi上下载的GDAL Python包,安装后是没有[ogr2ogr](https://www.osgeo.cn/GDAL /prOGRams/ogr2ogr.html)的,需前往[Download — GDAL documentation](https://GDAL .org/download.html#current-release),根据系统下载预编译好的版本,在系统高级设置->环境变量中设置好相关OGR变量后,即可在命令行工具中全局调用。
格式转换:gdb ->shapefile
命令行格式
ogr2ogr -f out_put_format output_directory path_to_gdb
output_format:输出格式
output_directory:输出保存路径
path_to_gdb:指定的gdb文件保存路径
示例
ogr2ogr -f "ESRI Shapefile" D:\extracted D:\Geodatabase.gdb
执行该命令后,D:\Geodatabase.gdb
中的所有图层都会被提取并以 Shapefile 格式保存在指定的输出路径中。
ogr2ogr裁剪
使用
下面是ogr2ogr裁剪命令格式
ogr2ogr -clipsrc clip_file output_file input_file
- clip_file:用于裁剪的图层所在路径
- output_file:结果输出路径
- input_file:被裁剪图层所在路径
示例
ogr2ogr -clipsrc E:\Desktop\assets\waters.shp E:\Desktop\clip_test\out.shp E:\Desktop\clip_test\M5_120_0_0000__2D_Zones.shp
注:上述三者文件类型均为shapefile
潜在问题
如果使用上述命令行过程中出现ERROR 6: GEOS support not enabled错误提示,并且在裁剪输出图层内无任何要素,
ERROR 6: GEOS support not enabled
原因是OGR并没有对特征之间的空间关系计算提供原生支持,该功能依赖另一个C++库——GEOS。GEOS 是一个用于地理空间计算的库,它在 GDAL 中负责处理几何计算任务,包括裁剪。如果在编译 GDAL 时未启用 GEOS 支持,将无法使用裁剪功能。
要解决此问题,可以尝试重新编译和安装 GDAL ,并确保在配置和编译过程中启用 GEOS 支持,或者前往GISInternals Support Site下载预编译好,且支持GEOS的GDAL 版本。
批量裁剪
ogr2ogr
待处理文件数量很多,如何批量裁剪,目前我能使出三种方法,第一种还是使用ogr2ogr,但它作为命令行工具,如何才能批量获取shapefile文件,然后批量调用它执行裁剪呢,既然是命令行,肯定是写bat脚了,但是那玩意我不会,所以借助Python,代码如下:
import subprocess import os def batch_clip(): out_root = r'output_file_root_path' clip_file = r'path_to_clip_shapefile' for root, folders, files in os.walk(r'all_shapefile_path'): for f in files: if f.endswith('.shp'): input_file = root + '\\' + f output_file = out_root + '\\' + f print("{:-^30}".format(f)) # 构造 ogr2ogr 命令 command = ['ogr2ogr', '-clipsrc', clip_file, output_file, input_file] # 执行命令 subprocess.call(command) print('============完成==============')
ArcPy
ArcGIS中,几乎每一个地理人估计都有的工具,它执行空间裁剪相当简单,其批处理程序具备解决重复性操作的能力,但是数量上升到成百上千时,效率不高,故借用它的Python开发包进行定制批处理操作。
注意:ArcGIS内置Python2,版本太老,语法和3不兼容,诸如文件路径出现中文等小细节会导致各种各样的问题,弃用,选择它的理由也就是地理人人手一个,唾手可得,不用再装其他依赖。
# -*- coding: cp936 -*- import os import arcpy clip = r"path_to_clip_shapefile" for root, folders, files in os.walk(r'shapefile_path'): for f in files: if f.endswith('.shp'): print('----%s' % f) input_file = root + '\\' + f output_file = root + '\\' + 'clipped' + '\\' + f # 读取裁剪范围的几何对象 with arcpy.da.SearchCursor(clip, ["SHAPE@"]) as cursor: for row in cursor: clip_geometry = row[0] # 执行裁剪 arcpy.Clip_analysis(input_file, clip_geometry, output_file)
定义输入输出变量如下,
input_file
:输出图层的文件路径
output_file
: 和裁剪范围的文件路径
clip_feature
:裁剪图层。
然后,使用 arcpy.da.SearchCursor()
从裁剪范围图层中获取几何对象。在此示例中,假设裁剪范围图层只有一个要素,因此使用了 for
循环来遍历记录并获取 SHAPE@
字段的值(即几何对象)。
最后,使用 arcpy.Clip_analysis()
函数执行裁剪操作,将输入图层按照裁剪范围进行裁剪,并保存到输出图层中。
请注意替换示例代码中的文件路径为实际的文件路径,并确保已经设置好 ArcGIS 环境以及相关的许可证。
Geopandas
Geopandas是一个新兴的,基于Python的地理空间处理包,(看名字就知道它和pandas有瓜葛)同样支持空间裁剪操作。
import geopandas as gpd import os for root, folders, files in os.walk(r'shapefile_path'): for f in files: if f.endswith('.shp'): print('----%s' % f) # 读取输入图层和裁剪图层 input_layer = gpd.read_file(root + '\\' + f) clip_layer = gpd.read_file(r'path_to_clip_layer') output_name =root + r"\cliped_" + f # 进行裁剪操作 output_layer = gpd.clip(input_layer, clip_layer) # 将裁剪结果保存为新图层 output_layer.to_file(output_name) print('============完成==============')
这样,空间运算的部分结束,之后再统一遍历所有的shapefile,统计目标信息就行了。
参考
- GDAL (OGR)库与GEOS库的协作 开源地理空间基金会中文分会 开放地理空间实验室 (osgeo.cn)
- 小助手(一个基于chatGPT的微信机器人)
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 实操Deepseek接入个人知识库
· 易语言 —— 开山篇
· Trae初体验