目录

• • • 初识ogr2ogr
      • 格式转换：gdb ->shapefile
    • ogr2ogr裁剪
      • 使用
      • 潜在问题
    • 批量裁剪
      • ogr2ogr
      • ArcPy
      • Geopandas
  • 参考

[TOC]

近些天，有一个项目需要对大量的地理数据进行统计，然后在前端展示，作为组里为数不多拥有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，统计目标信息就行了。

参考

1. GDAL (OGR)库与GEOS库的协作 开源地理空间基金会中文分会 开放地理空间实验室 (osgeo.cn)
2. 小助手（一个基于chatGPT的微信机器人)