1. 地理数据格式转换的核心痛点与方案选型搞GIS开发或者数据处理的朋友大概率都遇到过这样的场景从某个开放数据平台下载了一批GeoJSON格式的矢量数据兴冲冲地拖进桌面GIS软件里准备大干一场结果发现软件版本老旧不支持或者数据量一大就卡得让人想砸键盘。又或者你负责的一个后端服务需要把前端上传的GeoJSON统一转成Shapefile格式入库手动操作显然不现实必须写脚本自动化处理。这时候Python生态里的GDAL/OGR库就成了绕不开的选择。GeoJSON和Shapefile这两种格式表面上看都是存矢量地理数据的但骨子里的差异非常大。GeoJSON基于JSON文本格式一个文件就能装下点线面各种几何类型属性字段灵活读写方便特别适合网络传输和前后端交互。Shapefile则是Esri早年定义的一套格式实际上是一组文件的集合至少包含.shp、.shx、.dbf三个文件分别存几何、索引和属性。它的字段名有10个字符的长度限制字段类型也相对固定不支持嵌套结构。正因为这些差异转换过程不是简单的“另存为”中间有不少细节需要处理。这篇文章面向的读者很明确有Python基础、了解GIS基本概念、需要批量或自动化处理地理数据格式转换的开发者。我会从GDAL/OGR的安装配置讲起逐步拆解GeoJSON转Shapefile的完整实现过程包括坐标系处理、字段映射、编码问题、批量转换等核心环节最后分享一些我在实际项目中踩过的坑和总结出来的排查技巧。整套方案基于Python的gdal和ogr模块不依赖任何桌面软件可以在服务器端无人值守运行。2. 环境准备与GDAL/OGR基础认知2.1 GDAL/OGR到底是什么关系很多刚接触的朋友会被GDAL和OGR这两个名字搞晕。简单来说GDAL全称Geospatial Data Abstraction Library最初专注于栅格数据处理OGR则是它旗下的矢量数据处理模块。后来两者合并到同一个代码仓库统称GDAL但矢量相关的API仍然保留OGR的命名。在Python里你通过from osgeo import gdal, ogr, osr导入这三个模块分别对应栅格、矢量、空间参考系统。做GeoJSON转Shapefile主要用的是ogr和osr。理解这一点很重要因为你在搜索资料时会发现有人写gdal.OpenEx()有人写ogr.Open()有人用ogr2ogr命令行本质上都是同一套底层库的不同调用方式。Python绑定的优势在于可以把转换逻辑嵌入到更复杂的业务流程里比如从数据库读数据、转换后直接上传到对象存储这些是命令行工具做不到的。2.2 安装GDAL的正确姿势GDAL的安装是新手遇到的第一道坎。它底层是C写的Python绑定需要编译好的二进制包。在Windows上最省心的方式是用conda安装conda install -c conda-forge gdalconda-forge渠道的GDAL包维护得很好版本更新也及时。如果你坚持用pip在Linux上通常没问题pip install GDAL$(gdal-config --version)注意这里要匹配系统已安装的GDAL版本否则编译会报错。Windows上pip安装GDAL比较折腾需要预编译的whl文件而且要和Python版本、系统架构严格对应。我的建议是除非你有特殊需求否则直接用conda管理地理空间相关的Python环境能省下大量排查依赖的时间。安装完成后验证一下from osgeo import gdal, ogr, osr print(gdal.__version__)能正常打印版本号就说明环境没问题了。如果报ImportError: DLL load failed大概率是环境变量PATH里缺少GDAL的bin目录或者conda环境没激活。2.3 为什么选GDAL/OGR而不是其他库Python里处理地理数据的库不少比如geopandas、fiona、shapely。geopandas确实好用但它底层读写文件还是靠fionafiona底层又是GDAL。也就是说你绕了一圈还是在用GDAL只是多了一层封装。对于简单的转换任务geopandas几行代码就能搞定但一旦遇到字段类型复杂、编码特殊、需要精细控制转换参数的情况直接操作GDAL/OGR反而更灵活。另一个考虑是部署环境。geopandas依赖pandas、shapely、fiona等一系列库打包体积大在一些资源受限的环境里不太友好。GDAL/OGR虽然安装麻烦一点但一旦装好运行时依赖很干净。我在一个边缘计算项目里就因为这个原因放弃了geopandas改用纯GDAL/OGR实现最终打包出来的可执行文件小了很多。3. GeoJSON与Shapefile格式差异深度解析3.1 几何类型的表达差异GeoJSON用一套统一的JSON结构表达所有几何类型一个FeatureCollection里可以混合Point、LineString、Polygon、MultiPoint、MultiLineString、MultiPolygon甚至GeometryCollection。Shapefile则要求一个文件里所有要素的几何类型必须一致要么全是点要么全是线要么全是面。这就带来一个直接问题如果你的GeoJSON里混合了几何类型转成单个Shapefile是行不通的必须按几何类型拆分。我在实际项目中遇到过一种情况从某个API获取的POI数据大部分是点但偶尔有几个被标记成了Polygon比如一个园区用面表示。直接转换会报错解决方案是先遍历所有要素按几何类型分组分别输出成不同的Shapefile。这个逻辑后面会给出代码实现。3.2 属性字段的映射规则Shapefile的.dbf文件对字段有严格限制字段名最长10个字符只支持Character、Date、Numeric、Logical等几种类型不支持数组、对象等复杂结构。GeoJSON的properties是自由的JSON对象字段名可以很长值可以是字符串、数字、布尔、null甚至嵌套对象。转换时OGR会自动做映射但有几个坑需要注意。第一超过10字符的字段名会被截断如果两个字段前10个字符相同就会冲突。第二中文属性值需要确保编码正确Shapefile默认用UTF-8或者系统编码处理不当会乱码。第三null值在Shapefile里通常转成空字符串或0具体行为取决于字段类型。这些细节在批量转换前最好先做一轮数据清洗。3.3 坐标系与投影处理GeoJSON规范推荐使用WGS84地理坐标系EPSG:4326但实际数据不一定遵守。Shapefile本身不强制坐标系但配套的.prj文件记录了投影信息。转换时如果源数据有坐标系定义OGR会自动写入.prj如果没有生成的Shapefile就没有投影信息后续在GIS软件里加载时可能位置不对。我处理过一批从某地图服务导出的GeoJSON坐标是Web墨卡托EPSG:3857但文件里没有任何坐标系声明。直接转成Shapefile后在QGIS里加载发现位置偏到了南极。解决办法是在转换时显式指定源坐标系和目标坐标系用osr模块创建SpatialReference对象配合CoordinateTransformation做投影转换。这个环节是新手最容易忽略的也是导致“数据转出来了但位置不对”这类问题的根本原因。4. 完整转换流程与核心代码实现4.1 基础转换单个GeoJSON文件转Shapefile先来看最基础的场景一个GeoJSON文件转成一个Shapefile。核心思路是用ogr.Open()打开源文件获取图层然后创建目标数据源把图层复制过去。from osgeo import ogr, osr import os def geojson_to_shp_basic(geojson_path, shp_path): # 打开源文件 src_ds ogr.Open(geojson_path, 0) # 0表示只读 if src_ds is None: raise FileNotFoundError(f无法打开文件: {geojson_path}) src_layer src_ds.GetLayer(0) src_srs src_layer.GetSpatialRef() # 创建目标数据源ESRI Shapefile驱动 driver ogr.GetDriverByName(ESRI Shapefile) if os.path.exists(shp_path): driver.DeleteDataSource(shp_path) dst_ds driver.CreateDataSource(shp_path) # 创建目标图层几何类型与源图层一致 geom_type src_layer.GetGeomType() dst_layer dst_ds.CreateLayer( os.path.basename(shp_path).replace(.shp, ), src_srs, geom_type ) # 复制字段定义 src_layer_defn src_layer.GetLayerDefn() for i in range(src_layer_defn.GetFieldCount()): field_defn src_layer_defn.GetFieldDefn(i) dst_layer.CreateField(field_defn) # 逐要素复制 dst_layer_defn dst_layer.GetLayerDefn() for src_feature in src_layer: dst_feature ogr.Feature(dst_layer_defn) dst_feature.SetGeometry(src_feature.GetGeometryRef().Clone()) for i in range(dst_layer_defn.GetFieldCount()): dst_feature.SetField(i, src_feature.GetField(i)) dst_layer.CreateFeature(dst_feature) dst_feature None src_ds None dst_ds None print(f转换完成: {shp_path})这段代码看起来不复杂但有几个关键点值得展开。ogr.Open()的第二个参数0表示只读模式对于转换任务来说足够了避免误写源文件。CreateDataSource()如果目标路径已存在会报错所以先判断再删除或者用DeleteDataSource()清理。字段复制时直接复用源字段定义OGR会自动处理类型映射和名称截断。4.2 处理混合几何类型的拆分逻辑前面提到GeoJSON里混合几何类型时不能直接转成单个Shapefile。解决方案是先扫描所有要素按几何类型分组然后每组输出一个Shapefile。def geojson_to_shp_split(geojson_path, output_dir): src_ds ogr.Open(geojson_path, 0) src_layer src_ds.GetLayer(0) src_srs src_layer.GetSpatialRef() # 按几何类型分组收集要素 geom_groups {} for feature in src_layer: geom feature.GetGeometryRef() if geom is None: continue geom_type geom.GetGeometryType() if geom_type not in geom_groups: geom_groups[geom_type] [] geom_groups[geom_type].append(feature.Clone()) driver ogr.GetDriverByName(ESRI Shapefile) base_name os.path.splitext(os.path.basename(geojson_path))[0] for geom_type, features in geom_groups.items(): type_name ogr.GeometryTypeToName(geom_type).replace( , _) shp_path os.path.join(output_dir, f{base_name}_{type_name}.shp) if os.path.exists(shp_path): driver.DeleteDataSource(shp_path) dst_ds driver.CreateDataSource(shp_path) dst_layer dst_ds.CreateLayer( f{base_name}_{type_name}, src_srs, geom_type ) # 用第一个要素的字段定义创建字段 src_defn src_layer.GetLayerDefn() for i in range(src_defn.GetFieldCount()): dst_layer.CreateField(src_defn.GetFieldDefn(i)) dst_defn dst_layer.GetLayerDefn() for feat in features: dst_feat ogr.Feature(dst_defn) dst_feat.SetGeometry(feat.GetGeometryRef().Clone()) for i in range(dst_defn.GetFieldCount()): dst_feat.SetField(i, feat.GetField(i)) dst_layer.CreateFeature(dst_feat) dst_feat None dst_ds None print(f输出: {shp_path}, 要素数: {len(features)}) src_ds None这里用feature.Clone()把要素复制到内存列表里避免遍历完源图层后要素失效。ogr.GeometryTypeToName()把几何类型枚举转成可读名称方便拼文件名。实际使用时如果数据量特别大全部加载到内存可能吃不消可以改成两遍扫描第一遍统计类型第二遍分别写入。4.3 坐标系转换的完整实现当源数据的坐标系与目标需求不一致时需要在转换过程中做投影变换。下面这个例子把WGS84的GeoJSON转成Web墨卡托的Shapefiledef geojson_to_shp_reproject(geojson_path, shp_path, src_epsg4326, dst_epsg3857): src_ds ogr.Open(geojson_path, 0) src_layer src_ds.GetLayer(0) # 定义源和目标坐标系 src_srs osr.SpatialReference() src_srs.ImportFromEPSG(src_epsg) src_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) dst_srs osr.SpatialReference() dst_srs.ImportFromEPSG(dst_epsg) dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) # 创建坐标转换对象 transform osr.CoordinateTransformation(src_srs, dst_srs) driver ogr.GetDriverByName(ESRI Shapefile) if os.path.exists(shp_path): driver.DeleteDataSource(shp_path) dst_ds driver.CreateDataSource(shp_path) geom_type src_layer.GetGeomType() dst_layer dst_ds.CreateLayer( os.path.basename(shp_path).replace(.shp, ), dst_srs, geom_type ) src_defn src_layer.GetLayerDefn() for i in range(src_defn.GetFieldCount()): dst_layer.CreateField(src_defn.GetFieldDefn(i)) dst_defn dst_layer.GetLayerDefn() for src_feat in src_layer: geom src_feat.GetGeometryRef().Clone() geom.Transform(transform) dst_feat ogr.Feature(dst_defn) dst_feat.SetGeometry(geom) for i in range(dst_defn.GetFieldCount()): dst_feat.SetField(i, src_feat.GetField(i)) dst_layer.CreateFeature(dst_feat) dst_feat None src_ds None dst_ds NoneSetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)这行很关键。GDAL 3.0之后默认遵循EPSG定义的坐标轴顺序而传统GIS软件习惯经度在前纬度在后。不加这行转换出来的坐标可能会颠倒。这个坑我在升级GDAL版本后遇到过排查了半天才发现是轴顺序的问题。4.4 批量转换与性能优化实际项目里往往需要处理成百上千个文件逐个转换效率太低。下面是一个批量转换的框架加入了简单的并行处理import glob from concurrent.futures import ProcessPoolExecutor, as_completed def batch_convert(input_dir, output_dir, pattern*.geojson, workers4): os.makedirs(output_dir, exist_okTrue) files glob.glob(os.path.join(input_dir, pattern)) if not files: print(未找到匹配的文件) return def convert_one(fpath): base os.path.splitext(os.path.basename(fpath))[0] out_path os.path.join(output_dir, base .shp) try: geojson_to_shp_basic(fpath, out_path) return (fpath, True, None) except Exception as e: return (fpath, False, str(e)) with ProcessPoolExecutor(max_workersworkers) as executor: futures {executor.submit(convert_one, f): f for f in files} for future in as_completed(futures): fpath, success, error future.result() if success: print(f成功: {fpath}) else: print(f失败: {fpath}, 原因: {error})用ProcessPoolExecutor而不是线程池是因为GDAL的Python绑定在释放GIL方面表现一般多进程能更好地利用多核。workers数量建议设为CPU核心数的一半到全部具体取决于文件大小和磁盘IO。如果文件特别大并行反而会因为磁盘争抢变慢这时候串行加进度条可能更合适。5. 常见问题排查与避坑经验实录5.1 中文乱码问题的根源与解决中文乱码是Shapefile转换里最常遇到的问题。表现是属性表里的中文变成问号或者乱码。根源在于.dbf文件的编码声明。OGR创建Shapefile时默认使用系统编码在中文Windows上可能是GBK在Linux上可能是UTF-8。如果源GeoJSON是UTF-8目标环境期望GBK就会出问题。解决方案是在创建数据源时显式指定编码driver ogr.GetDriverByName(ESRI Shapefile) dst_ds driver.CreateDataSource(shp_path) # 设置编码选项 dst_ds.SetMetadataItem(ENCODING, UTF-8)或者在创建图层时通过选项传入dst_layer dst_ds.CreateLayer(name, srs, geom_type, options[ENCODINGUTF-8])实测下来最稳妥的做法是统一用UTF-8然后在读取端也确保用UTF-8解析。如果下游系统强制要求GBK可以在写入前把字符串转成GBK编码的bytes再设置但这样比较绕不如在数据流转的早期就统一编码标准。5.2 字段名截断与冲突的处理Shapefile字段名10字符限制导致的截断和冲突往往在转换完成后才发现。比如population_density和population_total都会被截断成population后者覆盖前者数据就丢了。预防措施是在转换前先检查字段名对超长的做重命名。可以写一个映射表def sanitize_field_names(layer_defn): name_map {} used_names set() for i in range(layer_defn.GetFieldCount()): original layer_defn.GetFieldDefn(i).GetName() # 截断到10字符 truncated original[:10] # 处理冲突 if truncated in used_names: suffix 1 while f{truncated[:8]}{suffix} in used_names: suffix 1 truncated f{truncated[:8]}{suffix} used_names.add(truncated) name_map[original] truncated return name_map然后在创建目标字段时用映射后的名字。这个逻辑虽然简单但能避免很多数据丢失的隐患。我在一个项目中就因为没做这个处理导致一批人口数据转换后字段互相覆盖最后不得不重新跑一遍。5.3 大文件内存溢出的应对策略GeoJSON文件动辄几百MB甚至上GB一次性加载到内存容易OOM。OGR的图层遍历本身是流式的但如果你像前面拆分几何类型那样把所有要素Clone到列表里内存就会暴涨。应对策略是分块处理。对于不需要拆分的场景直接流式遍历写入不要缓存要素。对于需要拆分的场景可以先用GetFeatureCount()和GetGeomType()快速统计如果发现只有一种几何类型就走流式路径如果有多种再考虑分块读取。另外GDAL有一些配置选项可以控制内存使用from osgeo import gdal gdal.SetConfigOption(OGR_GEOJSON_MAX_OBJ_SIZE, 0) # 不限制单个对象大小 gdal.SetConfigOption(OGR_GEOJSON_FEATURE_LIMIT, 0) # 不限制要素数这些选项在读取超大GeoJSON时有用但也要注意别把限制设得太宽松导致内存失控。实际项目中我一般会先检查文件大小超过500MB的就考虑用流式解析库如ijson预处理或者用ogr2ogr命令行工具做初步转换。5.4 坐标系缺失导致的定位偏移前面提到过坐标系问题这里再展开说一下排查方法。转换完成后如果发现数据在GIS软件里位置不对先检查生成的.prj文件是否存在、内容是否正确。可以用ogrinfo命令查看ogrinfo -al -so output.shp输出里会显示图层的坐标系信息。如果没有说明源数据没有坐标系定义需要在转换时手动指定。如果.prj内容不对可能是源数据的坐标系声明有误需要根据数据实际情况修正。一个实用的技巧是拿到任何GeoJSON先用文本编辑器打开看看有没有crs字段。虽然GeoJSON规范推荐WGS84但很多数据源会显式声明其他坐标系。如果没有crs字段默认按WGS84处理但如果数据实际是其他坐标系就会出问题。这时候需要结合数据来源和坐标数值范围来判断比如经度超过180的肯定是投影坐标系。5.5 常见问题速查表问题现象可能原因排查方法解决方案转换报错“无法创建数据源”目标路径已存在同名文件检查输出目录先删除或重命名属性表中文乱码编码不一致查看.dbf编码统一用UTF-8字段值丢失字段名截断冲突检查字段名长度转换前重命名字段数据位置偏移坐标系缺失或错误检查.prj文件显式指定坐标系内存溢出一次性加载过多要素监控内存使用流式处理或分块几何类型不匹配GeoJSON混合几何类型检查要素几何类型按类型拆分输出转换速度慢单线程处理大量文件查看CPU利用率多进程并行6. 进阶技巧与生产环境实践6.1 用ogr2ogr命令行做快速验证虽然这篇文章主打Python实现但ogr2ogr命令行工具在调试阶段非常有用。它可以一行命令完成转换适合快速验证数据是否正常ogr2ogr -f ESRI Shapefile output.shp input.geojson如果命令行能转成功而Python脚本报错说明问题出在代码逻辑上而不是数据本身。反过来如果命令行也失败那就是数据格式的问题需要先清洗数据。我习惯在写正式脚本前先用ogr2ogr跑一遍确认数据没有明显问题。6.2 字段类型转换的精细控制OGR默认的字段类型映射有时候不符合需求。比如GeoJSON里的整数在Shapefile里可能被映射成Real浮点导致精度问题。可以在创建字段时手动指定类型from osgeo import ogr field_defn ogr.FieldDefn(population, ogr.OFTInteger) dst_layer.CreateField(field_defn)对于字符串字段可以指定宽度field_defn ogr.FieldDefn(name, ogr.OFTString) field_defn.SetWidth(100) dst_layer.CreateField(field_defn)宽度设置很重要Shapefile的字符串字段默认宽度可能只有80超长会被截断。如果属性值里有长文本一定要显式设置足够大的宽度。6.3 生产环境的日志与监控在服务器上跑批量转换任务日志和监控必不可少。建议记录每个文件的转换状态、耗时、要素数量方便出问题时定位。可以用Python的logging模块import logging logging.basicConfig( levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s, handlers[ logging.FileHandler(convert.log), logging.StreamHandler() ] )在转换函数里加入耗时统计import time start time.time() # ... 转换逻辑 ... elapsed time.time() - start logging.info(f转换完成: {shp_path}, 耗时: {elapsed:.2f}秒)如果任务量很大还可以接入更完善的监控系统记录成功率、失败原因分布等指标。这些数据在优化转换流程时很有参考价值。6.4 输出结果的验证方法转换完成后不能只看文件生成了就完事还要验证内容是否正确。最简单的验证是要素数量对比def verify_conversion(geojson_path, shp_path): src_ds ogr.Open(geojson_path, 0) dst_ds ogr.Open(shp_path, 0) src_count src_ds.GetLayer(0).GetFeatureCount() dst_count dst_ds.GetLayer(0).GetFeatureCount() if src_count dst_count: print(f要素数量一致: {src_count}) else: print(f数量不一致! 源: {src_count}, 目标: {dst_count}) src_ds None dst_ds None更严格的验证可以抽样对比几何坐标和属性值。我在一个项目中就遇到过要素数量一致但部分几何为空的情况原因是源数据里有无效几何转换时被跳过了。所以除了数量最好再检查一下几何有效性。7. 个人实操心得与建议这套GeoJSON转Shapefile的方案我在多个项目里反复用过从单文件转换到每天处理上万文件的批处理管道都跑过。最大的体会是格式转换本身不难难的是处理各种边界情况和数据质量问题。源数据的坐标系声明、字段命名规范、几何有效性、编码一致性这些才是决定转换能否顺利的关键。我的建议是在写转换脚本之前先花时间了解你的数据。用ogrinfo看看图层结构用文本编辑器抽查几个要素确认坐标系和字段情况。这一步花十分钟能省下后面几个小时的排查时间。另外转换脚本一定要有完善的错误处理和日志记录生产环境里最怕的就是静默失败数据没转成功但没人发现。最后分享一个小技巧如果目标系统对Shapefile的字段名有严格要求可以在转换前先用pandas或纯Python对GeoJSON的properties做一轮清洗把字段名规范化、把不需要的字段删掉、把类型统一。这样转换脚本本身可以保持简洁数据清洗的逻辑单独维护职责更清晰。