GDAL Python操作shp文件：读取、新建与更新实战指南

1. 从一次数据交付事故说起：为什么shp操作值得单独拎出来讲

去年帮一个做城市绿化的团队处理一批地块数据，对方发来一个压缩包，里面是十几个shp文件。我当时的想法很简单——用Python批量读一遍，做个面积汇总，再按区域字段更新一下属性表，半小时的事。结果打开第一个文件就卡住了：属性表里中文字段全是乱码，坐标系标注和实际数据对不上，还有一个文件的几何类型混着面和线。那天下午我花了整整四个小时才把数据理干净，真正写代码的时间不到二十分钟。

这件事让我意识到，shp这种格式看着简单，坑却全藏在细节里。它不像GeoJSON那样一个文件搞定所有事，shp是一套文件组合——.shp存几何、.shx存索引、.dbf存属性、.prj存坐标系，少一个都可能出问题。而GDAL作为地理数据处理领域的老牌底层库，Python绑定（也就是osgeo.ogr和osgeo.gdal）提供了对shp最完整的控制能力，但它的API设计带着明显的C++血统，新手直接上手容易懵。

这篇内容就是把我这些年用GDAL Python操作shp的经验整理出来，重点讲三件事：怎么正确读取shp（包括那些容易翻车的编码和坐标系问题）、怎么从零新建一个shp、怎么安全地更新已有shp的属性。适合已经会一点Python、但在地理数据处理上还没形成系统方法的朋友。如果你正在做GIS相关的数据处理、空间分析或者数据交付，这些内容应该能帮你省下不少排查时间。

提示：本文所有代码基于GDAL 3.x版本，Python 3.8以上环境。GDAL 2.x和3.x在部分API上有差异，如果你用的是旧版本，注意对照官方迁移文档调整。

2. 读取shp之前，先把这几个概念理清楚

2.1 shp文件组的真实构成与GDAL的打开逻辑

很多人第一次用GDAL打开shp时会困惑：为什么代码里写的是.shp路径，但实际读取时GDAL好像把整个文件夹都扫了一遍？这是因为shp本质上是一个数据集（Dataset），而不是单个文件。GDAL在打开.shp时，会自动去同目录下找同名的.shx、.dbf、.prj、.cpg等配套文件。

.cpg这个文件特别值得说。它记录的是.dbf属性表的编码方式，但很多数据生产工具压根不生成它。没有.cpg的时候，GDAL只能靠猜，猜错了就是满屏乱码。我见过太多人拿到乱码数据第一反应是"GDAL有bug"，其实问题出在数据本身缺少编码声明。

打开shp的标准写法：

from osgeo import ogr # 开启异常抛出，否则GDAL出错时只返回None，排查起来很痛苦 ogr.UseExceptions() driver = ogr.GetDriverByName("ESRI Shapefile") dataset = driver.Open("/path/to/your_data.shp", 0) # 0表示只读，1表示可写 if dataset is None: raise FileNotFoundError("shp文件打开失败，检查路径和配套文件是否完整") layer = dataset.GetLayer(0) # shp通常只有一个图层 print(f"要素数量: {layer.GetFeatureCount()}") print(f"几何类型: {layer.GetGeomType()}")

这里有个经验：永远先调ogr.UseExceptions()。GDAL默认的错误处理是静默返回None，你拿到None之后根本不知道是路径错了、文件损坏了还是权限不够。开启异常之后，报错信息会直接告诉你问题在哪。

2.2 编码问题：乱码的根源与三种解法

中文乱码是shp读取里最高频的问题。根源在于.dbf文件的编码没有统一标准，有的用GBK，有的用UTF-8，还有用Latin-1的。GDAL在打开时会尝试用.cpg文件里的声明，没有.cpg就默认按Latin-1处理，中文自然全乱。

解法分三种情况：

第一种，数据自带.cpg且声明正确。这种情况直接读就行，GDAL会自动处理。你可以用下面这行代码确认当前编码：

# 查看GDAL识别到的编码 print(layer.GetMetadataItem("ENCODING"))

第二种，没有.cpg但你知道实际编码。可以通过配置选项强制指定：

from osgeo import gdal # 强制按GBK读取 gdal.SetConfigOption("SHAPE_ENCODING", "GBK") dataset = driver.Open("/path/to/your_data.shp", 0)

注意这个配置是全局的，读完记得重置：gdal.SetConfigOption("SHAPE_ENCODING", None)，否则会影响后续其他文件的读取。

第三种，编码不确定。我的做法是先读一条记录出来看，用feature.GetField()拿到原始字节，然后分别用GBK和UTF-8试解码，哪个能正常显示中文就用哪个。这个方法虽然笨，但在处理来源不明的数据时最可靠。

注意：如果.dbf里同时存在中文和特殊符号，GBK和UTF-8都可能部分失败。这种情况建议先用其他工具（比如QGIS）转成UTF-8再处理，不要在代码里硬扛。

2.3 坐标系识别：别信.prj文件的一面之词

.prj文件存的是WKT格式的坐标系描述，但它的可靠性参差不齐。我遇到过.prj写着WGS84但实际数据是CGCS2000的，也遇到过.prj文件直接是空的。GDAL读取坐标系的方式：

spatial_ref = layer.GetSpatialRef() if spatial_ref is not None: print(spatial_ref.ExportToWkt()) # 获取EPSG代码（如果能识别的话） print(spatial_ref.GetAuthorityCode(None)) else: print("该图层没有坐标系信息")

如果GetAuthorityCode返回None，说明GDAL无法自动匹配到标准EPSG代码。这时候你需要根据数据来源判断实际坐标系，然后用ImportFromEPSG手动设置。坐标系搞错的后果很严重——面积计算偏差可能达到千米级别，空间叠加分析会完全错位。

一个实用的判断技巧：拿几个已知地物的坐标点，和在线地图上的实际位置对比。如果偏差在几十米内，大概率是投影坐标系；如果偏差在经纬度量级（比如经度116左右、纬度39左右），那就是地理坐标系。

3. 逐要素读取：从几何到属性的完整提取

3.1 遍历要素的正确姿势与性能考量

读取shp的核心操作就是遍历要素。最基础的写法：

layer.ResetReading() # 重置游标到开头 for feature in layer: geom = feature.GetGeometryRef() name = feature.GetField("NAME") area = feature.GetField("AREA") print(f"{name}: {area}")

这里有个容易忽略的点：layer对象是有游标的。如果你之前已经遍历过一遍，再遍历时需要先ResetReading()，否则会直接从上次结束的位置继续，拿不到数据。我第一次踩这个坑的时候，以为是数据丢了，排查了半天。

性能方面，shp的读取速度主要受两个因素影响：要素数量和字段数量。一万条以内的数据，直接遍历完全没问题。如果超过十万条，建议用GetNextFeature()手动控制循环，配合批量处理逻辑，避免一次性把所有要素加载到内存。

layer.ResetReading() while True: feature = layer.GetNextFeature() if feature is None: break # 处理逻辑 feature = None # 及时释放，减少内存占用

feature = None这行看起来多余，但在处理大量数据时能明显降低内存峰值。GDAL的Python绑定底层是C++对象，Python的垃圾回收不会立即释放，手动置None是触发释放的有效方式。

3.2 几何对象的读取与常用操作

拿到geom之后，能做很多事情。常用的几个操作：

geom = feature.GetGeometryRef() # 几何类型 geom_type = geom.GetGeometryName() # 返回 "POLYGON", "POINT", "LINESTRING" 等 # 面积和长度（注意：单位取决于坐标系） area = geom.GetArea() length = geom.GetLength() # 获取坐标点 if geom_type == "POINT": x, y = geom.GetX(), geom.GetY() elif geom_type == "POLYGON": # 外环 ring = geom.GetGeometryRef(0) points = [(ring.GetX(i), ring.GetY(i)) for i in range(ring.GetPointCount())] # 边界框 minx, maxx, miny, maxy = geom.GetEnvelope()

面积计算有个大坑：如果数据是地理坐标系（比如WGS84），GetArea()返回的是平方度，这个值没有任何实际意义。要得到平方米，必须先把几何投影到投影坐标系，或者用OGRCoordinateTransformation做动态投影。我见过有人直接拿平方度当平方米用，结果算出来的面积差了十几个数量级。

正确的做法是：

from osgeo import osr # 假设源数据是WGS84，目标用Web Mercator source_srs = layer.GetSpatialRef() target_srs = osr.SpatialReference() target_srs.ImportFromEPSG(3857) transform = osr.CoordinateTransformation(source_srs, target_srs) geom_projected = geom.Clone() geom_projected.Transform(transform) area_sqm = geom_projected.GetArea()

注意Clone()这一步不能省，直接Transform会修改原始几何对象，影响后续操作。

3.3 属性字段的读取与类型处理

GetField()返回的Python类型取决于.dbf里的字段类型。整型返回int，浮点返回float，字符串返回str，日期返回的是字符串（GDAL不会自动转成datetime）。有个细节：如果字段值为空，GetField()返回None，直接参与计算会报TypeError。

安全的读取方式：

def safe_get_field(feature, field_name, default=0): value = feature.GetField(field_name) return value if value is not None else default area = safe_get_field(feature, "AREA", 0.0) name = safe_get_field(feature, "NAME", "")

字段名的大小写也需要注意。shp的.dbf字段名在Windows下不区分大小写，但在Linux下区分。如果你的代码要跨平台运行，建议先用layer.GetLayerDefn()拿到字段定义，确认实际字段名：

layer_defn = layer.GetLayerDefn() field_names = [layer_defn.GetFieldDefn(i).GetName() for i in range(layer_defn.GetFieldCount())] print(field_names)

4. 从零新建shp：字段定义与几何写入的完整流程

4.1 创建前的准备工作与驱动选择

新建shp和读取是两套逻辑。读取是打开已有数据集，新建是从驱动创建。GDAL里创建shp用的是同一个ESRI Shapefile驱动，但调用的是CreateDataSource方法：

from osgeo import ogr ogr.UseExceptions() driver = ogr.GetDriverByName("ESRI Shapefile") # 如果文件已存在，先删除（CreateDataSource不会覆盖） import os output_path = "/path/to/new_data.shp" if os.path.exists(output_path): driver.DeleteDataSource(output_path) dataset = driver.CreateDataSource(output_path)

这里有个关键点：CreateDataSource不会自动覆盖已有文件。如果目标路径已经存在同名shp，它会直接返回None或者报错。所以创建前先检查并删除是标准操作。DeleteDataSource会同时删掉配套的.shx、.dbf等文件，比手动os.remove更干净。

4.2 字段定义的细节：类型、宽度与精度

创建图层时需要指定几何类型和字段定义。字段定义是新手最容易出错的地方，因为.dbf的字段类型和Python类型不是一一对应的：

创建图层的完整代码：

# 创建图层，指定几何类型为多边形 layer = dataset.CreateLayer("new_layer", geom_type=ogr.wkbPolygon) # 添加字段 # 字符串字段必须指定宽度，否则默认宽度可能不够 layer.CreateField(ogr.FieldDefn("NAME", ogr.OFTString).SetWidth(100)) layer.CreateField(ogr.FieldDefn("AREA", ogr.OFTReal)) layer.CreateField(ogr.FieldDefn("COUNT", ogr.OFTInteger)) layer.CreateField(ogr.FieldDefn("UPDATE_DATE", ogr.OFTDate))

字符串宽度是个隐藏的坑。.dbf格式对字符串字段有宽度限制，默认是80个字符。如果你要存更长的文本（比如地址描述），必须显式设置SetWidth。超过254个字符的文本，.dbf存不下，需要考虑用其他格式或者拆分字段。

浮点字段还可以设置精度：SetPrecision(2)表示保留两位小数。但注意这个精度只影响显示，不影响实际存储的精度。

4.3 几何对象的构建与写入

创建几何对象有两种方式：一种是用ogr.Geometry手动构建，另一种是用WKT字符串创建。手动构建更灵活：

# 创建一个多边形 ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(116.0, 39.0) ring.AddPoint(116.1, 39.0) ring.AddPoint(116.1, 39.1) ring.AddPoint(116.0, 39.1) ring.AddPoint(116.0, 39.0) # 闭合环，首尾点必须相同 polygon = ogr.Geometry(ogr.wkbPolygon) polygon.AddGeometry(ring) # 创建要素并设置几何和属性 feature = ogr.Feature(layer.GetLayerDefn()) feature.SetGeometry(polygon) feature.SetField("NAME", "测试地块") feature.SetField("AREA", 12345.67) feature.SetField("COUNT", 10) feature.SetField("UPDATE_DATE", "2024-01-15") layer.CreateFeature(feature) feature = None # 释放

多边形环必须闭合，这是OGC几何规范的要求。首尾点坐标相同才算闭合，否则QGIS等软件打开时会报几何无效。我见过有人用循环生成点的时候忘了补最后一个点，结果数据能用但显示异常，排查起来很费时间。

写入完成后记得关闭数据集：

dataset = None # 触发flush和关闭

不关闭的话，数据可能还在内存缓冲区里，文件不完整。这是GDAL的C++底层机制决定的，Python的with语句对GDAL数据集不生效，必须手动置None。

4.4 批量写入的性能优化

如果要写入大量要素，逐条CreateFeature会比较慢。GDAL提供了事务机制来加速：

layer.StartTransaction() for i in range(10000): feature = ogr.Feature(layer.GetLayerDefn()) # ... 设置几何和属性 layer.CreateFeature(feature) feature = None layer.CommitTransaction()

事务能把一万条要素的写入时间从几十秒降到几秒。原理是减少了磁盘I/O次数，把多次小写入合并成一次大提交。但事务不是越大越好，如果数据量特别大（比如百万级），建议分批提交，每批一万条左右，避免内存占用过高。

5. 更新已有shp：属性修改与几何替换的边界

5.1 以可写模式打开与字段值的修改

更新shp的第一步是用可写模式打开：

dataset = driver.Open("/path/to/data.shp", 1) # 1表示可写 layer = dataset.GetLayer(0)

然后遍历要素，修改字段值，调用SetFeature保存：

layer.ResetReading() for feature in layer: old_value = feature.GetField("STATUS") if old_value == "待处理": feature.SetField("STATUS", "已处理") layer.SetFeature(feature) # 必须调用，否则修改不生效

SetFeature是必须的，很多人改了字段值之后发现文件没变化，就是因为漏了这一步。SetFeature会把修改后的要素写回数据集。

5.2 新增字段与删除字段的操作

给已有shp添加字段：

layer.CreateField(ogr.FieldDefn("NEW_FIELD", ogr.OFTString).SetWidth(50))

注意CreateField在可写模式下才能调用，而且添加字段后，已有要素的新字段值为空（None）。如果需要填充默认值，得再遍历一遍。

删除字段用DeleteField：

layer.DeleteField(layer.FindFieldIndex("OLD_FIELD", True))

FindFieldIndex的第二个参数True表示不区分大小写。删除字段是不可逆操作，建议先备份数据。

注意：shp格式对字段数量有限制，最多255个字段。而且字段名最长10个字符（超过会被截断）。如果你的字段名比较长，创建时就要注意，否则后续按名字查找会失败。

5.3 几何更新的注意事项

更新几何和更新属性类似，但有个额外步骤——修改几何后需要重新设置：

for feature in layer: geom = feature.GetGeometryRef() # 对几何做某种变换，比如平移 geom.Transform(osr.CoordinateTransformation(source_srs, target_srs)) feature.SetGeometry(geom) layer.SetFeature(feature)

这里有个陷阱：GetGeometryRef()返回的是内部几何的引用，直接修改它会影响原始数据。如果你只是想读取几何做判断，不想修改，应该用geom.Clone()复制一份再操作。

另外，shp不支持混合几何类型。一个shp文件里只能有一种几何类型（全是点、全是线或全是面）。如果你试图往面图层里写入点几何，GDAL会报错或者静默丢弃。这是shp格式本身的限制，不是GDAL的问题。

5.4 更新操作中的事务与回滚

批量更新时同样可以用事务：

layer.StartTransaction() try: for feature in layer: # 更新逻辑 layer.SetFeature(feature) layer.CommitTransaction() except Exception as e: layer.RollbackTransaction() print(f"更新失败，已回滚: {e}")

RollbackTransaction能在出错时撤销所有未提交的修改，保证数据一致性。这个机制在处理关键数据时特别有用，避免改了一半出错导致数据处于中间状态。

6. 那些文档里不会写的实战经验

6.1 中文路径与跨平台兼容问题

GDAL在Windows下对中文路径的支持一直不太稳定。我遇到过路径里有中文时Open返回None的情况，换成英文路径就正常。如果你必须在中文路径下工作，有两个方案：一是用os.path做路径编码转换，二是先把文件复制到临时英文目录处理完再复制回去。

跨平台方面，Windows和Linux对路径分隔符的处理不同，建议统一用os.path.join拼接路径，不要手写/或\。

6.2 大文件处理的内存控制

处理几十万条要素的shp时，内存是个绕不开的问题。除了前面提到的GetNextFeature逐条读取，还可以用layer.SetAttributeFilter()做属性过滤，只加载符合条件的要素：

layer.SetAttributeFilter("AREA > 1000")

这个过滤是在GDAL底层执行的，比在Python里遍历判断效率高得多。用完记得layer.SetAttributeFilter(None)清除过滤条件。

空间过滤用layer.SetSpatialFilter(geom)，传入一个几何对象，只返回与之相交的要素。做区域裁剪的时候特别有用。

6.3 常见报错与排查思路

geom.IsValid()是个很实用的方法，能在写入前检查几何是否合法。对于从其他来源拼接的几何，建议写入前都过一遍这个检查。

6.4 与其他工具的配合使用

GDAL Python不是万能的。有些操作，比如复杂的几何修复、拓扑检查，用QGIS的Processing工具或者PostGIS更高效。我的习惯是：GDAL负责批量读写和格式转换，复杂空间分析交给专业工具。两者配合的方式是GDAL输出中间结果，QGIS处理后再用GDAL读回来。

另外，如果项目里已经用了geopandas，简单的shp读写可以直接用gpd.read_file()，底层也是GDAL，但API更友好。只有在需要精细控制（比如指定编码、处理异常几何）时，才需要直接调GDAL。

7. 一个完整的实战案例：地块数据批量更新

假设有一个地块shp，包含字段NAME（地块名）、AREA（面积）、STATUS（状态），现在需要做三件事：把所有STATUS为"待更新"的记录改成"已更新"，给面积大于5000的地块添加一个LEVEL字段标记为"大型"，并输出一个新的shp文件。

from osgeo import ogr import os ogr.UseExceptions() driver = ogr.GetDriverByName("ESRI Shapefile") input_path = "/data/plots.shp" output_path = "/data/plots_updated.shp" # 删除已存在的输出文件 if os.path.exists(output_path): driver.DeleteDataSource(output_path) # 打开源数据 src_dataset = driver.Open(input_path, 0) src_layer = src_dataset.GetLayer(0) # 创建输出数据 out_dataset = driver.CreateDataSource(output_path) out_layer = out_dataset.CreateLayer("plots", geom_type=src_layer.GetGeomType()) # 复制字段定义 src_defn = src_layer.GetLayerDefn() for i in range(src_defn.GetFieldCount()): out_layer.CreateField(src_defn.GetFieldDefn(i)) # 添加新字段 out_layer.CreateField(ogr.FieldDefn("LEVEL", ogr.OFTString).SetWidth(20)) # 遍历并处理 src_layer.ResetReading() out_layer.StartTransaction() for feature in src_layer: new_feature = ogr.Feature(out_layer.GetLayerDefn()) new_feature.SetGeometry(feature.GetGeometryRef().Clone()) # 复制原有字段 for i in range(src_defn.GetFieldCount()): field_name = src_defn.GetFieldDefn(i).GetName() new_feature.SetField(field_name, feature.GetField(field_name)) # 更新状态 if feature.GetField("STATUS") == "待更新": new_feature.SetField("STATUS", "已更新") # 标记大型地块 area = feature.GetField("AREA") if area is not None and area > 5000: new_feature.SetField("LEVEL", "大型") else: new_feature.SetField("LEVEL", "普通") out_layer.CreateFeature(new_feature) new_feature = None out_layer.CommitTransaction() # 关闭数据集 src_dataset = None out_dataset = None print("处理完成")

这个案例覆盖了读取、新建、更新三个核心操作。实际跑的时候，注意SetField的字段名要和CreateField时一致，大小写敏感。另外Clone()几何那一步不能省，否则源数据的几何对象会被输出图层引用，关闭时可能出问题。

8. 关于坐标系转换的一个补充技巧

前面提到面积计算需要投影转换，这里再补充一个动态转换的实用方法。如果你经常需要在不同坐标系之间切换，可以封装一个转换函数：

def transform_geometry(geom, source_epsg, target_epsg): source_srs = osr.SpatialReference() source_srs.ImportFromEPSG(source_epsg) target_srs = osr.SpatialReference() target_srs.ImportFromEPSG(target_epsg) # 保持坐标轴顺序一致，避免经纬度颠倒 source_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) transform = osr.CoordinateTransformation(source_srs, target_srs) new_geom = geom.Clone() new_geom.Transform(transform) return new_geom

SetAxisMappingStrategy这行很关键。GDAL 3.x默认遵循EPSG定义的坐标轴顺序，有些坐标系是纬度在前、经度在后，和传统GIS习惯相反。不设置这个的话，转换出来的坐标可能经纬度颠倒，数据跑到地球另一边去。这个坑我在处理某批数据时踩过，排查了一整天才定位到原因。

实际用的时候，如果源数据已经有坐标系信息，可以直接从layer.GetSpatialRef()获取，不用手动指定EPSG。只有在源数据坐标系缺失或错误时，才需要手动指定。

9. 写在最后的一些个人习惯

用GDAL Python处理shp这些年，我养成了几个习惯，分享出来供参考。

第一，任何写操作之前先备份。shp的更新是原地修改，一旦写错很难恢复。我通常会在代码里加一行shutil.copy把原文件复制一份带时间戳的备份，成本极低但能救命。

第二，处理完的数据一定要用QGIS打开看一眼。代码跑通不代表数据正确，坐标系、编码、几何有效性这些问题，可视化检查是最直观的。我见过太多次代码没报错但数据打开是空白的情况。

第三，字段名尽量用英文。虽然shp支持中文字段名，但跨平台、跨软件的时候容易出问题。用英文加注释的方式，可维护性更好。

第四，善用ogr.UseExceptions()。这个前面提过，但值得再强调一次。GDAL的静默失败机制是新手最大的敌人，开启异常能让问题暴露得更早。

GDAL的Python绑定确实不算友好，API设计带着浓厚的C++风格，文档也比较分散。但它的能力是实打实的，尤其是在批量处理和精细控制方面，目前还没有能完全替代它的Python库。把读取、新建、更新这三个核心操作吃透，大部分shp相关的数据处理任务都能覆盖。剩下的就是多练、多踩坑、多总结，慢慢就形成自己的套路了。