1. 程式人生 > >(資料科學學習手札65)利用Python實現Shp格式向GeoJSON的轉換

(資料科學學習手札65)利用Python實現Shp格式向GeoJSON的轉換

一、簡介

  Shp格式是GIS中非常重要的資料格式,主要在Arcgis中使用,但在進行很多基於網頁的空間資料視覺化時,通常只接受GeoJSON格式的資料,眾所周知JSON(JavaScript Object Nonation)是利用鍵值對+巢狀來表示資料的一種格式,以其輕量、易解析的優點,被廣泛使用與各種領域,而GeoJSON就是指在一套規定的語法規則下用JSON格式儲存向量資料,本文就將針對GeoJSON的語法規則,以及如何利用Python完成Shp格式到GeoJSON格式的轉換進行介紹。

 

二、Shp轉GeoJSON

2.1 GeoJSON格式說明

  GeoJSON本質依舊是JSON,其基本格式如下:

{
  "type": "FeatureCollection",
  "features": []
}

  一個完整的GeoJSON檔案最外層為一個字典,把整個GeoJSON檔案看做自頂向下的樹狀結構的話,其根目錄包含鍵值對"type":"FeaturesCollection",以及存放所有要素的鍵值對"features":[],所有向量要素都存放在這個列表中,每個要素都是一個字典,下面我們來認識一下各種向量要素在GeoJSON中的規範格式:

點要素(Point):

  對於單個點要素,其格式如下:

{"type":"Feature",
    "properties":{value1,value2},
    "geometry":{
        "type":"Point",
        "coordinates":[經度,緯度]
        }
    }

  其中properties對應的值為這個要素對應的屬性表中按順序存放的值,geometry對應的值中type指明瞭要素型別,coordinates傳入一個包含兩個元素的列表,第一個元素代表經度,第二個元素代表緯度。

多點要素(MultiPoint):

  多點要素是點要素的特殊情況,其geometry下的type屬性傳入"MultiPoint",其coordinates屬性傳入的是一個二維列表,其最內層列表定義了每個點的經緯度,如下:

{"type":"Feature",
    "properties":{value1,value2},
    "geometry":{
        "type":"MultiPoint",
        "coordinates":[[經度1,緯度1],
                [經度2,緯度2]
            ]
        }
        }

線要素(LineString):

  線要素記錄的是一條線上所有折點的經緯度資訊,只需要按順序連線這些折點就可以還原一條線的形態,在GeoJSON中線要素與多點要素在coordinates屬性上格式相同,區別在於geometry屬性需要傳入"LineString",如下:

{"type":"Feature",
    "properties":{value1,value2},
    "geometry":{
        "type":"LineString",
        "coordinates":[[經度1,緯度1],
        [經度2,緯度2],
        [經度3,緯度3],
        [經度4,,緯度4]]
        }
    }

多線要素(MultiLineString):

  多線要素是多個線要素的組合,因此其coordinates傳入三維列表,來組合多條線,對應的geometry下type屬性為"MultiLineString",如下:

{"type":"Feature",
    "properties":{value1,value2},
    "geometry":{
        "type":"MultiLineString",
        "coordinates":
        [
            [
                [經度1,緯度1],
                [經度2,緯度2],
                [經度3,緯度3],
                [經度4,緯度4]
            ],
            [
                [經度5,緯度5],
                [經度6,緯度6]
            ]
        ]
                }
    }

多邊形要素(Polygon):

  多邊形要素記錄了構成一個多邊形所有邊緣折點的經緯度資訊,其coordinates屬性傳入"Polygon",其geometry下type屬性格式為三維列表,其第三層列表中巢狀的所有列表記錄的經緯度按順序連線即構成了一個多邊形,但需要注意的是,多邊形頭尾折點的經緯度需要相同,才能構成一個閉合的多邊形,如下:

{"type":"Feature",
    "properties":{value1,value2},
    "geometry":{
        "type":"Polygon",
        "coordinates":[
                        [
                          [經度1,緯度1],
                          [經度2,緯度2],
                          [經度3,緯度3],
                          [經度4,緯度4],
                          [經度1,緯度1]
                        ]
                      ]
        }
    }

多多邊形要素(MultiPolygon):

  多多邊形的格式為四維列表,其geometry下type屬性傳入"MultiPloygon",由於多多邊形要素中存在幾種特殊情況,下面我們在geojson.io中進行對應GeoJSON資料的視覺化以便於理解:

  互不重疊的兩個多邊形:

  下面是互不重疊的兩個多邊形的示例:

  對應的GeoJSON資料如下:

{
  "type": "Feature",
  "properties": {},
  "geometry": {
  "type": "MultiPolygon",
  "coordinates":
    [ 
        [
            [
                [102.74414062499999,36.217687122250574],
                [102.7001953125,35.585851593232356],
                [104.8590087890625,35.496456056584165],
                [104.96337890625,36.24427318493909],
                [102.74414062499999,36.217687122250574]
            ]
        ],
        [
            [
                [102.6397705078125,35.074964853989556],
                [103.0352783203125,34.23905366851639],
                [105.00732421875,34.24813554589752],
                [105.3973388671875,35.77771427205079],
                [104.556884765625,35.05698043137265],
                [102.711181640625,35.16931803601131],
                [102.6397705078125,35.074964853989556]
            ]
        ]
    ]
             }
}

  可以看到在多個多邊形不重疊時,直接將兩個多邊形要素對應的三維列表存放在最外層列表下即可。

  互有重疊的兩個多邊形:

  互有重疊的多個多邊形要素格式同多個不重疊的多邊形,效果如下:

  對應的GeoJSON資料如下:

{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "MultiPolygon",
    "coordinates": [
      [
        [
          [101.6455078125,27.68352808378776],
          [114.78515624999999,27.68352808378776],
          [114.78515624999999, 35.209721645221386],
          [101.6455078125,35.209721645221386],
          [101.6455078125,27.68352808378776]
        ]
      ],
      [
        [
          [104.2822265625,30.107117887092357],
          [108.896484375,30.107117887092357],
          [108.896484375,33.76088200086917],
          [104.2822265625,33.76088200086917],
          [104.2822265625,30.107117887092357]
        ]
      ]
    ]
  }
}

  有孔的多邊形:

  有孔的多邊形在類別上也是歸類到MultiPolygon,下面是一個示例:

  對應的GeoJSON資料如下,可以看出其與多個重疊的多邊形的區別在於多邊形向量資訊巢狀在第二層列表中:

{
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPolygon",
        "coordinates":
    [ 
        [
            [
                [101.6455078125,27.68352808378776],
                [114.78515624999999,27.68352808378776],
                [114.78515624999999,35.209721645221386],
                [101.6455078125,35.209721645221386],
                [101.6455078125,27.68352808378776]
            ],
            [
                [104.2822265625,30.107117887092357],
                [108.896484375,30.107117887092357],
                [108.896484375,33.76088200086917],
                [104.2822265625,33.76088200086917],
                [104.2822265625,30.107117887092357]
            ]
        ]
    ]
  }
}

 


2.2 將Shp格式轉換為GeoJSON

  在2.1中我們較為詳細的瞭解到向量資料在GeoJSON資料中具體的表現形式,通過下面的自編函式,以Shp檔名稱(去除檔案拓展名)、Shp檔案編碼、GeoJSON檔案編碼為輸入引數:

def Shp2JSON(filename,shp_encoding='utf-8',json_encoding='utf-8'):
    '''
    這個函式用於將shp檔案轉換為GeoJSON檔案
    :param filename: shp檔案對應的檔名(去除檔案拓展名)
    :return:
    '''

    '''建立shp IO連線'''
    reader = shapefile.Reader(filename,encoding=shp_encoding)

    '''提取所有field部分內容'''
    fields = reader.fields[1:]

    '''提取所有field的名稱'''
    field_names = [field[0] for field in fields]

    '''初始化要素列表'''
    buffer = []

    for sr in tqdm(reader.shapeRecords()):
        '''提取每一個向量物件對應的屬性值'''
        record = sr.record

        '''屬性轉換為列表'''
        record = [r.decode('gb2312','ignore') if isinstance(r, bytes)
                  else r for r in record]

        '''對齊屬性與對應數值的鍵值對'''
        atr = dict(zip(field_names, record))

        '''獲取當前向量物件的型別及向量資訊'''
        geom = sr.shape.__geo_interface__

        '''向要素列表追加新物件'''
        buffer.append(dict(type="Feature",
                           geometry=geom,
                           properties=atr))

    '''寫出GeoJSON檔案'''
    geojson = codecs.open(filename + "-geo.json","w", encoding=json_encoding)
    geojson.write(json.dumps({"type":"FeatureCollection",
                              "features":buffer}) + '\n')
    geojson.close()
    print('轉換成功!')

  下面我們通過一個示例來展示實際轉換效果,使用到的Shp資料為中國省份資料,在arcgis中效果如下:

import shapefile
import json
import codecs

def Shp2JSON(filename,shp_encoding='utf-8',json_encoding='utf-8'):
    '''
    這個函式用於將shp檔案轉換為GeoJSON檔案
    :param filename: shp檔案對應的檔名(去除檔案拓展名)
    :return:
    '''

    '''建立shp IO連線'''
    reader = shapefile.Reader(filename,encoding=shp_encoding)

    '''提取所有field部分內容'''
    fields = reader.fields[1:]

    '''提取所有field的名稱'''
    field_names = [field[0] for field in fields]

    '''初始化要素列表'''
    buffer = []

    for sr in tqdm(reader.shapeRecords()):
        '''提取每一個向量物件對應的屬性值'''
        record = sr.record

        '''屬性轉換為列表'''
        record = [r.decode('gb2312','ignore') if isinstance(r, bytes)
                  else r for r in record]

        '''對齊屬性與對應數值的鍵值對'''
        atr = dict(zip(field_names, record))

        '''獲取當前向量物件的型別及向量資訊'''
        geom = sr.shape.__geo_interface__

        '''向要素列表追加新物件'''
        buffer.append(dict(type="Feature",
                           geometry=geom,
                           properties=atr))

    '''寫出GeoJSON檔案'''
    geojson = codecs.open(filename + "-geo.json","w", encoding=json_encoding)
    geojson.write(json.dumps({"type":"FeatureCollection",
                              "features":buffer}) + '\n')
    geojson.close()
    print('轉換成功!')


if __name__ == '__main__':
    import os
    os.chdir(r'C:\Users\hp\Desktop\飛線圖素材')
    Shp2JSON(filename='bou2_4p.shp',
             shp_encoding='gbk',
             json_encoding='utf-8')

  執行之後同一目錄下出現對應的json檔案:

  匯入到Kepler.gl中進行視覺化:

from keplergl import KeplerGl
import json

with open('bou2_4p.shp-geo.json') as b:
    data = json.load(b)

map1 = KeplerGl(height=700,data={'layer1':data});map1

  

 

  以上就是本文的全部內容,如有筆誤望指出!

&n