如何使用OGR提取OpenStreetMap Shapefile中指定要素（如自行车道）

嗨，这事儿好办！咱们可以用OGR的SetAttributeFilter()方法直接按fclass字段筛选要素，最后把筛选结果导出成新的Shapefile，完美实现减小文件体积的需求。下面是完整的代码示例和步骤说明：

完整代码实现 #

from osgeo import ogr

# 注册所有OGR驱动（确保能正常读写Shapefile）
ogr.RegisterAll()

# 打开原始道路Shapefile
input_path = 'gis.osm_roads_free_1.shp'
input_ds = ogr.Open(input_path, 0)  # 0代表只读模式
if not input_ds:
    raise Exception(f"无法打开目标文件：{input_path}")

# 获取原始图层
input_layer = input_ds.GetLayer(0)

# 设置属性过滤器：筛选fclass为自行车道的要素
# 常见自行车道的fclass值是'cycleway'，如果需要包含多类可以用IN语法，比如"fclass IN ('cycleway', 'path')"
filter_query = "fclass = 'cycleway'"
input_layer.SetAttributeFilter(filter_query)

# 创建输出Shapefile（保存过滤后的结果）
output_path = 'filtered_cycleways.shp'
# 获取Shapefile驱动
driver = ogr.GetDriverByName('ESRI Shapefile')
# 如果输出文件已存在，先删除旧文件
if driver.DeleteDataSource(output_path) == 0:
    print(f"已移除原有输出文件：{output_path}")

# 创建新的数据源
output_ds = driver.CreateDataSource(output_path)
if not output_ds:
    raise Exception(f"无法创建输出文件：{output_path}")

# 复制原始图层的空间参考和字段结构
srs = input_layer.GetSpatialRef()
output_layer = output_ds.CreateLayer('cycleways', srs=srs, geom_type=input_layer.GetGeomType())

# 把原始图层的所有字段定义复制到新图层
layer_def = input_layer.GetLayerDefn()
for i in range(layer_def.GetFieldCount()):
    field_def = layer_def.GetFieldDefn(i)
    output_layer.CreateField(field_def)

# 遍历筛选后的要素，逐个复制到新图层
feature = input_layer.GetNextFeature()
while feature:
    output_layer.CreateFeature(feature)
    feature = input_layer.GetNextFeature()

# 释放资源
input_ds = None
output_ds = None

print(f"筛选完成！已将自行车道要素保存到：{output_path}")

关键步骤说明 #

• 注册驱动：ogr.RegisterAll()是个稳妥的操作，确保OGR能识别Shapefile格式，避免出现读写失败的情况。
• 属性过滤：SetAttributeFilter()接收类似SQL WHERE子句的字符串，能精准匹配fclass字段的值。你可以根据实际需求调整过滤条件，比如要包含其他自行车相关道路类型，就用IN语法扩展范围。
• 导出新文件：直接导出过滤后的要素是减小文件体积的核心，新文件只会保留你需要的内容，比原文件精简很多。
• 资源清理：最后把数据源设为None，能及时释放文件句柄，避免资源占用。

如果想确认筛选结果的数量，可以在设置过滤后加一行代码：print(f"共筛选到{input_layer.GetFeatureCount()}个自行车道要素")，方便快速验证结果是否符合预期。

内容的提问来源于stack exchange，提问作者Yorian