在最近的項目中，完成了許多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文件內容如下：
