python怎么读shape文件?
一、Shapefile基础概念
Shapefile是GIS领域标准矢量格式,由ESRI开发,包含三个核心文件:
- .shp:存储几何要素(点/线/面)
- .shx:几何索引文件
- .dbf:属性数据表
需保证同名文件在同一目录(如roads.shp
,roads.shx
,roads.dbf
)。
二、Python读取方案详解
1. pyshp库(轻量级首选)
import shapefilesf = shapefile.Reader("data/rivers.shp") # 读取文件# 获取元数据
print(sf.shapeType) # 几何类型(1=点,3=线,5=面)
print(sf.bbox) # 地理范围 [minX, minY, maxX, maxY]# 遍历几何要素
shapes = sf.shapes()
for i, shape in enumerate(shapes):points = shape.points # 坐标列表 [(x1,y1), (x2,y2)...]parts = shape.parts # 多部件索引(如岛屿多边形)# 示例:打印第一条线的首点坐标if i == 0: print("首点坐标:", points[0])# 读取属性表
records = sf.records()
for rec in records:print(rec["NAME"], rec["LENGTH_KM"]) # 字段名需实际存在
2. geopandas库(数据分析推荐)
import geopandas as gpdgdf = gpd.read_file("data/countries.shp") # 自动解析几何+属性# 核心操作
print(gdf.head()) # 查看前5行数据
print(gdf.crs) # 坐标系(如EPSG:4326)
gdf.plot() # 自动绘制地图# 属性查询
asia = gdf[gdf["CONTINENT"] == "Asia"] # 筛选亚洲国家
asia.to_file("asia.gpkg", driver="GPKG") # 导出为GeoPackage
3. GDAL/OGR库(高性能专业级)
from osgeo import ogrds = ogr.Open("data/lakes.shp")
layer = ds.GetLayer(0) # 获取第一层# 遍历要素
for feat in layer:geom = feat.GetGeometryRef() # 几何对象name = feat.GetField("NAME") if geom.GetGeometryType() == ogr.wkbPolygon:area = geom.Area() # 计算多边形面积print(f"{name}: {area:.2f} km²")# 读取空间参考
spatial_ref = layer.GetSpatialRef()
print(spatial_ref.ExportToWkt()) # 输出WKT格式坐标系
三、关键问题解决方案
-
中文路径/乱码处理
# pyshp sf = shapefile.Reader("中文/道路.shp", encoding="gbk")# GDAL gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
-
批量读取大型文件
# 分块读取降低内存 with shapefile.Reader("large_file.shp") as sf:for i in range(0, len(sf), 1000): # 每次处理1000条batch = sf.shapeRecords()[i:i+1000]for rec in batch:process(rec.shape, rec.record)
-
坐标系转换
# geopandas示例 gdf = gdf.to_crs("EPSG:3857") # 转Web墨卡托投影
四、工具对比与选型建议
库 | 安装命令 | 适用场景 | 优势 |
---|---|---|---|
pyshp | pip install pyshp | 轻量读取/基础操作 | 无依赖、API简单 |
geopandas | pip install geopandas | 空间分析/数据清洗 | 集成pandas,支持空间运算 |
GDAL | pip install GDAL | 专业GIS开发/坐标系转换 | 支持格式多,功能强大 |
💡 选型建议:
- 快速查看数据 → pyshp
- 数据分析/可视化 → geopandas
- 坐标转换/格式处理 → GDAL
五、进阶技巧
-
写入Shapefile
# pyshp写入示例 w = shapefile.Writer("new_data.shp", shapeType=shapefile.POLYLINE) w.field("ROAD_ID", "N") # 添加数字字段 w.line([[[102, 35], [103, 36]]]) # 添加线段 w.record(101) # 对应属性 w.close()
-
空间计算
# geopandas计算缓冲区 rivers = gpd.read_file("rivers.shp") buffered = rivers.geometry.buffer(0.01) # 创建1km缓冲区
-
可视化强化
import matplotlib.pyplot as plt fig, ax = plt.subplots() gdf.plot(ax=ax, column="GDP", legend=True, cmap="viridis") plt.savefig("map.png", dpi=300)