多次元配列からのフィーチャのインデックスの抽出

Question

3次元空間における雲水濃度の値を表す配列を持っています。雲水濃度がある閾値を上回る場所では、私は雲があると言います（以下の断面図を参照）。ドメインの大半は乾燥していますが、ドメインの大部分には約400メートルの基盤を持つ成層雲があります。

私は何をしたいの（x、y、z）は、クラウドベースおよびクラウドトップの位置座標を抽出しています。次に、風速の垂直成分を表す別の3次元配列上でこれらの座標を使用して、クラウドベースで上昇気流を取得します。

私が今やっていることはうまくいくが、遅いです。私はそれをスピードアップするためにNumPyを利用する方法がなければならないと感じています。

これは私が今やっているものです：

# 3d array representing cloud water at a particular timestep t 
qc = QC(t) 

# get the coordinates where there is cloud 
cloud_coords = argwhere(qc > qc_thresh) 

# Arrays to hold the z values of cloud base (cb) and cloud top (ct) 
zcb = zeros((nx,ny)) 
zct = zeros((nx,ny)) 

# Since each coordinate (x,y) will in general have multiple z values 
# for cloud I have to loop over all (x,y) and 
# pull out max and min height for each point (x,y) 
for x in range(nx): 
    # Pull out all the coordinates with a given x value 
    xslice = cloud_coords[ where(cloud_coords[:,0] == x) ] 

    for y in range(ny):  
     # for the given x value select a particular y value 
     column = xslice[ where(xslice[:,1] == y) ] 

     try: 
      zcb[x,y] = min(column[:,2]) 
      zct[x,y] = max(column[:,2]) 
     except: 
      # Because there may not be any cloud at all 
      # (a "hole") we fill the array with an average value 
      zcb[x,y] = mean(zcb[zcb.nonzero()]) 
      zct[x,y] = mean(zct[zct.nonzero()]) 


# Because I intend to use these as indices I need them to be ints 
zcb = array(zcb, dtype='int') 
zct = array(zct, dtype='int') 

出力は、クラウドベース（トップ）のZ座標を含む2次元配列である

私は、これらを使用します別の配列のインデックスを使ってクラウドベースの風速などの変数を取得する：

wind = W(t) 
j,i = meshgrid(arange(ny),arange(nx)) 
wind_base = wind[i,j,zcb] 

私はシミュレーションの多くのタイムステップでこれを行い、最も遅い部分はすべての（x、y）座標に対するPythonループです。 NumPyを使用してこれらの値をより速く抽出することについての助けがあれば、大いに感謝します。

Answer 1

あなたの問題でnumpyをうまく使うことができるという疑いが正しいです。実際には、明示的にあなたがのほとんどを行うことができます。3.

Pythonで複雑なオブジェクトがあるintのdtypeと新しい終わりnp.array()を使用してアレイ、および1つのインスタンス生成のためにあなたがやっている複数の非効率性があり、いくつかのベクトル化されたnumpyの行で動作します。アイデアは、雲が出現するインデックス、または雲が終わるインデックス（高さ軸に沿った）を見つけるだけで十分です。 numpy.argmaxを使用してベクトル化された方法で行うことができます。それは本当に効率的なソリューションの心臓部です：

import numpy as np 
import matplotlib.pyplot as plt 

# generate dummy data 
qc_thresh = 0.6 
nx,ny,nz = 400,400,100 
qc = np.zeros((nx,ny,nz)) 
# insert random cloud layer 
qc[...,50:80] = np.random.rand(nx,ny,30) 
# insert holes in clouds for completeness 
qc[np.random.randint(nx,size=2*nx),np.random.randint(ny,size=2*nx),:] = 0 

def compute_cloud_boundaries(): 
    cloud_arr = qc > qc_thresh 

    # find boundaries by making use of np.argmax returning first maximum 
    zcb = np.argmax(cloud_arr,axis=-1) 
    zct = nz - 1 - np.argmax(cloud_arr[...,::-1],axis=-1) 

    # logical (nx,ny)-shaped array where there's a cloud 
    cloud_inds = (zcb | (zct!=nz-1)).astype(bool) 
    # this is short for `(zcb==0) | (zct!=nz-1)` 

    # fill the rest with the mean 
    zcb[np.logical_not(cloud_inds)] = zcb[cloud_inds].mean() 
    zct[np.logical_not(cloud_inds)] = zct[cloud_inds].mean() 

    return zcb,zct 

私はあなたのアプローチに対して、（対応する小さな例で完了）上記を確認し、それが正確に同じ結果を与えます。私が言ったように、考え方は、cloud_arr = qc > qc_threshは、湿度が雲に適しているかどうかを示す論理的な配列です。次に、最後（高さ）の軸に沿ってこの（本質的にバイナリ！）行列の最大値を調べます。 np.argmaxを呼び出すと、各平面2Dインデックスの最初（最下部）の高さの値がわかります。雲の頂点を得るためには、配列 を逆にして、反対側から同じことをする必要があります（結果のインデックスを変換して処理します）。配列を逆にすると、コピーではなくビューが作成されるため、効率的です。最後に、雲がない点を修正します。より良い制約の代わりに、argmaxによって返される最も高いインデックスがエッジポイントに対応する箇所を確認します。現実の気象データを考慮すると、最下部と最上部の測定値がでないことが確認できます。は雲に対応しているため、これは安全な基準である必要があります。

ここショーのためのダミーデータの断面は次のとおり

In [24]: %timeit compute_cloud_boundaries() 
10 loops, best of 3: 29.1 ms per loop 

In [25]: %timeit orig() # original loopy version from the question 
1 loop, best of 3: 9.37 s per loop 

300以上であると思われる：上記400x400x100場合の

非代表タイミング倍のスピードアップ。もちろん、実際のユースケースはこのアプローチの適切なテストになりますが、うまくいくはずです。

---

インデックス作成のステップでは、インデックス用にオープングリッドを使用し、配列ブロードキャストを使用してメモリを確保することができます。少しもこの段階スピードアップするかもしれない追加の(nx,ny)字型の配列を割り当てることがない：

wind = W(t) 
i,j = np.ogrid[:nx,:ny] 
wind_base = wind[i,j,zcb] 

あなたが見ることができるように、np.ogridはmeshgridと同等であるものに一緒に放送形状(nx,1)と(1,ny)のオープングリッドを作成しますコール。