🙌

シェープファイル(shpファイル)をPython+Geopandasで扱った

2021/08/21に公開

Python

gis

GeoPandas

tech

背景

都道府県などの基礎的な地理情報は国土地理院などが公式のデータを配布していることがあります．そのようなデータを扱うとき，シェープファイル(*.shp)やGeoJSONなどのファイルを扱うことがあります．今日はPython+Geopandas+Shapelyを利用した手順を自分のためにまとめました (これらの処理のpure Julia版について知っている方は教えてください)．

これらの地理情報ファイルの背景知識や関係性について，例えば

のような記事を見ていただければ関連技術について理解が深まると思います (私は全然詳しくないです)．

利用するデータ

今日はこちらのリンク先で国土交通省が公開している「国土数値情報データ > 行政区域データ」を利用します．特に理由はないですが長野県のデータを使ってみようと思うので，適当な年度の（ここでは令和3年度）長野県データ（N03-20210101_20_GML.zip）をダウンロードしました．

ディレクトリ構造です (note.ipynbに使い捨てのスクリプトを書きます)．shpファイル以外のデータも必要なのでそのまま置いておいてください．

> tree
.
├── N03-20210101_20_GML
│   ├── KS-META-N03-21_20_210101.xml
│   ├── N03-21_20_210101.dbf
│   ├── N03-21_20_210101.geojson
│   ├── N03-21_20_210101.prj
│   ├── N03-21_20_210101.shp
│   ├── N03-21_20_210101.shx
│   └── N03-21_20_210101.xml
├── N03-20210101_20_GML.zip
└── note.ipynb

データの読み込みと可視化

Python+Geopandas+Shapelyを利用した処理をします．

Geopandasを利用した読み込みと可視化

Pandasのデータフレームのような感じで読み込みします．

import geopandas as gpd
import matplotlib.pyplot as plt

path_shp = "./N03-20210101_20_GML/N03-21_20_210101.shp"
gdf = gpd.read_file(path_shp)
gdf.head()

次のようなデータフレームが得られます（都道府県情報，市町村名，ID，ポリゴン情報などです）．

具体的なポリゴンを描画してみます．0行目が長野市らしいので，0行目のポリゴンをmatplotlibで描画します（ただし私は長野市の形を知らないので合っているかどうかは見ても分からないですが，Wikipeidaの長野市のページに貼られている図を見ると雰囲気合っていそうですね）．またipython notebookを利用しているときは，ポリゴンを単純に実行させるだけでセルに表示されるようです．ここではmatplotlibに渡すためにPolygonのexteriorからxy座標列を取得しています．

f = plt.figure(figsize=(6, 6))
a = f.gca()
a.plot(*gdf.iloc[0].geometry.exterior.xy)

次にデータフレームの全てのポリゴンを描画してみます．ここでは手抜きして gdf.plot() を呼び出してnotebookから結果を確認します．これは完全な長野ですね．

点がポリゴンに含まれているか判定する

Geopandasで読み込んだポリゴンの情報はshapely.geometry.polygon.Polygon型のオブジェクトになっています．そのため，点が含まれているかどうかや外面 (上で利用したexterior)，面積などの情報を取得できます．ここではある点 $(lat, lon)$ がポリゴンに含まれているかどうかを判定したいと思います．まずはインターネットで長野市の緯度経度情報を取得してきます．長野市のポリゴンがこの地点を含んでいるかどうかを判定します．

poly0 = gdf.iloc[0].geometry
p_nagano = shapely.geometry.point.Point(138.18111, 36.65139) # 県庁所在地
print(poly0.contains(p_nagano)) # > True

当然のようにTrueが出力されます．

ポリゴンを結合(union)する

国土地理院の行政区画データは各市区町村（？）について用意されていて便利ですが，ここで県を1つにまとめたオブジェクトとして「点が長野県に入っているか？」みたいな処理を実行したいときを考えます．このとき国土地理院のデータから取得したポリゴンをくっつけて1つの県ポリゴンを用意すると便利そうです．このようなポリゴンの和集合処理を実行するには cascaded_union/unary_union という名前の技術を利用します．

一番簡単な適用方法としては，Geopandasが読み込んだデータフレームのgeometryの列がGeoSeriesになっているため，この列に対して直接unary_unionを取得します．

nagano = gdf["geometry"].unary_union

似たような形でプロットさせてみるコードです．

nagano = gdf["geometry"].unary_union
f = plt.figure(figsize=(6, 10))
f.gca().plot(*nagano.exterior.xy)

無事外周が取れました．

その他

何か思い出したら/別の操作が必要になったら書きます．

オマケ | ランダムな点が長野県に含まれるかどうかをプロットした

境界線の $(x, y)$ 情報が得られているため， $x, y$ の範囲が分かっています (min/max)．ここから乱数を生成させ，作成した点 $(x_i, y_i)$ が含まれていれば赤色で，含まれていなければ青色で散布図を描きます．プログラムです．

x, y = nagano.exterior.xy
rangeX = min(x), max(x)
rangeY = min(y), max(y)

trial = 5000
pos, neg = [], []
for n in range(trial):
    x = np.random.uniform(*rangeX)
    y = np.random.uniform(*rangeY)
    flag_isin = nagano.contains(shapely.geometry.point.Point(x, y))
    if flag_isin:
        pos.append([x, y])
    else:
        neg.append([x, y])

pos = np.array(pos)
neg = np.array(neg)

# 描画
f = plt.figure(figsize=(6, 10))
a = f.gca()
a.plot(*nagano.exterior.xy)
a.scatter(pos[:, 0], pos[:, 1], color="r", alpha=0.5)
a.scatter(neg[:, 0], neg[:, 1], color="b", alpha=0.5)
plt.tight_layout()