読者です 読者をやめる 読者になる 読者になる

Met Post

Python、プログラミングのメモ、気象、海洋のこととか

GeoPandasで国土数値情報を加工してみた 2

引き続き、地理情報データをPandasで扱うための拡張ライブラリgeopandasを利用して国土数値情報データを加工してみた。

metpost.hatenablog.com

引き続き、国土数値情報データのダウンロードページの土地利用区分の3次メッシュデータを利用する。今回は、2次メッシュに変換したときの土地利用区分を求めてみる。

まずは、データの読み込みを行う。

# 処理したいshapeファイルのリスト
shpfiles = [shpfile for shpfile in glob.glob('./*/*jgd*/*.shp')]

# shapeファイルを読み込んで結合させる
df = pd.concat([gpd.read_file(shpfile, encoding='shift-jis') for shpfile in shpfiles])

# 後の処理のため3次メッシュを明示的に整数にしておく
df['メッシュ'] = df['メッシュ'].apply(int)

前回は、3次メッシュのまま表示させていたので気にならなかったが、このままだと3次メッシュを結合させる際にPolygonの境界が小数点以下でずれていることが問題となる。元データのせいなのか、Pythonでの読み込みの際の型変換が問題なのかはわからないが、隣接格子で微妙に座標がずれてしまいPolygonを連結させることができなかった。

f:id:ttmych:20170420005129p:plain

これを修正するために、Polygonの座標を四捨五入して小数点以下を整えることにする。以下のエントリを参考にした。

gis.stackexchange.com

from shapely.geometry import shape, mapping
import numpy as np

def set_precision_shape(shape_in, precision):
    try:
        geojson = mapping(shape_in)
        coords = geojson['coordinates']
        geojson['coordinates'] = np.round(np.array(geojson['coordinates']), precision)
    except: 
        return shape_in
    return shape(geojson)

df['geometry'] = df['geometry'].apply(lambda g: set_precision_shape(g, 6))

f:id:ttmych:20170420005804p:plain

次に2次メッシュでの最大値を求める。geometryの集約については、geopandas.dissolveで処理できる。

# 2次メッシュに変換
df['メッシュ'] = df['メッシュ'].apply(lambda x: x//100)

# geometryを集約
df_geo = df[['メッシュ','geometry']].dissolve(by='メッシュ')
# geometry以外の要素の合計を求める
df_attr = df.groupby('メッシュ').sum()
# メッシュコードをキーにして結合させる
df_mesh2 = pd.concat([df_geo, df_attr], axis=1).reset_index()

f:id:ttmych:20170420010057p:plain

あとは前回と同じく、メッシュ毎の最大の利用区分を求めてshapeファイルに出力する。

cols = [e for e in df_mesh2.columns.tolist() if e not in ['geometry', 'メッシュ']]
df_mesh2['landtype'] = df_mesh2[cols].apply(lambda x: x.argmax(), axis=1)
df_mesh2 = df_mesh2[['メッシュ','landtype','geometry']]

QGISで表示すると以下のようになった。2次メッシュへの変換はうまく計算できているようだ。

f:id:ttmych:20170420012838p:plain:w500