for文を減らし，ユークリッド距離の計算速度を向上させたい

前提・実現したいこと

ユークリッド距離の計算を高速化したいと考えています．
しかし，for文を多用する方法でしか実現できていません．
この結果，コードが煩雑になり，計算速度も出ません．
そこで，for文を減らし，計算の高速化が実現できる方法をご教授願いたいです．

###詳細
やりたいことは，陸地の各地点から，海までの最短距離を求めることです．
最初に，標高データを作成しています．その後の処理では，最初の2回のfor文で，
標高データから陸地を探し，その陸地の座標を獲得します．次の2回のfor文で，その陸地から海に該当する全地点の距離を
求めます．求めた距離を距離格納用配列に入れます．その後，距離格納用配列の最小値を
「海までの最短距離」として得ています．

このような規模の標高データですと，計算時間はあまり問題になりませんが，
1000*1000の標高データだと，計算に時間がかかってしまいます．

発生している問題・エラーメッセージ

以下に記載したコードでは，エラーは発生しません．
しかし，扱うデータが大きくなった場合，計算時間がかかるので，
その高速化を実現したいと考えています．

該当のソースコード

Python
1#ライブラリ・モジュールの読み込み
2import numpy as np
3
4#標高データ作成
5dem = np.array([(-10, -5, 0, 5, 10), 
6                (-5, 0, 10, 12, 15), 
7                (-5, -5, 10, 20, 20), 
8                (0, -5, 0, 10, 25), 
9                (10, 15, 20, 25, 30)])
10
11#距離計算用の配列準備
12d1 = np.zeros((5, 5))
13d1[:,:] = np.nan
14
15#距離格納用の配列準備
16d2 = np.zeros((5, 5))
17
18#各地点における，海までの最短距離を算出
19dem_interval = 15　#標高データが格納されている間隔
20i = 5
21j = 5
22
23for k in range(i):
24    
25    for l in range(j):
26        
27        if DEM[k, l] >= 0:
28            
29            x = np.array([k, l]) #該当DEMの座標を格納
30        
31            for m in range(i):
32            
33                for n in range(j):
34                    
35                    if DEM[m, n] < 0:
36                        
37                        y = np.array([m, n]) #海のDEMの座標を格納
38                        d1[m, n] = np.linalg.norm(x - y) #該当DEMと海DEMが異なる場合は，DEM間の距離計測
39            
40            d2[k, l] = np.nanmin(d1)
41            
42        else:
43            
44            d2[k, l] = 0
45
46d2 = dem_interval*d2     
47print(d2)
48

試したこと

Numpy配列では，for文をベクトル化すれば計算速度の高速化が見込めるとネットには書いてありました．
しかし，自分が未熟なため，その情報を自分のスクリプトに落とし込めませんでした．

行動規範の内容に同意します

回答1件

ベストアンサー

scipyを使ってみてはどうでしょうか。
ndimage.distance_transform_edt を使うと各地点から最寄りのバックグラウンド(値が0)までのユークリッド距離を計算することができます。
最初のdemデータの標高が0以上の部分を1、それ以外を0にしたデータを作成し、この関数で処理するだけです。

python
1import numpy as np
2from scipy import ndimage
3
4def solve(dem, dem_interval=15):
5    field = np.where(dem >= 0, 1, 0)
6    res = ndimage.distance_transform_edt(field)
7    return res * dem_interval
8
9dem = np.array([(-10, -5,  0,  5, 10),
10                (-5,   0, 10, 12, 15),
11                (-5,  -5, 10, 20, 20),
12                ( 0,  -5,  0, 10, 25),
13                (10,  15, 20, 25, 30)])
14
15res = solve(dem)
16print(res)
17
18#[[ 0.          0.         15.         30.         45.        ]
19# [ 0.         15.         21.21320344 33.54101966 47.4341649 ]
20# [ 0.          0.         15.         30.         45.        ]
21# [15.          0.         15.         30.         45.        ]
22# [21.21320344 15.         21.21320344 33.54101966 47.4341649 ]]