1. GIS数据格式转换的核心场景与痛点在土地管理、城乡规划、自然资源调查等实际工作中GIS数据流转就像不同语言国家之间的外交官会谈——需要专业翻译官完成格式转换。我处理过某省国土三调项目就遇到过县级单位提交的SHP文件需要批量转成省级平台要求的GDB格式同时还要生成Excel统计报表的典型场景。坐标系问题是最常见的拦路虎。去年有个项目甲方提供的CAD图纸采用北京54坐标系而我们的分析系统要求2000国家大地坐标系。直接转换会导致地块位置偏移300多米最后通过EPSG代码4524→4490的转换才解决。建议转换前先用gdalsrsinfo命令检查原始数据坐标系gdalsrsinfo input.shp -o proj4批量处理效率是另一个痛点。传统ArcGIS手动操作处理500个SHP文件需要整天时间而用Python脚本配合GDAL库20分钟就能完成。我曾写过一个自动化脚本主要用到以下核心函数import os from osgeo import ogr def batch_shp_to_gdb(input_folder, output_gdb): driver ogr.GetDriverByName(FileGDB) gdb_ds driver.CreateDataSource(output_gdb) for shp in os.listdir(input_folder): if shp.endswith(.shp): shp_ds ogr.Open(os.path.join(input_folder, shp)) layer shp_ds.GetLayer() gdb_ds.CopyLayer(layer, shp[:-4])2. SHP与文本格式的深度互转技巧2.1 SHP转TXT的字段映射实战土地报备场景中经常需要将SHP属性转为固定格式的TXT。有个关键细节字段顺序决定输出结构。比如耕地占补平衡项目要求前三个字段必须是地块编号、地块面积、地类编码。我推荐先用QGIS查看字段元数据layer iface.activeLayer() for field in layer.fields(): print(field.name(), field.typeName())编码问题曾让我踩过大坑。某次转换后的TXT在Windows记事本显示乱码最后发现是UTF-8 BOM编码问题。解决方案是在Python写入时明确编码格式with open(output.txt, w, encodingutf-8-sig) as f: f.write(\ufeff) # 添加BOM头2.2 TXT转SHP的坐标处理玄机测绘行业提供的TXT通常使用Y,X坐标顺序即纬度在前而GIS软件默认是X,Y。有次转换后地块跑到非洲去了就是因为没注意这个细节。建议转换时用以下参数校验# GDAL创建图层时显式指定坐标系 srs osr.SpatialReference() srs.ImportFromEPSG(4524) # 以CGCS2000为例 out_layer out_ds.CreateLayer( parcels, srs, ogr.wkbPolygon, [ENCODINGUTF-8] # 指定中文编码 )对于带号处理8位坐标值前两位就是带号。比如38512345.678,3456789.123表示38度带。可以用PyProj进行带号提取和投影转换from pyproj import Transformer transformer Transformer.from_crs(EPSG:4524, EPSG:4490) x, y transformer.transform(38512345.678, 3456789.123)3. CAD与GIS的跨界转换方案3.1 DXF到SHP的图层保留策略建筑设计院提供的CAD图纸通常包含几十个图层转换时需要选择性过滤。我习惯先用FME Workbench预览DXF结构再用以下OGR命令提取特定图层ogr2ogr -f ESRI Shapefile output.shp input.dxf -sql SELECT * FROM layers WHERE OGR_GEOMETRYLINESTRING文字标注处理是个技术活。某次转换后所有中文标注变成问号原因是CAD使用了特殊字体。解决方案是先用AutoCAD将文字导出为多行文本属性再在转换时指定编码options [ENCODINGGB18030] # 支持生僻字的编码 driver.CreateDataSource(output.shp, optionsoptions)3.2 坐标系匹配的黄金法则遇到过最棘手的情况是CAD文件没有坐标系信息但甲方坚持说肯定是2000坐标系。后来通过对比图纸上的控制点坐标发现其实是地方独立坐标系。现在我的工作流程必定包含以下步骤用CAD的LIST命令查询至少三个已知点的图纸坐标获取这些点的实际大地坐标通过RTK测量或已有数据使用Ground Control Points进行配准# 使用GDAL进行GCP校正 gcp_points [ (图纸X1, 图纸Y1, 大地X1, 大地Y1), (图纸X2, 图纸Y2, 大地X2, 大地Y2) ] ds.SetGCPs([gdal.GCP(*p) for p in gcp_points], srs.ExportToWkt())4. 批量处理与性能优化之道4.1 海量SHP合并的避坑指南合并上千个SHP时内存管理至关重要。有次合并3000个宗地图直接让32GB内存的服务器崩溃后来改用分块处理方案# 分批次合并策略 chunk_size 100 for i in range(0, len(shp_files), chunk_size): batch_files shp_files[i:ichunk_size] vrt gdal.BuildVRT(ftemp_{i}.vrt, batch_files) gdal.Translate(fbatch_{i}.shp, vrt) del vrt # 显式释放内存属性字段统一是另一个常见问题。合并时建议先用OGR检查所有文件的字段结构ogrmerge.py -o merged.gdb -single -overwrite_ds input/*.shp -field_strategy FirstLayer -src_layer_field_name source_file4.2 Excel与GDB的高效互转国土部门的报表经常要求Excel与GDB双向同步。对于大数据量导出传统ArcGIS的Excel驱动有6万行限制可以用OpenPyXL突破from openpyxl import Workbook wb Workbook(write_onlyTrue) ws wb.create_sheet() for row in cursor: # 使用游标分批读取 ws.append(row) wb.save(big_data.xlsx)几何对象处理需要特别注意。将带坐标的Excel转回SHP时WKT格式比分开的X/Y列更可靠# 生成WKT格式的几何字段 df[geometry] df.apply( lambda r: fPOLYGON(({r.x1} {r.y1},{r.x2} {r.y2},...)), axis1 )在某个智慧城市项目中我们开发了自动化质检模块用PyShp库直接读取SHP文件进行拓扑检查import shapefile sf shapefile.Reader(parcels.shp) for shape in sf.shapes(): if not shape.isValid: print(f无效图形在记录{sf.record(shape.oid)})