ラスターデータの操作

Pythonで扱う地理空間データ

Joris Van den Bossche

Open source software developer and teacher, GeoPandas maintainer

値のグリッドの図

画像出典: QGIS ドキュメント

ラスターデータ

降雨確率マップ

複数バンドのラスターデータ

rasterio パッケージ

import rasterio

GDAL への「Pythonic」バインディング
ラスターファイルの読み書き
処理ツール（マスク、再投影、リサンプリングなど）

https://rasterio.readthedocs.io/en/latest/

ラスターファイルを開く

import rasterio

src = rasterio.open("DEM_world.tif")

メタデータ:

src.count

src.width, src.height

(4320, 2160)

ラスターデータ = numpy配列

array = src.read()

標準のnumpy配列:

array

array([[[-4290, -4290, -4290, ..., -4290, -4290, -4290],
        [-4278, -4278, -4278, ..., -4278, -4278, -4278],
        [-4269, -4269, -4269, ..., -4269, -4269, -4269],
        ...,
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804]]], dtype=int16)

ラスターデータセットのプロット

rasterio.plot.show()を使用:

import rasterio.plot

rasterio.plot.show(src, cmap='terrain')

ベクターデータに基づく情報抽出

rasterstats: ベクター形状に基づき、地理空間ラスターデータの要約統計を算出 (https://github.com/perrygeo/python-rasterstats)

rasterstatsでラスタ値を抽出

ポイントベクター:

rasterstats.point_query(geometries, "path/to/raster", 
                        interpolation='nearest'|'bilinear')

ポリゴンベクター:

rasterstats.zonal_stats(geometries, "path/to/raster",
                        stats=['min', 'mean', 'max'])

rasterstatsでラスタ値を抽出

result = rasterstats.zonal_stats(countries.geometry, "DEM_gworld.tif", 
                                 stats=['mean'])

countries['mean_elevation'] = pd.DataFrame(result)

countries.sort_values('mean_elevation', ascending=False).head()

            name    continent                     geometry  mean_elevation
157   Tajikistan         Asia  POLYGON ((74.98 37.41, ...      3103.231105
85    Kyrgyzstan         Asia  POLYGON ((80.25 42.34, ...      2867.717142
24        Bhutan         Asia  POLYGON ((91.69 27.77, ...      2573.559846
119        Nepal         Asia  POLYGON ((81.11 30.18, ...      2408.907816
6     Antarctica   Antarctica  (POLYGON ((-59.57 -80.04...     2374.075028
..           ...          ...                          ...             ...

練習しましょう！

Pythonで扱う地理空間データ