np.gradientを用いて地形の傾斜を計算したい

Question

np.gradientを用いて地形の傾斜を計算したい

現在，地形の傾斜のデータをPythonを用いて取得しようとしています．
計算の高速化のため，for文を使用せず，np.gradientを用いて，
一括で傾斜する計算する方法を考えているのですが，
それがわからないため質問させていただきました．

実現したいこと

以下の傾斜を求める手法をnp.gradientで実現したいと考えています．
これはppにおける傾斜を該当グリッドの周囲8方向から取得するものです．

この方法を実現するために，以下のようなfor文を入れたコードを作成しましたが，
for文を用いずに傾斜を求める手法をご教授頂ければと思います．

該当のソースコード

#モジュールのインポート
import numpy as np
import math

#標高データの準備
current_dem = np.array([[1,2,3,4], [5,7,9,10], [8,15,12,9], [8,5,3,1]])

#傾斜の格納要配列の準備
value_slope = np.zeros((4, 4))
value_slope[:,:] = np.nan

#標高（グリッドデータ)のx,y方向の間隔
DX = 1
DY = 1

#傾斜の計算
for i in range(4):
    for j in range(4):
        if i != 0 and j != 0 and i != 3 and j != 3:

            pp = current_dem[i][j]

            d1 = current_dem[i+1][j-1]

            d2 = current_dem[i+1][j]

            d3 = current_dem[i+1][j+1]

            d4 = current_dem[i][j-1]

            d6 = current_dem[i][j+1]

            d7 = current_dem[i-1][j-1]

            d8 = current_dem[i-1][j]

            d9 = current_dem[i-1][j+1]

            xx = (d1+d4+d7-(d3+d6+d9))/(6*DX)
            yy = (d7+d8+d9-(d1+d2+d3))/(6*DY)

            value_slope[i, j] = math.sqrt(xx*xx + yy*yy)

print(value_slope)
#[[       nan        nan        nan        nan]
 #[       nan 5.11262055 4.50308536        nan]
 #[       nan 0.97182532 3.06412939        nan]
 #[       nan        nan        nan        nan]]

Comments

[jbpb0]

参考
https://algorithm.joho.info/image-processing/prewitt-filter/
https://algorithm.joho.info/programming/python/opencv-prewitt-filter-py/

Answer 1

今回のような計算で、forループをなくして高速化をするのが目的なら畳み込みを使うのがいいと思います。numpyには二次元の畳み込みの関数はないので、scipyのconvolve2dを使います。
(scipyなしでもできますが、scipyを使ったほうが速度面でも有利じゃないかと)

import numpy as np
from scipy.signal import convolve2d

#標高データの準備
current_dem = np.array([[1,2,3,4], [5,7,9,10], [8,15,12,9], [8,5,3,1]])

#標高（グリッドデータ)のx,y方向の間隔
DX = 1
DY = 1

grad_x = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]]) / (6 * DX)
grad_y = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]]) / (6 * DY)

xx = convolve2d(current_dem, grad_x, mode='same', boundary='fill', fillvalue=np.nan)
yy = convolve2d(current_dem, grad_y, mode='same', boundary='fill', fillvalue=np.nan)

value_slope = np.sqrt(xx * xx + yy * yy)
print(value_slope)

numpy.gradient を使ってやるなら、こんな感じではないでしょうか。
(最終的な計算結果には影響しないですが、上の畳み込みの場合とxx, yyの符号が違っているので注意ください)

ey, ex = np.gradient(current_dem, DY, DX)

xx = (ex[:-2, 1:-1] + ex[1:-1, 1:-1] + ex[2:, 1:-1]) / 3
yy = (ey[1:-1, :-2] + ey[1:-1, 1:-1] + ey[1:-1, 2:]) / 3

# 上の部分は sliding_window_view を使って下記のようにも書けます
#from numpy.lib.stride_tricks import sliding_window_view
#xx = np.mean(sliding_window_view(ex[:, 1:-1], 3, axis=0), axis=-1)
#yy = np.mean(sliding_window_view(ey[1:-1, :], 3, axis=1), axis=-1)

value_slope = np.full(current_dem.shape, np.nan)
value_slope[1:-1, 1:-1] = np.sqrt(xx * xx + yy * yy)
print(value_slope)

Comments

[gantakun]

回答ありがとうございます．参考にさせていただきます．scipyは使用した経験がなかったのですが，
これを機に学んでみたいと思います．

Answer 2

実質的に for ループと等価なので参考までにどうぞ。

#傾斜の計算
def calc_slope(i, j):
    xx = np.sum(current_dem[i-1,j-1:j+2] - current_dem[i+1,j-1:j+2])/(6*DX)
    yy = np.sum(current_dem[i-1:i+2,j-1] - current_dem[i-1:i+2,j+1])/(6*DX)
    value_slope[i, j] = np.sqrt(xx**2 + yy**2)

r, c = current_dem.shape
prd = np.dstack(np.meshgrid(range(1,r-1), range(1,c-1))).reshape(-1, 2).T
np.vectorize(calc_slope)(*prd)

Comments

[gantakun]

回答ありがとうございます．参考にさせていただきます．