Pythonで3次元配列を補間する

Question

私の質問は、ここに表示されているコード応答を拡張しています。Interpolating a 3d array in Python. How to avoid for loops?。関連元の溶液のコードは以下の通りです：

import numpy as np 
from scipy.interpolate import interp1d 
array = np.random.randint(0, 9, size=(100, 100, 100)) 
x = np.linspace(0, 100, 100) 
x_new = np.linspace(0, 100, 1000) 
new_array = interp1d(x, array, axis=0)(x_new) 
new_array.shape # -> (1000, 100, 100) 

x_newが一定の1次元配列であるが、私のx_new が一定の1次元配列ではなく、インデックスに依存してどのような場合時には、上記のアプローチは素晴らしい作品他の3次元配列の緯度/経度次元の値。私のx_newは、サイズが355x195x192（時間x緯x長）で、今は緯度と経度の寸法を二重にしています。 x_newは緯度/経度のペアごとに異なるので、以下のようにループを回避するにはどうすればよいですか？ループ処理には数時間かかります...

x_new=(np.argsort(np.argsort(modell, 0), 0).astype(float) + 1)/np.size(modell, 0) 
## x_new is shape 355x195x192 
## pobs is shape 355x1 
## prism_aligned_tmax_sorted is shape 355x195x192 
interp_func = interpolate.interp1d(pobs, prism_aligned_tmax_sorted,axis=0) 
tmaxmod = np.empty((355, 195, 192,)) 
tmaxmod[:] = np.NAN          
for latt in range(0, 195): 
    for lonn in range(0, 192): 
     temp = interp_func(x_new[:,latt,lonn]) 
     tmaxmod[:,latt,lonn] = temp[:,latt,lonn] 

ありがとうございました！

Answer 1

私はあなたがそれらのループを取り除く方法を知っていますが、あなたはそれを気に入らないでしょう。

問題がinterp1dのこの使用は、あなたが本質的に1Dドメイン、各スカラーxためあなたは2次元アレイ状Fを有する即ちF(x)機能上補間行列値関数を与えることです。あなたがしようとしている補間はむしろこれです。それぞれの(lat,lon)ペアのための個々の補間を作成することです。これはF_(lat,lon)(x)の行に沿っています。

この問題が発生する理由は、ユースケースでは、クエリポイントごとに行列値のF(x)を計算していますが、1つの要素を除くすべての行列要素を破棄することになりますこのペアに対応するクエリポイントの場合は[lat,lon]）。したがって、無関係な関数値をすべて計算する不必要な計算をしています。問題は、より効率的な方法があるかどうかわからないということです。

あなたのユースケースは、背中の裏側に適切なメモリで固定することができます。あなたのループが何時間も実行されているという事実は、これが実際にあなたのユースケースでは可能ではないことを示唆していますが、とにかく私はこの解決方法を示します。アイデアは、あなたの3d配列を2次元に変え、この形の補間を行い、補間された結果の有効2d部分空間に沿って対角要素を取ることです。クエリポイントごとに無関係の行列要素をすべて計算しますが、配列をループする代わりに、単一の索引付け手順で関連する行列要素を抽出することができます。これのコストは、はるかに大きな補助配列を作成することであり、使用可能なRAMに収まらない可能性が最も高いでしょう。

import numpy as np 
from scipy.interpolate import interp1d 
arr = np.random.randint(0, 9, size=(3, 4, 5)) 
x = np.linspace(0, 10, 3) 
x_new = np.random.rand(6,4,5)*10 

## x is shape 3 
## arr is shape 3x4x5 
## x_new is shape 6x4x5 

# original, loopy approach 
interp_func = interp1d(x, arr, axis=0) 
res = np.empty((6, 4, 5)) 
for lat in range(res.shape[1]): 
    for lon in range(res.shape[2]): 
     temp = interp_func(x_new[:,lat,lon]) # shape (6,4,5) each iteration 
     res[:,lat,lon] = temp[:,lat,lon] 

# new, vectorized approach 
arr2 = arr.reshape(arr.shape[0],-1) # shape (3,20) 
interp_func2 = interp1d(x,arr2,axis=0) 
x_new2 = x_new.reshape(x_new.shape[0],-1) # shape (6,20) 
temp = interp_func2(x_new2) # shape (6,20,20): 20 larger than original! 
s = x_new2.shape[1] # 20, used for fancy indexing ranges 
res2 = temp[:,range(s),range(s)].reshape(res.shape) # shape (6,20) -> (6,4,5) 

結果resとres2配列がまったく同じなので、アプローチはおそらく動作します：

とにかく、ここ1であなたの現在のアプローチを比較し、アクションでトリックです。しかし、私が言ったように、形状のクエリー配列(nx,nlat,nlon)に対しては、通常は膨大なメモリが必要な形状(nx,nlat*nlon,nlat*nlon)の補助配列が必要です。

---

私はちょうど1つによって、あなたの1D補間を実行していると考えることができる唯一の厳格な選択肢：ダブルループでnlat*nlon補間を定義します。これは補間器を作成するオーバーヘッドが大きくなりますが、一方で、不要な作業を補間して配列値を計算することはありません。

最後に、多変量補間を使用することを検討してください（私はinterpolate.interpndまたはinterpolate.griddataと考えています）。あなたの関数が緯度と経度の関数として滑らかであると仮定すると、より高い次元で完全なデータセットを補間するのが理にかなっているかもしれません。このようにすれば、一度インターポレータを作成し、必要なポイントで正確にクエリを実行する必要があります。

---

、あなたの現在の実装にこだわってしまう場合は、おそらく大幅に最後の位置にあなたの補間軸を移動することで、パフォーマンスを向上させることができます。このように、ベクトル化された各演算は、連続したメモリブロック（デフォルトのCメモリ順序を前提とします）で動作し、これは "1d配列のコレクション"の考え方によく合います。だから、あなたは元の順序を復元する必要がある場合は、あなたが最後にres.transpose(2,0,1)バックオーダーを移調することができます

arr = arr.transpose(1,2,0) # shape (4,5,3) 
interp_func = interp1d(x, arr, axis=-1) 
... 
for lat ...: 
    for lon ...: 
     res[lat,lon,:] = temp[lat,lon,:] # shape (4,5,6) 

の線に沿って何かをする必要があります。