在最近的项目中，完成了许多python处理矢量数据的算法程序，比如缓冲区分析、叠置分析、统计分析等，主要用到的是GDAL/OGR库，很多功能都参照了此链接中的示例：http://pcjericks.github.io/py-gdalogr-cookbook/#
由于完成的功能在上面链接中都能找到实现的代码，所以原来不打算写博客的，但是今天写的这个SHP转GeoJSON的算法把我坑得不轻，着实郁闷了好久才解决，因此就写一下权当纪念吧。有时候写程序真的是一件充满玄学的事情，原理上怎么看怎么对，但是运行起来却偏偏不是那回事，就算调试也给你不报错直接闪退，特别是处理空间数据时属性的编码问题，更是搞得人头皮发麻。这种情况下，能不能调的出来全靠运气，搞不好一个小问题卡一个星期都是正常，这可能也是程序员普遍发量稀少的原因吧。更加玄学的是，解决问题的方式又常常充满戏剧性，折腾了一个早上加半个下午没找到原因，起来喝杯水却突然灵感乍现，三四行代码调整一下顺序，一顿骚操作，点击运行，突然就能运行成功了，再看一下结果，还真是对的！山穷水尽疑无路，柳暗花明又一村，这种感觉真是刺激。
好了，将Shapefile转化为GeoJSON的完整代码如下：
1 # -*- coding: utf-8 -*- 2 from osgeo import ogr 3 import gdal 4 import sys 5 import os 6 7 def ChangeToJson(vector, output): 8 print("Starting........") 9 #打开矢量图层 10 gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") 11 gdal.SetConfigOption("SHAPE_ENCODING", "GBK") 12 shp_ds = ogr.Open(vector) 13 shp_lyr = shp_ds.GetLayer(0) 14 15 # 创建结果Geojson 16 baseName = os.path.basename(output) 17 out_driver = ogr.GetDriverByName('GeoJSON') 18 out_ds = out_driver.CreateDataSource(output) 19 if out_ds.GetLayer(baseName): 20 out_ds.DeleteLayer(baseName) 21 out_lyr = out_ds.CreateLayer(baseName, shp_lyr.GetSpatialRef()) 22 out_lyr.CreateFields(shp_lyr.schema) 23 out_feat = ogr.Feature(out_lyr.GetLayerDefn()) 24 25 #生成结果文件 26 for feature in shp_lyr: 27 out_feat.SetGeometry(feature.geometry()) 28 for j in range(feature.GetFieldCount()): 29 out_feat.SetField(j, feature.GetField(j)) 30 out_lyr.CreateFeature(out_feat) 31 32 del out_ds 33 del shp_ds 34 print("Success........") 35 36 if __name__ == '__main__': 37 shapefile = 'D:/Python37Work/HainanAOI/affected_bridge.shp' 38 out = 'D:/Python37Work/HainanAOI/affected_bridge.json' 39 ChangeToJson(shapefile, out)
有以下几个注意事项：
1、gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")和gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
这两行用于支持中文，前一句用于支持文件名称及路径内的中文，后一句用于支持属性字段中的中文，最好将其放在最前面，否则读取矢量文件时不能按GBK编码方式读取。
2、out_lyr.CreateFields(shp_lyr.schema)和out_feat = ogr.Feature(out_lyr.GetLayerDefn())
这两句的顺序一定不能错，否则会导致部分字段无法写入，应该是先对out_lyr创建相应的字段，然后再由out_feat来一个个生成Feature。
最终生成的.json文件内容如下：
