Python pyshp库实战：Shapefile文件读写与GIS数据处理全解析

1. 项目概述：为什么Shapefile依然是GIS的“活化石”？

如果你在地理信息系统（GIS）、城市规划、环境科学或者数据分析领域工作过，哪怕只是浅尝辄止，也一定绕不开一个文件格式：Shapefile。这个由ESRI公司在90年代初推出的数据格式，以其简单的结构和广泛的兼容性，统治了地理数据交换领域近三十年。时至今日，尽管有GeoJSON、GPKG（GeoPackage）等更现代格式的挑战，Shapefile凭借其“行业普通话”的地位，依然是数据交换、项目协作中最常见、最稳妥的选择。

“pyshp读写shapefile”这个标题，指向的正是用Python处理这个经典格式的核心技能。pyshp（官方库名shapefile）是一个纯Python库，它不依赖GDAL/OGR等庞大的C/C++库，却能完整地读写Shapefile。这意味着你可以在任何Python环境中，轻装上阵地操作.shp（几何图形）、.shx（几何图形索引）、.dbf（属性数据）这一系列文件。对于数据分析师、自动化脚本开发者、以及需要将地理处理流程嵌入Web后端或轻量级应用的工程师来说，pyshp提供了极大的便利。

掌握pyshp，你就能在Python生态里自由地：从政府开放数据平台下载的行政区划数据中提取特定城市的边界；将业务数据（如销售网点、物流轨迹）转换为空间数据进行分析；或者将处理好的地理数据导出，供QGIS、ArcGIS等专业软件进行可视化。它解决的是地理数据“进得来、出得去、处理得了”的基础问题，是空间数据分析工作流中不可或缺的一环。

2. 核心原理与文件结构拆解：Shapefile的“三驾马车”

在动手写代码之前，我们必须先理解Shapefile到底是什么。它不是一个单一文件，而是一个由多个文件组成的集合，每个文件扮演着不同的角色。pyshp库的强大之处，就在于它用纯Python优雅地封装了对这组文件的操作。

2.1 Shapefile的组成文件与角色

一个完整的Shapefile至少包含三个核心文件，它们像三驾马车，共同承载了地理数据：

1. 主文件 (.shp): 存储地理要素的几何图形信息。例如，一个点要素的坐标、一条线的节点序列、一个多边形的环和顶点坐标。它是二进制格式，直接读取是乱码，需要专门的解析器。
2. 索引文件 (.shx): 这是.shp文件的索引。它记录了每个几何图形在.shp文件中的起始位置（偏移量）和长度。有了它，软件可以快速定位和读取特定的图形，而无需遍历整个文件，这对处理大型数据至关重要。
3. 属性文件 (.dbf): 以dBASE IV表格格式存储每个地理要素的属性数据。例如，一个代表城市的多边形，其.shp文件存储边界坐标，而对应的.dbf文件则存储城市名称、人口、GDP等字段和记录。.dbf是早期数据库格式，但因其简单，被广泛支持。

此外，常见的辅助文件还包括：

• .prj: 存储坐标系统信息（如WGS84, CGCS2000）。非常重要，没有它，你的数据只是一堆没有意义的数字坐标。pyshp可以读写此文件，但本身不进行坐标转换。
• .cpg: 可选，用于指定.dbf文件的字符编码（如UTF-8），解决中文乱码问题。
• .sbn/.sbx: 空间索引文件，加速空间查询，通常由GIS软件生成。

pyshp在读取时，你只需要提供主文件名（如'counties.shp'），它会自动寻找同名的其他文件。在写入时，它会一次性生成所有必要的文件。

2.2 pyshp的工作模式：Reader与Writer

pyshp的API设计非常直观，主要围绕两个核心类展开，这与Shapefile的读写分离特性完美对应：

• shapefile.Reader: 用于读取已存在的Shapefile。你可以通过它遍历所有要素（shapeRecords()或iterShapeRecords()），分别获取几何图形（shape）和属性记录（record），也可以读取文件头信息（bbox,shapeType等）。
• shapefile.Writer: 用于创建新的Shapefile。你需要先定义几何类型（shapeType）和属性字段（field），然后通过shape()和record()方法依次添加图形和属性，最后调用save()生成所有文件。

这种设计模式清晰地将数据消费（读）和生产（写）分开，符合大多数数据处理流程的直觉。

注意：一个常见的误解是认为.shp文件包含了所有信息。实际上，.dbf文件同样重要。在pyshp中，几何和属性是紧密关联但分别处理的。当你删除一个要素时，需要确保从图形列表和属性记录列表中同步删除对应的条目，否则会导致数据错位。

3. 从零开始：使用pyshp读取Shapefile全流程

让我们从一个具体的例子开始。假设你从统计部门拿到了一个名为city_boundaries.shp的文件，里面包含了多个城市的边界多边形及其名称、代码。我们的目标是读取它，并筛选出特定人口规模的城市。

3.1 环境准备与库安装

首先，确保你的Python环境（建议3.7以上）已经就绪。安装pyshp非常简单，因为它没有任何二进制依赖：

pip install pyshp

安装完成后，你可以在Python中导入它。库的名称是shapefile，但通常我们为其设置一个简短的别名sf，以方便编码：

import shapefile as sf

3.2 基础读取与数据探查

读取一个Shapefile的第一步是创建Reader对象。

# 假设Shapefile文件位于当前目录，无需添加后缀 reader = sf.Reader('city_boundaries')

创建好reader对象后，我们可以先探查一下这个数据的基本情况，这就像拿到一份新数据先看“元数据”。

# 1. 查看几何类型 print(f"几何类型代码: {reader.shapeType}") # 输出如 5 (代表多边形 Polygon) # shapeType代码含义: 1=点，3=线，5=多边形，8=多点，11=点Z，13=线Z，15=多边形Z等 # 2. 查看空间范围 (边界框) bbox = reader.bbox print(f"数据边界框: {bbox}") # 格式: [最小经度, 最小纬度, 最大经度, 最大纬度] # 3. 查看属性字段定义 fields = reader.fields[1:] # fields的第一个元素是删除标记，通常跳过 for field in fields: print(f"字段名: {field[0]}, 字段类型: {field[1]}, 长度: {field[2]}, 精度: {field[3]}") # 字段类型示例: 'C'表示字符型，'N'表示数值型，'F'表示浮点型，'D'表示日期型 # 4. 查看要素总数 print(f"要素总数: {len(reader)}")

这些信息对于后续处理至关重要。例如，知道了shapeType，你才能正确地理解几何数据；知道了字段定义，你才知道如何正确地提取属性。

3.3 遍历要素与提取数据

最常用的方法是遍历每一个要素，同时获取其几何图形和属性。iterShapeRecords()方法是一个生成器，适合处理大型文件，因为它不会一次性将所有数据加载到内存。

# 用于存储目标城市的信息 target_cities = [] for shape_record in reader.iterShapeRecords(): # shape_record是一个对象，包含 .shape 和 .record 属性 geom = shape_record.shape # 几何对象 attr = shape_record.record # 属性列表，顺序与fields定义一致 # 假设字段定义是: ['CITY_NAME', 'CITY_CODE', 'POPULATION'] city_name = attr[0] population = attr[2] # 注意索引从0开始 # 进行业务逻辑判断，例如筛选人口大于500万的城市 if population and population > 5000000: # 提取几何信息。对于多边形，points包含所有环的所有顶点 # shape.points 返回顶点列表 [[x1,y1], [x2,y2], ...] # shape.parts 指明每个环的起始顶点在points列表中的索引 city_boundary_points = geom.points target_cities.append({ 'name': city_name, 'population': population, 'geometry': city_boundary_points }) print(f"找到大城市: {city_name}, 人口: {population}") # 处理完成后，关闭reader（虽然不是必须，但是好习惯） reader.close()

对于简单的需求，你也可以使用shapeRecords()方法一次性获取所有要素的列表，但请注意数据量过大时可能的内存压力。

3.4 处理常见读取问题：中文乱码与复杂几何

中文乱码问题：这是处理中文数据时最常见的“坑”。Shapefile的.dbf文件默认编码通常是系统本地编码（如gbk），而现代环境多用UTF-8。如果读取时出现乱码，你需要指定编码。

# 方法一：在创建Reader时指定编码（如果存在.cpg文件，pyshp可能会自动识别） try: reader = sf.Reader('city_boundaries', encoding='gbk') # 尝试用gbk编码 except UnicodeDecodeError: reader = sf.Reader('city_boundaries', encoding='utf-8') # 尝试用utf-8编码 # 方法二：更稳妥的方式是，在读取属性后对字符串字段进行解码 for shape_record in reader.iterShapeRecords(): attr = shape_record.record # 假设第一个字段是城市名，是字符串类型 city_name_raw = attr[0] if isinstance(city_name_raw, bytes): # 尝试解码 try: city_name = city_name_raw.decode('gbk') except: city_name = city_name_raw.decode('utf-8', errors='ignore') else: city_name = str(city_name_raw)

复杂几何类型：Shapefile支持带Z值（高程）或M值（测量值）的几何类型（如PointZ, PolyLineM）。pyshp会将这些值存储在shape.z或shape.m列表中。在处理3D数据或路径测量数据时，需要额外关注这些数组。

if reader.shapeType in [11, 13, 15, 18]: # 这些是带Z值的类型 for shape_record in reader.iterShapeRecords(): geom = shape_record.shape points = geom.points # 二维坐标 [ [x,y], ... ] z_values = geom.z # 对应的高程值 [z1, z2, ...] # 处理三维数据...

4. 实战进阶：使用pyshp创建与编辑Shapefile

读懂了数据，下一步就是创造数据。假设我们需要根据业务数据，生成一个全国零售店网点的Shapefile，包含店名、地址和日销售额属性。

4.1 创建新的Shapefile：定义结构与添加数据

创建过程是一个“先搭架子，再填内容”的过程。

import shapefile as sf # 1. 创建Writer对象，并指定几何类型。1代表点（Point） writer = sf.Writer('retail_stores', shapeType=1) # 2. 定义属性字段。field方法的参数：字段名、字段类型、最大长度、小数精度 # 字段类型：'C'=字符，'N'=整数/小数，'F'=浮点，'D'=日期 writer.field('STORE_NAME', 'C', 50) # 店名，字符型，最大50长度 writer.field('ADDRESS', 'C', 100) # 地址，字符型，最大100长度 writer.field('SALES', 'N', 12, 2) # 销售额，数值型，总长12位，小数2位 writer.field('OPEN_DATE', 'D') # 开业日期，日期型 # 3. 添加数据（假设stores_data是一个字典列表） stores_data = [ {'name': '中心旗舰店', 'address': '人民路1号', 'sales': 125000.50, 'date': '2023-05-01'}, {'name': '东区分店', 'address': '创业大道88号', 'sales': 89000.00, 'date': '2022-11-15'}, ] for store in stores_data: # 添加几何图形：.point(x, y, [z], [m]) # 这里需要真实的经纬度坐标，示例中使用虚构值 writer.point(116.4074, 39.9042) # 假设是北京的坐标 # 添加属性记录：.record(*args)，参数的顺序必须与field定义的顺序严格一致！ writer.record(store['name'], store['address'], store['sales'], store['date']) # 4. 保存文件。这一步会生成 retail_stores.shp, .shx, .dbf 等文件 writer.save() print("Shapefile保存成功！")

重要提示：writer.record()的参数顺序必须与之前调用writer.field()的顺序完全一致，否则会导致属性数据错位，这是新手最容易出错的地方之一。建议使用变量名来明确对应关系，或者将数据整理成与字段定义同序的列表。

4.2 设置投影信息（.prj文件）

生成的Shapefile默认没有投影信息。为了让它在GIS软件中正确显示，我们必须创建.prj文件。这需要你知道数据的坐标系统WKID（Well-Known ID）或WKT（Well-Known Text）字符串。

# 方法一：使用epsg.io的代码（推荐，最常用） # 例如，为WGS84经纬度坐标创建.prj文件 prj_content = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]' with open('retail_stores.prj', 'w') as f: f.write(prj_content) # 方法二：使用pyproj库动态生成（更专业） # 首先安装 pip install pyproj from pyproj import CRS crs = CRS.from_epsg(4326) # WGS84 with open('retail_stores.prj', 'w') as f: f.write(crs.to_wkt())

4.3 编辑现有Shapefile：修改与删除

pyshp没有提供直接的“编辑模式”。编辑的思路是：读取 -> 在内存中修改数据 -> 写入一个新文件。这是函数式数据处理中常见的模式。

场景：删除销售额低于某个阈值的店铺，并为剩余店铺添加一个“等级”字段。

import shapefile as sf # 1. 读取原始文件 reader = sf.Reader('retail_stores') shapeType = reader.shapeType fields = reader.fields records = reader.records() shapes = reader.shapes() # 2. 准备新的Writer，并复制原有字段定义 writer = sf.Writer('retail_stores_updated', shapeType=shapeType) for field in fields[1:]: # 跳过第一个删除标记字段 writer.field(*field) # 3. 添加一个新字段 writer.field('RANK', 'C', 10) # 4. 遍历，筛选，并添加新数据 new_records = [] new_shapes = [] for i, (shape, record) in enumerate(zip(shapes, records)): sales = record[2] # 假设销售额是第三个字段 if sales >= 100000: # 筛选条件 new_shapes.append(shape) # 构建新的属性记录：旧字段 + 新字段值 new_record = list(record) if sales > 200000: new_record.append('A') # 添加等级 else: new_record.append('B') new_records.append(new_record) # 5. 将筛选和修改后的数据写入Writer for shape, record in zip(new_shapes, new_records): writer.shape(shape) writer.record(*record) # 6. 保存新文件 writer.save() reader.close()

这种方法本质上是创建了一个全新的数据集。对于大型数据，需要注意内存使用。对于更复杂的编辑（如修改某个图形的顶点），你可以直接操作shape.points列表，然后再用writer.shape()添加。

5. 性能优化与高级技巧：处理大规模数据

当Shapefile包含数十万甚至上百万个要素时，简单的遍历操作可能会变得缓慢。以下是一些提升效率的实战技巧。

5.1 使用迭代器与分块处理

始终优先使用iterShapeRecords()或iterShapes()和iterRecords()，避免一次性将shapes()和records()全部读入内存。

# 好的做法：迭代处理 with sf.Reader('huge_data') as reader: # 使用上下文管理器，确保文件关闭 for sr in reader.iterShapeRecords(): # 处理每个要素 process_feature(sr) # 可以每处理一定数量就保存或输出一次，减少内存峰值

5.2 利用NumPy进行批量几何计算

如果需要对所有图形的坐标进行数学运算（如平移、缩放），将坐标数据转换为NumPy数组会极大提升速度。

import numpy as np import shapefile as sf reader = sf.Reader('data') points_list = [] # 收集所有点图形的坐标 for shape in reader.iterShapes(): if shape.shapeType == 1: # 点 # shape.points 是 [[x, y]] 列表 points_list.append(shape.points[0]) # 转换为NumPy数组 (n_points, 2) points_array = np.array(points_list) # 进行批量运算，例如将所有点向东平移0.01度 points_array[:, 0] += 0.01 # 再写回新的Shapefile（这是一个简化的例子，实际需重建图形对象） writer = sf.Writer('shifted_points', shapeType=1) writer.field('ID', 'N') for i, pt in enumerate(points_array): writer.point(pt[0], pt[1]) writer.record(i) writer.save()

5.3 空间过滤：使用边界框预筛选

如果你只关心某个矩形区域内的数据，可以先利用reader.bbox和每个shape的bbox属性进行快速粗筛，避免对每个图形进行复杂的几何计算。

target_bbox = [115.0, 38.0, 118.0, 41.0] # 目标区域边界框 for shape_record in reader.iterShapeRecords(): shape = shape_record.shape # 图形边界框与目标边界框是否相交（快速判断） if not (shape.bbox[2] < target_bbox[0] or # 图形最右 < 目标最左 shape.bbox[0] > target_bbox[2] or # 图形最左 > 目标最右 shape.bbox[3] < target_bbox[1] or # 图形最上 < 目标最下 shape.bbox[1] > target_bbox[3]): # 图形最下 > 目标最上 # 再进行精确的几何判断（如点是否在多边形内） if precise_intersection_check(shape, target_bbox): process_feature(shape_record)

6. 避坑指南与常见问题排查

在实际使用pyshp的过程中，你肯定会遇到一些意想不到的问题。下面是我从大量实践中总结出的“血泪教训”。

6.1 文件锁定与权限问题

在Windows系统上，如果你用Reader打开了一个文件，但没有关闭它，再去写入或删除这个文件，可能会遇到“权限被占用”的错误。

解决方案：

• 使用上下文管理器：这是最推荐的方式。

  with sf.Reader('data.shp') as reader: # 在此块内操作reader data = list(reader.iterShapeRecords()) # 退出块后，文件自动关闭

• 显式关闭：在 finally 块或处理完成后手动关闭。

  reader = sf.Reader('data.shp') try: # 操作 finally: reader.close()

• 写入时注意：Writer.save()之后，Writer对象的工作就完成了。如果需要再次写入，应创建新的Writer实例。

6.2 几何类型不匹配错误

尝试将线（ShapeType=3）添加到点（ShapeType=1）类型的Writer中，会引发错误。

排查步骤：

1. 打印reader.shapeType确认源数据的几何类型。
2. 创建Writer时，确保shapeType参数与你要写入的数据类型一致。如果你要写入多种类型（通常不建议，Shapefile标准规定一个文件一种类型），需要统一为最复杂的类型（如将点和线都存为“多点”或“多线”，但这会破坏属性关联的直观性）。

6.3 属性数据错位或丢失

这是最高频的问题，症状是：在GIS软件中打开，图形和属性对不上；或者某个字段的值全部显示为None。

原因与解决：

1. 字段顺序不一致：writer.record(a, b, c)中的a, b, c必须与之前writer.field()定义的字段顺序、数量、类型完全匹配。建议使用列表或元组来传递记录值，避免手动输入时出错。

   field_names = ['Name', 'Value'] record_values = ['Test', 100] # 确保field_names和record_values的顺序逻辑一致 writer.record(*record_values)

2. 字段长度不足：定义字段writer.field('NAME', 'C', 5)时，最大长度为5。如果实际字符串'北京市'长度超过5，写入时会被截断或导致错误。在定义字段时，预留足够的长度。
3. 数据类型不匹配：尝试将字符串写入'N'（数值）字段，或反之。确保写入的数据类型与字段定义相符。日期字段需要传入datetime.date对象或符合特定格式的字符串。

6.4 生成的Shapefile在GIS软件中无法打开或显示异常

1. 缺少必要文件：确保.shp,.shx,.dbf三个文件在同一目录下，且主文件名相同。pyshp的save()方法会生成它们。
2. 投影问题：数据在GIS软件中显示的位置不对（如跑到非洲或北极）。检查并正确创建.prj文件。用文本编辑器打开.prj文件，确认其内容是正确的WKT字符串。
3. 几何错误：某些GIS软件对几何图形的有效性检查很严格。例如，多边形不能自相交，环的顶点顺序（外环逆时针、内环顺时针）需符合规范。pyshp本身不检查这些，它“忠实”地记录你给它的顶点。如果遇到显示问题，可能需要用更专业的库（如shapely）进行几何验证和修复。
4. 编码问题：属性中的中文显示为乱码。确保写入时字符串是str类型（Python 3默认unicode）。如果从其他源（如GBK编码的CSV）读取数据，先将其解码为unicode。写入后，可以尝试手动创建一个.cpg文件，里面只写一行UTF-8，并与其他文件放在一起，提示GIS软件使用UTF-8编码打开。

6.5 性能瓶颈排查

当处理速度很慢时：

• 检查循环内部：避免在遍历十万级要素的循环内部进行复杂的文件I/O操作（如频繁打开小文件、打印日志到控制台）。
• 使用分析工具：用Python的cProfile模块分析代码，找到耗时最长的函数。
• 考虑升级方案：对于超大规模（千万级点）的数据，纯Python的pyshp可能力不从心。此时应考虑使用基于C/C++的GDAL/OGR库（通过fiona或ogrPython绑定），或者将数据导入空间数据库（如PostGIS）进行处理。

最后，一个最朴素的建议：在处理重要数据前，先在小样本（如前10个要素）上完整跑通你的读写逻辑，并用QGIS或ArcGIS快速打开检查一下。这能帮你提前发现大部分几何、属性和投影问题，避免批量处理后的返工。pyshp就像一把精准的螺丝刀，在理解Shapefile这套老式但稳固的体系后，它能帮你高效地完成大多数地理数据的基础操作。