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())

変換行列の設定

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に変換する．関数内の使用可能な変数についてはリンクを参照のこと．下記では，基本的なアイテム（shapes，out_shape，transform）に加え，ポリゴン外の配列値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')

各地域の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')

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

先述のアイテムに加えて，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())

あとは全く同様．

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

まず，新たなカラムを加える．

world['new_id'] = 0
print(world.head())

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

for i, name in enumerate(world['continent'].unique()):
    print(i, name)
    world.loc[world['continent'] == name, 'new_id'] = i + 1
print(world.head())

あとは同様．

南→北向きにラスタ化

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

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

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

ポリゴンとラスタの描画

ここからは，作成したラスタデータと元のベクトルデータとを描画するコードを示す．
ここでも，東西南北の順番に注意．

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

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

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

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