NumPy配列の要素を、別の配列を参照して移動させたい

Question

前提・実現したいこと

PythonとGISを使用して地形発達シミュレーションを行おうとしています。
そこで、河川での侵食・堆積を表現したいです。

現在の状況

元となるデータは5mDEMを使用しています。形式はGeotiffのラスターファイルです。
まず、SAGA GISのChannel network and drainage basinsというものによって、流れる先を示すラスターデータを作成してあります。
下の表のように数字と方角が対応しています（ピクセルに０とあった場合には北へ流れる）。
8の場合は移動しません。

　|7|0|1|
|6| |2|
|5|4|3|

また、侵食量は式があり、求めることが出来ます。

考えたこと

このようにやればよいのではないかという予測
0. DEM及び流下先ラスターデータをnumpy配列に変換する

1. DEM配列の各要素ごとに侵食量を計算し、新たな配列（土砂配列）とする
2. 土砂配列を流下先配列に従って移動させる
3. 移動後の要素

わからないこと

やりたいことを画像で示すと、以下のようになります。

2×2の配列があった場合、移動後の土砂の配列が分かれば侵食・堆積後の標高は簡単な計算で求まります。

---

しかし、赤線で囲った部分の表現の仕方が見当もつきません。

調べた結果numpyの場合、cumsum関数を使えば隣の要素との計算が、スライシング機能によって要素の抽出ができると
分かりました。
また、roll関数によって位置をずらすことができると分かりました。
ただ、今回のように要素によってバラバラな動きをする方法は調べてもわかりませんでした。

教えていただけると幸いです。

補足情報（FW/ツールのバージョンなど）

ひとまず、移動距離は1マスとしたいです。

conda version          : 4.8.4
conda-build version  : 3.20.0
python version         : 3.8.5.final.0

Answer 1

移動後の土砂を求めるのに、
「下に移動する土砂」「上に移動する土砂」等、方向別に個別に求めるようにすると、numpy と相性が良くなりますよ。

t = np.zeros_lize(s) # 堆積する予定の行列・・・今は全て0
m = np.array(((4,4),(2,8))) #動き行列
s = np.array(((1,1),(1,1))) #浸食量行列

#以下を8方向分処理する・・下記は下向きの例
m_down = (m==4) #下に移動する場所は True でそれ以外はFalse の行列
s_down = s * m_down #浸食されて下に移動する量の行列
t[1:,:] += s_down[:-1,:] #下に移動した浸食量を堆積に追加

これを8方向について行えば赤枠の中の移動後の土砂が t に求まります。

Comments

[lawa]

回答ありがとうございます！
m==5のときも試してみたところ、無事上手くいきました。

スライシング機能によって移動を実現しているんですね。
スライシング機能をそのように使うとは思いつきもしませんでした。

Answer 2

流下先配列の特定の座標に着目したとき、その座標の値が「4」ならば、一つ下の座標に侵食量の値がプラスされると、いうことかと思います。
逆に言えば、ある特定座標の上の座標の値が「4」ならば、その座標に上の座標の侵食量の値がプラスされるということになります。

したがって配列を畳み込み演算のように3*3のブロックごとに処理していけば良いです。

from skimage.util import view_as_windows

# 移動マップ
d = np.array([[7, 0, 1],
              [6, 9, 2],
              [5, 4, 3]])[::-1, ::-1]

m = np.array([[4, 4], [2, 8]])  # 流下先配列
s = np.array([[1, 1], [1, 1]])  # 侵食量配列

# 移動後の土砂配列の計算
m_block = view_as_windows(np.pad(m, 1, constant_values=8), (3, 3))
# [[8 8 8 8]
#  [8 4 4 8]       [[888]  [[888]   [[844]   [[448]
#  [8 2 8 8]     [  [844]   [448]    [828]    [288]  ]
#  [8 8 8 8]] ->    [828]], [288]],  [888]],  [888]] 
s_block = view_as_windows(np.pad(s, 1, constant_values=0), (3, 3))
c = np.sum(((m_block == d) * s_block), (2, 3))

c
# array([[0, 0],
#        [1, 2]])

np.random.seed(0)
m = np.random.randint(0, 9, (5, 5))  # 流下先配列
# [[5 0 3 3 7]
#  [3 5 2 4 7]
#  [6 8 8 1 6]
#  [7 7 8 1 5]
#  [8 4 3 0 3]]
s = np.random.randint(1, 3, (5, 5))  # 侵食量配列
# [[2 1 2 2 1]
#  [1 2 1 2 2]
#  [2 2 2 1 2]
#  [1 2 2 2 2]
#  [1 2 1 1 2]]

# 移動後の土砂配列の計算
m_block = view_as_windows(np.pad(m, 1, constant_values=8), (3, 3))
s_block = view_as_windows(np.pad(s, 1, constant_values=0), (3, 3))
c = np.sum(((m_block == d) * s_block), (2, 3))
# [[0 0 0 2 0]
#  [0 0 0 3 3]
#  [4 1 0 4 2]
#  [0 0 0 1 0]
#  [0 0 0 2 0]]

---

（以下は以前の解答）

全体的な処理フローは専門外のためよくわかってないですが（すみません）、
とりあえず「赤枠内でやりたいこと」というのはこういうことでしょうか……？

from skimage.util import view_as_windows

d = np.array([[7, 0, 1],
              [6, 9, 2],
              [5, 4, 3]])
m = view_as_windows(d, (2, 2))

v = np.array([[4, 4],
              [2, 8]])

np.sum(m == v[::-1, ::-1], (2, 3))
# array([[0, 0],
#        [1, 2]])

Comments

[lawa]

素早い回答ありがとうございます！
あまりよくわかっていないので、確認しつつ質問させていただきたいのですが、
・view_as_windowsで区切ったのは、元の2×2の配列の大きさに合わせるためですか
・私が画像の例で出した侵食量の配列はどこに反映されているのですか
・v[::-1, ::-1]は反転のような方法を取っているということですか
・本番では何千×何千の配列になるのですが、そのときも2×2に区切って処理をした方がよいのでしょうか

質問ばかりで申し訳ありませんが、よろしくお願いいたします。

[lawa]

返信が遅くなってしまい、申し訳ありません。

「どこへ移動させるか」ではなく、「どこから移動してくるか」という逆の発想なんですね。
詳しく教えていただき、ありがとうございます！