pythonを用いたshapefileやgeojsonの読込および描画

[Sep. 2019] 古い情報をいろいろと更新

初めに

ポリゴンや点群などの地理情報関連データをpython環境で扱うために，これまでに様々なモジュールが開発されてきた．モジュールごとに機能や出力形式などは様々であり，統一性はほとんどとられていない．ゆえに，理解が難しい．そこで，この記事を書くに至った．本記事はファイルの読込から描画までを主に対象としている．shapefileの新規作成や書込は対象としていない．

開発環境はpython2.7系およびpython3系．

分類

本記事では，各モジュールを次のように分類してみた．

• 読込
  • json，geojson
  • pyshp (shapefile)
  • fiona
  • cartopy (←pyshp，fionaの力を借りている)
  • geopandas
• 処理
  • shapely
  • geopandas (←shapelyの力を借りている)
• 描画
  • matplotlib + descartes
  • geopandas (+ geoplot)
  • cartopy (+ geopandas)
  • folium (+ json, pandas)

--------------------------------------------------------------------------------------------

json，geojson（読込）

(geo)json形式のテキストファイルの読込．geojsonというモジュールもあるが，こちらは主にgeojsonファイル生成に強く，読込だけを考えればjsonモジュールだけで十分．

読込

import json
with open("file.json") as f:
    src = json.load(f) # 辞書型で取得

中身の確認（一例．中身の構造に依存する）

print(src.keys())               # 辞書インデックス確認（['type', 'features']）
ftr = src.get('features')[n]    # n番目のフィーチャー
print(len(src.get('features'))) # フィーチャーの総数
print(ftr.keys())               # 辞書インデックス確認（['geometry', 'type', 'properties']）
geom = ftr.get('geometry')      # geometry情報の取得
prop = ftr.get('properties')    # property情報の取得
print(prop.keys())              # propertyタグ一覧

--------------------------------------------------------------------------------------------

pyshp (shapefile)（読込）

とにかくシンプルである．しかし，表示させることのできる情報も最低限のもののみである．読込で用いるReader()関数にはencodingオプションが付随しており，'SHIFT-JIS'で記述された日本語データの読込も可能（python2系のルールとして，print関数の使用時はutf-8でencodeする必要がある）．

読込

開く

import shapefile

# ファイル読込（encoding='SHIFT-JIS'とすることで日本語データにも対応）
src = shapefile.Reader('shapefile.shp',encoding='SHIFT-JIS')

フィーチャ選択

shp = src.shape(n)              # n番目のフィーチャー（geometry座標タプルのリスト）
rec = src.record(n)             # n番目のフィーチャー（property情報リスト）
print(rec0[i])                  # i番目のproperty情報. py2系で.encode('utf-8')必要

# リストとして取り出すこともできる
shp = src.shapes()[n]
rec = src.records()[n]

# 2つを合成したものもある
shprecs = src.shapeRecords()
shp = shprecs[n].shape
rec = shprecs[n].record

# 一度に全部を読み込んでデータが重くなる場合，イテレータ機能を用いる．
for shprec in src.iterShapeRecords():
    shp = shprec.shape
    rec = shprec.record
:   :   :

中身の確認

print(src.bbox)             # 全体の緯度経度矩形情報（タプル）
print(src.numRecords)       # フィーチャーの総数
print(shp.points)           # フィーチャー境界線(経度緯度タプルのリスト)
print(shp.shapeTypeName)    # フィーチャーの種別

GeoJSON形式で取得

print(shp.__geo_interface__)

--------------------------------------------------------------------------------------------

fiona（読込）

関数等がシンプルな点で，pyshpのshapefileと似ている．fionaの大きな特徴として，open()がイテレータとして機能する．そのため，next(src)，もしくはsrc.next()とすることによって，頭から次々とフィーチャーを取得することができる．フィルタリングをするような際も，一般にfiona.open()で開いたソースをforループに入れて処理を行うことが推奨されているようである．

読込

開く

import fiona

# ファイル読込（encoding='SHIFT-JIS'とすることで日本語データにも対応）
src = fiona.open('shapefile.shp',encoding='SHIFT-JIS')

中身の確認

# 注:一度len(src)などを実施するとイテレータが末尾に到達してしまう．その際はsrcを再読込．
print(src.bounds)                   # 全体の緯度経度矩形情報（タプル）
print(len(src))                     # フィーチャーの総数
shp0 = next(src)                    # 1つ出力される（geojson形式の辞書型）
# shp0 = src.next()                 # 同上
geom = shp0.get('geometry')         # geometry情報の取得
prop = shp0.get('properties')       # property情報の取得
print(prop.keys())                  # propertyタグ一覧

フィルタリング

src.filter(bbox={rect_tuple})   # 空間フィルタリング

shapelyの構造で抽出（フィルタリング込み）

from shapely.geometry import Polygon, MultiPolygon, shape

mp = MultiPolygon(
    [shape(ftr['geometry']) for ftr in fiona.open('shapefile.shp')
    if ftr['properties']['property'] == 'item'])

--------------------------------------------------------------------------------------------

geopandas（読込，処理，描画）

csvのようなテーブル形式のデータ処理を行うpandasモジュールベースのパッケージ．
pandas.DataFrame形式で読込を行うことができ，pandasモジュールで定義される便利なメソッドを活用できる．後述するshapely.geometryモジュールの仕様が組み込まれており，ポリゴンの面積や重心などを表示させることができる．基礎的な描画にも対応しており，非常に高度な技術を求めない限り，本モジュール一つで十分な機能が備わっている．

読込

開く

import geopandas as gpd

# ファイル読込（encoding='SHIFT-JIS'とすることで日本語データにも対応）
df_shp = gpd.read_file('shapefile.shp',encoding='SHIFT-JIS')
df_gjs = gpd.read_file('gjs.geojson',encoding='SHIFT-JIS')

フィーチャ選択

ftr = df_shp.iloc[n]                    # n番目のフィーチャー
print(ftr.{property}.encode('utf-8'))   # propertyの値確認
geom = ftr.geometry                     # geometry情報の取出

中身の確認

print(len(df_shp.index))    # フィーチャー総数
print(df_shp.head())        # 上位数行表示
print(df_shp.columns        # propertyタグ一覧(geometry含む)
print(df_shp.total_bounds)  # 全体の緯度経度矩形情報（numpy）
print(df_shp.count())       # 各propertyごとの，値が入っているフィーチャー総数
# (下記メソッドは上のgeomに対しても適用化)
print(df_shp.geom_type)     # 各フィーチャーの種別
print(df_shp.area)          # 各フィーチャーの面積
print(df_shp.length)        # 各フィーチャーの周長
print(df_shp.centroid)      # 各フィーチャーの重心

フィルタリング

df_valid = df_shp[df_shp.{property}.notnull()]  # nullの排除
filtered = df_shp[df_shp.{property}=='item']    # フィルタリング(詳細はpandasを参照)
filtered = df_shp[df_shp.geometry.within(geom)] # 空間フィルタリング

処理

GeoJSON形式で取得

print(df.to_json())

shapelyの構造で抽出

print(type(df.geometry[n])) # ポリゴン<shapely.geometry.polygon.Polygon>
print(type(df.exterior[n])) # 境界線<shapely.geometry.polygon.LinearRing>

描画

# axはmatplotlib.axesインスタンス．指定しない場合は新規figureに描画される．
ax = plt.axes()
df.plot(ax=ax,column='property')        # propertyの値に基づくカラーリング
df.plot(ax=ax,color='r')                # 単色描画
df.boundary.plot(ax=ax,edgecolor='k')   # 境界線のみ描画

参照：geoplot（描画モジュール）

geopandasで開いたデータをより多様に描画する．→ リンク ←を参照．

--------------------------------------------------------------------------------------------

cartopy（読込，描画）

cartopyは，投影法などを指定して地図情報の描画が可能なモジュールである．
そのioモジュールを用いてshapefileを読み込むことができる．それ以外には，geopandasとcartopyを組み合わせて様々な投影法を用いた描画を行うといったこともできる．

shapereader()で読み込んだデータはgeometry及びrecordそれぞれに対するジェネレータを有しており，fionaと非常に似ている．というのも，読込のベースとしてpyshp系のshapefileまたはfionaを採用しているため当然である．しかしながら，cartopyパッケージ内での（特に描画に適した）用途のために，一部機能を制約・変更している．例えば，geometry情報で取り出されるのが，後述するshapely.geometry形式のポリゴンである点がfionaなどと異なる点である．また，読み込み時のencodingオプションを引き継いでいない．
どうしてもencodingをもちいたい，といった場合には，ソースコードを簡単に書きかえればよい．

（参考）encoding対応：site-packages/cartopy/io/shapereader.py

    :   :   :
class BasicReader(object):
    def __init__(self, filename, encoding='utf-8'): # ←追加点
        # Validate the filename/shapefile
        self._reader = reader = shapefile.Reader(filename, encoding=encoding) # ←追加点
    :   :   :
class FionaReader(object):
    def __init__(self, filename, bbox=None, encoding=None): # ←追加点
        self._data = []

        with fiona.open(filename, encoding=encoding) as f: # ←追加点
    :   :   :

読込

from cartopy.io import shapereader

src = shapereader.Reader('shapefile.shp')
# src = shapereader.Reader('shapefile.shp', encoding='SHIFT-JIS') ←上記変更実施時
for geom,rec in (src.geometries(), src.records()):
    :   :   :

描画

慣れていないとcartopyの描画については少し癖が強い．
投影法（crs）の設定は必須である．迷ったらとりあえずccrs.PlateCarree()を使う．
また，座標の値など，自動で描画してくれるアイテムが全くないので，下記コードを参考にアレンジするか，geoviewsのような描画サポートパッケージを使用する．

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.feature import ShapelyFeature

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent("""(tuple_rect)""", ccrs.PlateCarree())
ax.coastlines(resolution='10m')
gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
shp_ftr = ShapelyFeature(src.geometries(),
                                ccrs.PlateCarree(), facecolor='r')
ax.add_feature(shp_ftr)

--------------------------------------------------------------------------------------------

ogr (from osgeo)（読込）

地理情報系で有名なosgeoパッケージのモジュール．ogrはOpenGIS Simple Features Reference Implementationの略称らしい．図の描画はJSON形式の文字列などで出力することにより別のモジュールに任せることとなるが，shapefileなどのファイルの読込や書込のために非常に高度な処理を施すことができる点では非常に便利である．同じosgeo下のgdalと組み合わせれば，rasterデータへの変換などもできる（→事例←）．
※デフォルトのエンコーダがutf-8であるため，ExportToJson()などではSHIFT-JISで書かれた日本語データの読込がうまくいかない．

読込

開く

from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
src = driver.Open('shapefile.shp')  # 0は読取専用の意味
# src = ogr.Open('shapefile.shp') でも可
lyr  = src.GetLayer(0)              # 通常意識しなくてよい1文．

フィーチャ選択

ftr  = lyr.GetFeature(n)            # n番目のフィーチャー
print(ftr.keys())                   # 対象フィーチャーの各アイテムの変数
print(ftr.items())                  # 対象フィーチャーの各アイテム（辞書型）
print(ftr.GetField('property').decode('shift-jis').encode('utf-8')) # （2系）propertyの値確認（日本語であればdecodeする）
geom = ftr.geometry())              # geometry情報の取出
print(geom.ExportToJson())          # str型でJSON形式取得(Polygon)
# 下記，2系では日本語データがあるとUnicodeDecodeErrorが発生．
print(ftr.ExportToJson())           # str型でJSON形式取得(Polygon)

中身の確認

lyr_dfn = lyr.GetLayerDefn()
print([lyr_dfn.GetFieldDefn(i).GetName() for i in range(lyr_dfn.GetFieldCount())]) # 項目リスト表示(geometry除く)
print(lyr.GetExtent())                  # 全体の緯度経度矩形情報（タプル）
print(lyr.GetFeatureCount())            # フィーチャーの総数
print(geom.Area())                      # フィーチャー面積
print(geom.Boundary().ExportToJson())   # フィーチャー境界線(LineString)
print(geom.Centroid().ExportToJson())   # フィーチャーの重心(Point)

フィルタリング（注：lyr内のフィーチャー群が上書きされる）

lyr.SetAttributeFilter("{property} = 'item'")   # 項目でフィルタリング
lyr.SetSpatialFilter(geom)                      # 空間でフィルタリング

--------------------------------------------------------------------------------------------

shapely（格納，処理）

平面幾何図形を処理する代表的なモジュール．格納のためには，GeoJSON形式のデータから抽出したgeometry情報が必要になる．重心や面積の計算ができるようになるほか，他のフィーチャーとの結合や，union処理なども可能となる．処理についての説明は省略する．

格納

from shapely.geometry import shape
# from shapely.geometry import asShape
from shapely.geometry.polygon import Polygon

geom = geoj_file.get('geometry')
shp_geom = shape(geom)
# shp_geom = asShape(geom) でも可

中身の確認

print(shp_geom.geom_type)   # フィーチャーの種別
print(shp_geom.area)        # フィーチャーの面積
print(shp_geom.bounds)      # フィーチャーの矩形情報（タプル）
print(shp_geom.length)      # フィーチャーの周長
print(shp_geom.centroid.xy) # フィーチャーの重心
print(shp_geom.exterior.xy) # フィーチャーの外周のx座標リスト，y座標リスト

--------------------------------------------------------------------------------------------

matplotlib + descartes（描画）

descartesは，matplotlib.axesに対してpatchとしてポリゴンなどを貼り付けることを容易に可能とするモジュールである．maplotlibにおける描画の工夫もいろいろ考えられる．ここでは取り上げていないが，basemapを用いて既存のポリゴンの力を借りて描画することも大事である．

descartes（描画）

import matplotlib.pyplot as plt 
from descartes import PolygonPatch

geom = geoj_file.get('geometry')
poly = PolygonPatch(geom, fc='r', ec='k', alpha=0.5)

ax.add_patch(poly)

ax.plot()を用いて境界線（描画）

xs, ys = shp_geom.exterior.xy
ax.plot(xs,ys,c='k')

参考：PatchCollectionによる一括描画

from matplotlib.collections import PatchCollection
from descartes import PolygonPatch
from shapely.geometry import Polygon, MultiPolygon, shape

multipoly = MultiPolygon(shp_geom_list) # shapeのリストを用いてMultiPolygonの作成

def colorscheme(poly): # 描画色決定用関数
    :   :   :

patches = [] # patchのリスト
for poly in multipoly:
    fcolor = colorscheme(poly) #描画色指定
    polypatch = PolygonPatch(poly, fc=fcolor, ec='k', zorder=1)
    patches.append(polypatch)
patchcoll = PatchCollection(patches, match_original=True) # patchのリストを用いてPatchCollectionの作成
ax.add_collection(patchcoll)
ax.autoscale()

--------------------------------------------------------------------------------------------

folium （描画）

地理座標系データの描画ということで，最後に実際の地図への描画をするためのfoliumモジュールの紹介．geojsonファイルであればダイレクトに読み込んで図示することができる．

GeoJSONファイルを用いたシンプルな描画

下記コードでは単色のポリゴン描画がなされる．
別途geometryに対応するpropertyテーブルを用意して，style_functionを編集してやれば，propertyの値を用いた配色など，自由度が大きくなる．

import folium

m = folium.Map(
    location={coords}, 
    zoom_start=7, 
    width=800, 
    height=600, 
    tiles='openstreetmap') # マップ作成
folium.GeoJson('gjs.geojson', name='region_name',
              style_function = lambda x: {
                                        'fillOpacity': 0.5,
                                        'fillColor': 'Orange',
                                        'color': 'Red'
              }).add_to(m)

folium.LayerControl().add_to(m) # 地図表示切替ボタン
m.save('svname.html') # 保存

pandasモジュールと組み合わせた階級区分（Choropleth）図

dataでpandas.DataFrameを指定，columnsはのDataFrameの中の2列(key，値)をタプルで指定．key_onにはGeoJSONにおけるプロパティをfeature.id，もしくは'feature.properties.xxx'として入力する．2つのデータソースのキーがそろっていることが重要．

folium.Choropleth(geo_path='gjs.geojson', data=df,
             columns=['key', 'property'],
             key_on='feature.properties.xxx',
             threshold_scale=[100, 200, 300, 400, 500],
             fill_color='colormap').add_to(m)

終わりに

pandasに慣れている人であれば，基本的にgeopandasを選択しておけば事足りるはず．しかしながら，少し複雑な図を描きたいと思ったとき，上記のパッケージをいくつか組み合わせたような事例が散見されるから大変．参照しているコードが，どのパッケージのクラスを用いているのかに気をつけながら，自由度の高い描画を行えればよいと思う．