3

ポリゴンの値を細かい規則的なグリッドにビン化したい。たとえば、次の座標があります。

data = 2.353
data_lats = np.array([57.81000137,  58.15999985,  58.13000107,  57.77999878])
data_lons = np.array([148.67999268,  148.69999695,  148.47999573,  148.92999268])

私の通常のグリッドは次のようになります。

delta = 0.25
grid_lons = np.arange(-180, 180, delta)
grid_lats = np.arange(90, -90, -delta)
llx, lly = np.meshgrid( grid_lons, grid_lats )
rows = lly.shape[0]
cols = llx.shape[1]
grid = np.zeros((rows,cols))

これで、ポリゴンの中心に対応するグリッド ピクセルを非常に簡単に見つけることができます。

centerx, centery = np.mean(data_lons), np.mean(data_lats)
row = int(np.floor( centery/delta ) + (grid.shape[0]/2))
col = int(np.floor( centerx/delta ) + (grid.shape[1]/2))
grid[row,col] = data

ただし、まだポリゴンと交差しているグリッド ピクセルがいくつか存在する可能性があります。したがって、ポリゴン内に一連の座標 (data_lons、data_lats) を生成し、以前のように対応するグリッド ピクセルを見つけたいと思います。座標をランダムまたは体系的に生成することを提案しますか? 私は失敗しましたが、まだ努力しています。

注: 1 つのデータ セットには約 80000 個のポリゴンが含まれているため、非常に高速 (数秒) である必要があります。それが、オーバーラップの領域を考慮していないため、このアプローチを選択した理由でもあります...（私の以前の質問のようにデータビニング：非常に遅い通常のメッシュへの不規則ポリゴン）

4

2 に答える 2

1

角のピクセル間の座標を計算するだけで、簡単で汚いソリューションに取り組みました。見てみましょう：

dlats = np.zeros((data_lats.shape[0],4))+np.nan
dlons = np.zeros((data_lons.shape[0],4))+np.nan
idx = [0,1,3,2,0] #rearrange the corner pixels

for cc in range(4):
    dlats[:,cc] = np.mean((data_lats[:,idx[cc]],data_lats[:,idx[cc+1]]), axis=0)
    dlons[:,cc] = np.mean((data_lons[:,idx[cc]],data_lons[:,idx[cc+1]]), axis=0)

data_lats = np.column_stack(( data_lats, dlats ))
data_lons = np.column_stack(( data_lons, dlons ))

したがって、赤い点は元のコーナーを表し、青い点はそれらの間の中間ピクセルを表します。

ここに画像の説明を入力

これをもう一度実行して、中央のピクセル (geo[:,[4,9]]) を含めることができます。

dlats = np.zeros((data.shape[0],8))
dlons = np.zeros((data.shape[0],8))

for cc in range(8):
    dlats[:,cc] = np.mean((data_lats[:,cc], geo[:,4]), axis=0)
    dlons[:,cc] = np.mean((data_lons[:,cc], geo[:,9]), axis=0)

data_lats = np.column_stack(( data_lats, dlats, geo[:,4] ))
data_lons = np.column_stack(( data_lons, dlons, geo[:,9] ))

ここに画像の説明を入力

これは非常にうまく機能し、次のように各ポイントを対応するグリッド ピクセルに直接割り当てることができます。

row = np.floor( data_lats/delta ) + (llx.shape[0]/2)
col = np.floor( data_lons/delta ) + (llx.shape[1]/2)

ただし、最終的なビニングには約 7 秒かかります!!! このコードを高速化するにはどうすればよいですか:

for ii in np.arange(len(data)):
    for cc in np.arange(data_lats.shape[1]):
        final_grid[row[ii,cc],col[ii,cc]] += data[ii]
        final_grid_counts[row[ii,cc],col[ii,cc]] += 1
于 2013-03-06T09:43:30.947 に答える
1

次のアプローチをテストして、十分に高速かどうかを確認する必要があります。まず、すべての緯度と経度を次のように変更して、グリッドに (おそらく分数の) インデックスを作成する必要があります。

idx_lats = (data_lats - lat_grid_start) / lat_grid step
idx_lons = (data_lons - lon_grid_start) / lon_grid step

次に、ポリゴンを三角形に分割します。任意の凸多角形の場合、多角形の中心をすべての三角形の 1 つの頂点として取得し、次に多角形の頂点を連続したペアとして取得できます。しかし、多角形がすべて四角形の場合は、最初の頂点に 0、1、2、2 番目の頂点に 0、2、3 を使用して、それらを 2 つの三角形のみに分割する方が速くなります。

特定の点が三角形の内側にあるかどうかを知るために、ここで説明されている重心座標アプローチを使用します。この最初の関数は、一連の点が三角形の内側にあるかどうかをチェックします。

def check_in_triangle(x, y, x_tri, y_tri) :
    A = np.vstack((x_tri[0], y_tri[0]))
    lhs = np.vstack((x_tri[1:], y_tri[1:])) - A
    rhs = np.vstack((x, y)) - A
    uv = np.linalg.solve(lhs, rhs)
    # Equivalent to (uv[0] >= 0) & (uv[1] >= 0) & (uv[0] + uv[1] <= 1)
    return np.logical_and(uv >= 0, axis=0) & (np.sum(uv, axis=0) <= 1)

頂点によって三角形が与えられると、三角形のバウンディング ボックス内の格子点で上記の関数を実行することにより、その内部の格子点を取得できます。

def lattice_points_in_triangle(x_tri, y_tri) :
    x_grid = np.arange(np.ceil(np.min(x_tri)), np.floor(np.max(x_tri)) + 1)
    y_grid = np.arange(np.ceil(np.min(y_tri)), np.floor(np.max(y_tri)) + 1)
    x, y = np.meshgrid(x_grid, y_grid)
    x, y = x.reshape(-1), y.reshape(-1)
    idx = check_in_triangle(x, y, x_tri, y_tri)
    return x[idx], y[idx]

四角形の場合は、この最後の関数を 2 回呼び出すだけです。

def lattice_points_in_quadrilateral(x_quad, y_quad) :
    return map(np.concatenate,
               zip(lattice_points_in_triangle(x_quad[:3], y_quad[:3]),
                   lattice_points_in_triangle(x_quad[[0, 2, 3]],
                                              y_quad[[0, 2, 3]])))

サンプル データに対してこのコードを実行すると、2 つの空の配列が返されます。これは、四角形の点の順序が驚くべきものであるためです。インデックス 0 と 1 は一方の対角線を定義し、2 と 3 は他方を定義します。上記の関数は、頂点がポリゴンの周りに並べられることを期待していました。実際にこれとは別の方法で処理を行っている場合は、2 番目の呼び出しをlattice_points_in_triangleinsideに変更しlattice_points_in_quadrilateralて、使用されるインデックス[0, 1, 3]が[0, 2, 3].

そして今、その変更で：

>>> idx_lats = (data_lats - (-180) ) / 0.25
>>> idx_lons = (data_lons - (-90) ) / 0.25
>>> lattice_points_in_quadrilateral(idx_lats, idx_lons)
[array([952]), array([955])]

グリッドの解像度を 0.1 に変更した場合:

>>> idx_lats = (data_lats - (-180) ) / 0.1
>>> idx_lons = (data_lons - (-90) ) / 0.1
>>> lattice_points_in_quadrilateral(idx_lats, idx_lons)
[array([2381, 2380, 2381, 2379, 2380, 2381, 2378, 2379, 2378]),
 array([2385, 2386, 2386, 2387, 2387, 2387, 2388, 2388, 2389])]

タイミングに関しては、このアプローチは、私のシステムでは、ニーズに対して約 10 倍遅くなります。

In [8]: %timeit lattice_points_in_quadrilateral(idx_lats, idx_lons)
1000 loops, best of 3: 269 us per loop

つまり、20秒以上見ています。80,000 ポリゴンを処理します。

于 2013-03-05T22:34:21.060 に答える