概要
- GRASS GIS r.mapcalc.simpleでは周りのセルの値を相対のindexで取得できる
- ただ、GRASSはWindowsのQGISのインストーラーがMSIになってから使えていない。
- 2024/8/7追記：QGIS 3.32あたりからGRASSが使えるようになっています。 Fix GRASS7/SAGA UnicodeDecodeError in QGIS of Windows-Japanese or others #49226
- GDALラスタ計算機は機能豊富なnumpy関数が使えるからなんとかならんか という
GDALラスタ計算機
- A,B,C ... の文字に入力ラスタを割り当て、「A*2」のように式を入れる。
- A[0,0]でラスタの起点セルの値を取得できるが、絶対のindexなので、当該セルの周りの値は取得できない。
- 加減乗除+-*/は普通に書ける。
- 等号、不等号=<>は使えず「equal(A, 0)」のような関数で書く。
- ラスタ単位での処理なので、個別のセルに対するforループのような処理はできない。
- 論理関数 https://numpy.org/doc/stable/reference/routines.logic.html
- 数学関数 https://numpy.org/doc/stable/reference/routines.math.html
いろいろ試す
- indexとなる連続値のラスタを新たに作ればいいのでは、、と思って作ってみた。
- が、「A[B,C]」のような式はB,Cが当該の値ではなくラスタなので、使えなかった。
ラスタサイズの取得 shape
- 形状のタプルを返す
- https://numpy.org/doc/stable/reference/generated/numpy.shape.html#numpy.shape
- shape(A)[0] ->Y軸（高さ）
- shape(A)[1] ->X軸（幅）
同形状で値が0のラスタ zeros
- ゼロで埋めた配列を返す
- zeros(shape(A)) -> Aと同形状で0で埋めたラスタを返す
- 返す配列の形状が同じであれば、新たに作ってもよいらしい
起点から1セルごとに値が１増加するラスタ
- repeat([arange(shape(A)[1])],shape(A)[0] ,axis=0) ->X軸方向に連続値
- repeat([arange(shape(A)[0])],shape(A)[1] ,axis=0).T ->Y軸方向に連続値
発想の転換、周りの値を移動させる
roll 関数
- 値をシフトする関数で、左に1セルシフトすれば、1つ右側の値を取得できる。
- ただし、左1セルシフトすると、右端の値は左端の値になる。
- roll(A, (0, -1)) →左に1セルシフト（右側の値を取得）
- (roll(A, (0, -1)) + roll(A, (0, 1)))/2 → 左右のセル値の平均
- https://numpy.org/doc/stable/reference/generated/numpy.roll.html?highlight=roll#numpy.roll
ラスタの外周セルの削除をラスタ計算機でできないか
- 周囲１セルを計算に使用した場合、ラスタの外周の１セル分は意味のない値。
- 先ほどの軸のインデックスラスタを使って、外周の値の削除できるが、4辺あるので、煩雑
- たとえば、X軸の左端が0のラスタBを使って、Aの左端を0に
- not_equal(B,0) * A
外周を削除するのは、他のプロセッシングで切り抜いた方速い
- レイヤ範囲の抽出
- バッファ(マイナス１セル分)
- マスクレイヤで切り抜く（マスクレイヤの領域に切り抜き範囲を一致させる のチェックを外す）
おわりに
- QGISのGDALラスタ計算機で周囲の値を使って計算することはできないことはないが、ややトリッキーな気がする。
- もっと簡単な方法があったらコメントください。