QGISの通常クリギング

By K.Yoshimi

March 27, 2024

はじめに

本記事では、オープンソースの地理情報システムであるQGIS (LTR version 3.34.4)に付属したクリギング機能を使用して、地理情報に対して補間処理を行います。

地理情報データ

この記事では、クリギングに使用する地理情報データとして、GIS実習オープン教材「空間補間」¹ の実習用データ「interpolation」を使用します。

上記サイトからデータをダウンロードすると、富士山周辺の標高点データが入っています。

fuji.tif（GeoTIFF形式）：正解となる標高データ
elevation.shp（Shapefile形式）：fuji.tifからランダムに取り出した$5,087$点からなる標高データ

この記事では、このelevation.shpに含まれる標高データを測定データとして、均一な$30m$格子点（約$25$万点）上での標高をクリギングにより推測します。

QGISに、fuji.tifとelevation.shpをドラッグ&ドロップなどで読み込むと下図のようになります。

QGISのGUIの右下を見て分かる通り、座標系はEPSGコード:6676の平面直角座標系となっており、 fuji.tifのレイヤプロパティで確認すると、データの領域は

左下の点$(-6150.8, -95452.3)$、右上の点$(7548.3, -78765.2)$
x方向の幅：$13642.8m$
y方向の幅：$16405.2m$

となっていることを注意します。

fuji.tif

また、elevation.shpを読み込み、下図のようにランダムな標高点になっていることを確認できます。

ランダムな標高点

SAGAのクリギング機能

QGISには、モジュールライブラリとして「自動地球科学分析のためのシステム(System for Automated Geoscientific Analyses；SAGA GIS）」が同梱されており、それに含まれるクリギング機能を利用できます。

上記で紹介したGIS実習オープン教材「空間補間」では、SAGAのOrdinary Kriging（通常クリギング）を使用した例が計算されています。

手順を再掲しますと

プロセッシング > ツールボックス > SAGA

と進み

Spatial and Geostatistics - Kriging > Ordinary Kriging

をダブルクリックします。

SAGAのOridinary Kriging

Ordinary Krigingのパラメータに

Points: elevation [EPSG:6676]
Attribute: VALUE
Output extent: レイヤから計算 > elevation
Cellsize: $30.0$
Maximum Distance: $0.0$
Lag Distance Classes: $100$
Model: $a + b * x$

を設定して、実行ボタンクリックします。

❗ 測定点データが$5,087$点、推定点が約$25$万点あり、クリギング処理に時間がかかります。

実行すると、理論バリオグラムの推定とクリギングの推定が一度に行われますが、バージョンの違いも影響しているのか、標高の推定は上手くいっていないようです。

線形関数モデルによる推定結果

推定精度が悪い理由として、バリオグラムモデルに、デフォルトの「a + b * x」（線形モデル）を設定していますが、これが今回の標高データに対して適切かどうか不明である点が挙げられます。

経験バリオグラムの生成と理論バリオグラムの推定

以上で見たように、QGISに組み込まれたSAGAのOrdinary Krigingは、理論バリオグラムの推定から、クリギングによる補間までを一度に行え便利ですが、 elevation.shpに含まれるランダムな標高点データから経験バリオグラムを生成し、適切な理論バリオグラムの推定する部分がブラックボックス化され、クリギングの推定を改善するのが難しくなっています。

また、QGISの中には、理論バリオグラムを推定するための視覚的なツールが用意されておらず、経験バリオグラムと推定した理論バリオグラムの比較などを行うのが困難です。

そこで、SAGAのOrdinary Krigingを適切に行うためには、外部ツールを使用して、理論バリオグラムを推定しておくことが必要になってきます。

理論バリオグラムを推定するフリーのツールとして、R言語のライブラリ等が数多く存在しますが、ここでは、Pythonの地球統計バリオグラム推定ライブラリである scikit-gstatを使用します。

まず、elevation.shpから座標値と標高点データのCSVファイル：elevation.csvを作成しておきます。

elevation.zip

このelevation.csvを入力とする、scikit-gstatによる下記のPythonコードを使用して、経験バリオグラムとフィッティングした理論バリオグラムを確認します。

import argparse
import numpy as np
import skgstat as skg


def main():

    # ArgumentParserは、プログラムの引数をパースします。
    parser = argparse.ArgumentParser()
    parser.add_argument('--input-csv', default='elevation.csv')
    parser.add_argument('--model', default='exponential')
    parser.add_argument('--n-lags', type=int, default=20)
    parser.add_argument('--maxlag', type=float, default=12500)
    args = parser.parse_args()

    input_csv = args.input_csv
    model = args.model
    n_lags = args.n_lags
    maxlag = args.maxlag

    # numpyのgenfromtxt関数で指定したcsvファイルを読み込みます。
    data = np.genfromtxt(input_csv, dtype=None, names=True, deletechars='', delimiter=',', autostrip=True, encoding='utf_8_sig')

    x = data['x'].astype(np.float64)
    y = data['y'].astype(np.float64)
    elevs = data['VALUE'].astype(np.float64)
    coords = np.vstack((x,y)).T

    V = skg.Variogram(coordinates=coords, values=elevs)
    V.model = model
    V.n_lags = n_lags
    V.maxlag = maxlag

    print(V)

    fig = V.plot()
    fig.savefig('variogram.png')


if __name__ == '__main__':
    main()

手始めに、下記のパラメータを設定して、Pythonスクリプトを実行し、理論バリオグラムの推定を行います。

model: exponential
⬅️ 理論バリオグラムのモデル。ここでは指数関数モデル
n-lags: $25$
⬅️ 経験バリオグラムのラグ数
maxlag: $20000$ m
⬅️ 経験バリオグラムのラグの最大長。つまり、考慮する２測定点データ間の最大距離

実行：

> python skgstat_variogram.py --model exponential --n-lags 25 --maxlag 20000

実行結果：（理論バリオグラムの推定結果: maxlag=20000m）

exponential Variogram
---------------------
Estimator:         matheron
Effective Range:   6175.72
Sill:              33602.14
Nugget:            0.00

理論バリオグラムの推定: maxlag=20000m

経験バリオグラム(青点)は$12500m$付近まで上昇して、それ以降は相関を失いますので、パラメータを以下のように変更して、再度、pythonスクリプトを実行します。

model: exponential
n-lags: $20$
maxlag: $12500$ m

実行：

> python skgstat_variogram.py --model exponential --n-lags 20 --maxlag 12500

実行結果：（理論バリオグラムの推定結果）

exponential Variogram
---------------------
Estimator:         matheron
Effective Range:   11514.69
Sill:              44630.52
Nugget:            0.00

理論バリオグラムの推定: maxlag=12500m

今度は、ナゲット無しの指数関数モデルで経験バリオグラムがよくフィッティングされていることが確認できます。

❗ 実用的に適切な理論バリオグラムを推定するには、 n-lags（経験バリオグラムのラグ数）、maxlag（ラグの最大長）を調整したり、ナゲットの有無、複数のmodel(線形モデル、球関数モデル、指数関数モデル、ガウスモデル等）から、 AIC（赤池情報量基準）、BIC（ベイズ情報量基準）などの基準で適切なモデルを選択したりする必要があります。

なお、scikit-gstatにおける、指数関数モデルは次式で与えられます。

$$ \gamma = b + C_0 \times (1 - e^{-\frac{h}{a}}) $$

ここで

$b$（Nugget）: バリオグラムのナゲット
$C_0$（Sill）: バリオグラムのシル
$r$（Effective Range）: バリオグラムの実効範囲　($a=\frac{r}{3}$)

です。この式に、上記のPythonプログラムで推定した結果を代入すると

$$ \gamma = 44630.52*(1.0-\text{exp}(-x/3838.23)) $$

となります。

通常クリギングによる推定

以上で、理論バリオグラムの推定が完了しましたので、QGISのSAGAのOrdinary Krigingに、このナゲット無しの指数関数モデルを指定して、実行してみましょう。

Ordinary Krigingのパラメータとして

Points: elevation [EPSG:6676]
Attribute: VALUE
Output extent: レイヤから計算 > elevation
Cellsize: $30.0$
Maximum Distance: $12500.0$
Lag Distance Classes: $20$
Model: 44630.52*(1.0-exp(-x/3838.23))

を設定し、実行します。

❗ 測定点データが$5,087$点、推定点が約$25$万点あり、クリギング処理に時間がかかります。

推定結果は、下図のように大分改善されました。

指数関数モデルによる推定結果

また、推定誤差分散も出力され、下の右図のようになりました。

今回のように、理論バリオグラムのナゲットがゼロの場合は、観測点と一致する推定点において、推定値が観測値と一致するため、推定誤差分散も観測点付近では小さな値（黒い点）となります。したがって、ナゲット無しの場合は、観測点に近い位置ての推定値は滑らかでなくなる傾向があることに注意します。

指数関数モデルによる推定誤差分散

おわりに

本記事では、QGISにあるSAGAのクリギング機能を用いて、地理情報データの補間を試みました。

この機能は、経験バリオグラムの生成、理論バリオグラムの推定、クリギングによる推定を一度に行え、便利なのですが、実用的には、経験バリオグラムの生成、理論バリオグラムの推定は慎重に行う必要があり、この記事では外部ツールで対処しました。

このあたりをQGIS内で行えるとより便利なツールになるのではないかと思いました。

また、このSAGAのクリギング機能には、leave-one-out交差検証やk-分割交差検証の機能もありますので、適切な結果を得るために、これらをもっと活用できるかもしれません。

出典：GIS Open Educational Resources WG, CC BY-NC-SA 4.0 ↩︎

qgis
saga