Edited at

Pythonを用いたベクトル形式データのラスタ化と描画


はじめに

地理情報データは,ベクトル形式(ライン,ポリゴンなど)データとラスタ(グリッドベース)データに大別される.ベクトルデータを用いたラスタデータのマスキングやフィルタリングをおこなうに際し,ベクトルデータをラスタデータとしてあらかじめ変換しておくと処理が楽な場合がある.ということで,今回はrasterioモジュールを用いて,geopandasで読み込んだベクトルデータをラスタデータに変換する.

また,覚書的に,ベクトルデータとラスタデータを同じ図の上に描画するコードについてもまとめた.


使用モジュール

Python3.7系にて操作を行った.

import numpy as np

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from functools import partial
import pyproj
from shapely.ops import transform
import geopandas as gpd
import rasterio as rasio
import rasterio.features as rasioftr


世界地図ポリゴンの読込

今回は,geopandasモジュールで一般公開されている世界地図データを使用する.

変換する解像度は1度×1度(180グリッド×360グリッド)とした.

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

print(world.head())

v2r_table1.png


変換行列の設定

rasterioモジュールの変換ツールを用いて,緯度経度および配列サイズ情報を用いて変換行列を作成する.東西南北の順番および緯度経度長さの順番に注意.

nlatlon = (180, 360)

WSEN_ar = (-180, -90, 180, 90)
WSEN_trans = rasio.transform.from_bounds(*WSEN_ar, *nlatlon[::-1])
print(WSEN_trans)
'''
| 1.00, 0.00,-180.00|
| 0.00,-1.00, 90.00|
| 0.00, 0.00, 1.00|
'''


既存のカラム値を挿入してラスタ化

世界地図データには,各地域の人口とGDPのデータが付与されている.ラスタ化する際にいずれの値を保持するかを決定する必要がある.


各地域の人口

rasterio.featureクラス下のrasterize()関数を使用して,shapely形式のベクトルデータをnumpy arrayに変換する.関数内の使用可能な変数についてはリンクを参照のこと.下記では,基本的なアイテム(shapesout_shapetransform)に加え,ポリゴン外の配列値fill,境界部分の設定all_touchを設定している.ここでも,東西南北の順番および緯度経度長さの順番に注意.

item = 'pop_est'

shps = [(s, i) for i, s in zip(world[item], world['geometry'])]
result = rasioftr.rasterize(
shps,
out_shape=nlatlon,
transform=WSEN_trans,
fill=np.nan,
all_touched=True)
plt.imshow(result)
plt.colorbar(orientation='horizontal')

v2r_thumb_pop.png


各地域のGDP(百万ドル)

人口とまったく同様に変換を行う.

item = 'gdp_md_est'

shps = [(s, i) for i, s in zip(world[item], world['geometry'])]
result = rasioftr.rasterize(
shps,
out_shape=nlatlon,
transform=WSEN_trans,
fill=np.nan,
all_touched=True)
plt.imshow(result)
plt.colorbar(orientation='horizontal')

v2r_thumb_gdp.png


新規カラムを挿入してラスタ化

先述のアイテムに加えて,geopandas内の関数で計算できるポリゴン面積と,文字列データカラムを利用した大陸州ごとのラスタ化を実施した.


各地域の面積

この記事を基にして計算を行った.

def CalcArea(geom):

geom_swap = Swap_xy(geom)
area = geom_to_area(geom_swap)
return area * 1e-6 # m2 -> km2

world['area_km'] = world['geometry'].map(CalcArea)
print(world.head())

v2r_table2.png

あとは全く同様.

v2r_thumb_area.png


各地域の大陸州ごとに値を変えたラスタ化

まず,新たなカラムを加える.

world['new_id'] = 0

print(world.head())

v2r_table3.png

次に,既存のカラムcontinentの項目ごとに数値を割り振る.

for i, name in enumerate(world['continent'].unique()):

print(i, name)
world.loc[world['continent'] == name, 'new_id'] = i + 1
print(world.head())

v2r_table4.png

あとは同様.

v2r_thumb_cont.png


南→北向きにラスタ化

上記までのラスタ化では,緯度次元(0次元)が北から南向きとなっている.しかしながら,緯度が増える方向(南から北向き)のほうが都合のいいことが多く,その際は下記のように東西南北の順番を変更する必要がある.

WNES_ar = (-180,  90, 180, -90)

WNES_trans = rasio.transform.from_bounds(*WNES_ar, *nlatlon[::-1])

この変換行列を用いて描画した図は,上と全く同様にimshow()を用いて描画すると上下が反転する.関数内でorigin='lower'を付け加えることで南→北向きの正しい方向の図を描画できる.

v2r_thumb_cont_s2n.png


ポリゴンとラスタの描画

ここからは,作成したラスタデータと元のベクトルデータとを描画するコードを示す.

ここでも,東西南北の順番に注意.

WESN_ar = (-180, 180, -90,  90)

WESN_w = (-180, 180, -90, 90)
# WESN_w = ( 130, 150, 30, 50)

fig = plt.figure(figsize=(14, 7))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.set_extent(WESN_w, crs=ccrs.PlateCarree())
im = ax.imshow(result, extent=WESN_ar, origin='upper') #, origin='lower')
world.plot(ax=ax, facecolor='none', edgecolor='k')

v2r_glob_cont.png

上記のWESN_w(ax.set_extent()にて使用)を置き換えると,描画範囲を変更することができる.

v2r_japan_cont.png


グリッド中心を基準としたラスタ化

最後に,rasterize()においてall_point=Falseとした際の図を示す.グリッドの中心がポリゴン内に含まれるか否かで判定がなされるため,上図と比較してマスクされた範囲が多くなる.

v2r_japan_cont_untouch.png