Edited at

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


初めに

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

開発環境はpython2.7系.

※個人的にpython3系におけるパッケージインストールがうまくいっておらず,←解決.そちらでの挙動は未確認.でき次第,更新予定.


分類

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


  • 読込


    • json,geojson

    • pyshp (shapefile)

    • fiona

    • cartopy

    • geopandas



  • 処理


    • shapely

    • geopandas (←shapelyの力を借りている)



  • 描画


    • matplotlib + descartes

    • geopandas (+ geoplot)

    • cartopy

    • folium



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


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].encode('utf-8') # i番目のproperty情報

# リストとして取り出すこともできる
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'])

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


cartopy(読込,描画)

cartopyは,投影法などを指定して地図情報の描画が可能なモジュールである.そのioモジュールを用いてshapefileを読み込むことができる.

shapereader()で読み込んだデータはgeometry及びrecordそれぞれに対するジェネレータを有しており,fionaと非常に似ている.geometry情報で取り出されるのが,後述するshapely.geometry形式のポリゴンである点がfionaなどと異なる点である.

cartopyではこのshapely形式に対応した描画機能を有している.cartopyの豊富な描画機能およびデフォルトで使用できるポリゴンデータが使えるのは大きな利点である.

※デフォルトのエンコーダがutf-8であるため,src.records().next()においてSHIFT-JISで書かれた日本語データの読込がうまくいかない.


読込

from cartopy.io import shapereader

src = shapereader.Reader('shapefile.shp')
for geom,rec in (src.geometries(), src.records()):
: : :


描画

from cartopy.feature import ShapelyFeature

shp_ftr = ShapelyFeature(src.geometries())
ax.add_feature(shape_feature)

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


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に描画される.

df.plot(ax=ax,column='property') # propertyの値に基づくカラーリング
df.plot(ax=ax,color='r') # 単色描画
df.boundary.plot(ax=ax,edgecolor='k') # 境界線のみ描画


参照:geoplot(描画モジュール)

地図投影法に関する情報を持つcartopyというパッケージとの組合せを目的としたものらしい.geopandasで開いたデータをより多様に描画する.→ リンク ←を参照.

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


ogr (from osgeo)(読込)

地理情報系で有名なosgeoパッケージのモジュール.ogrはOpenGIS Simple Features Reference Implementationの略称らしい.図の描画はJSON形式の文字列などで出力することにより別のモジュールに任せることとなるが,shapefileなどのファイルの読込や書込のために非常に高度な処理を施すことができる点では非常に便利である.同じosgeo下のgdalと組み合わせれば,rasterデータへの変換などもできる(→事例←).

※自分の使用環境および使用データだと,下記でGetGeometryRef()を使用したオブジェクトの処理時にSegmentation Fault (core dump)を起こす現象が発生している.

※また,デフォルトのエンコーダが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.items() # 対象フィーチャーの各アイテム(辞書型)
print ftr.GetField('property').decode('shift-jis').encode('utf-8') # propertyの値確認(日本語であればdecodeする)
geom = ftr.geometry() # geometry情報の取出
print geom.ExportToJson() # str型でJSON形式取得(Polygon)
# 下記,日本語データがあるとUnicodeDecodeErrorが発生.
print ftr.ExportToJson() # str型でJSON形式取得(Polygon)


中身の確認

lyr_dfn = lyr.GetLayerDefn()

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


参考:確認(通常はできるはずだが個人的にCoreDump失敗事例)

# geom = lyr.GetFeature(n).GetGeometryRef()     # n番目のフィーチャー

# geom.Area().ExportToJson() # フィーチャーの面積
# geom.Centroid().ExportToJson() # フィーチャーの重心


フィルタリング(注: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)

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


folium (描画)

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


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

import folium

m = folium.Map(
location={coords},
zoom_start=7,
width=800,
height=600,
tiles='openstreetmap') # マップ作成
m.geo_json(
geo_path='gjs.geojson',
name='region_name'
fill_color='none',
line_color='Orange') # GeoJSON読込(境界線)
folium.GeoJson('gjs.geojson', name='region_name').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つのデータソースのキーがそろっていることが重要.

m.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')


終わりに

日本語が入っているデータを使用する際は,encodingの可否に注意する必要がある(python3での挙動は未確認).pandasに慣れている人であれば,基本的にgeopandasを選択しておけば事足りるはず.それ以外であれば,読込後の解析作業に応じて読込モジュールは使い分けるのがよいと思われる.