13

地球の小さな領域を表す何千ものポリゴンがテーブル形式で保存されています (4 つのコーナー座標が与えられます)。さらに、各ポリゴンにはデータ値があります。ファイルは、たとえば次のようになります。

lat1,  lat2,  lat3,  lat4,  lon1,   lon2,   lon3,   lon4,   data
57.27, 57.72, 57.68, 58.1,  151.58, 152.06, 150.27, 150.72, 13.45
56.96, 57.41, 57.36, 57.79, 151.24, 151.72, 149.95, 150.39, 56.24
57.33, 57.75, 57.69, 58.1,  150.06, 150.51, 148.82, 149.23, 24.52
56.65, 57.09, 57.05, 57.47, 150.91, 151.38, 149.63, 150.06, 38.24
57.01, 57.44, 57.38, 57.78, 149.74, 150.18, 148.5,  148.91, 84.25
...

ポリゴンの多くは交差またはオーバーラップしています。ここで、緯度 -90° から 90° まで、経度 -180° から 180° までの範囲の an*m 行列を、たとえば 0.25°x0.25° のステップで作成して、(面積加重) 平均データを保存したいと考えています。各ピクセル内にあるすべてのポリゴンの値。

したがって、通常のメッシュ内の 1 つのピクセルは、1 つまたは複数のポリゴンの平均値を取得します (または、ポリゴンがピクセルとオーバーラップしない場合はなし)。各ポリゴンは、このピクセル内の面積分率に応じて、この平均値に寄与する必要があります。

基本的に、通常のメッシュとポリゴンは次のようになります。

ここに画像の説明を入力

ピクセル 2 を見ると、2 つのポリゴンがこのピクセル内にあることがわかります。したがって、面積分数を考慮して、両方のポリゴンの平均データ値を取得する必要があります。結果は、通常のメッシュ ピクセルに格納する必要があります。

私はウェブを見回しましたが、これまでのところ満足のいくアプローチは見つかりませんでした。私は日常業務で Python/Numpy を使用しているので、それに固執したいと思います。これは可能ですか？パッケージは有望に見えますが、どこから始めればよいかわかりません...すべてをpostgisデータベースに移植するのは大変な労力であり、私の方法にはかなりの障害があると思います.

4

2 に答える 2

3

それを行う方法はたくさんありますが、はい、Shapely が役に立ちます。あなたのポリゴンは四角形のように見えますが、私がスケッチするアプローチはそれを考慮していません. shapely.geometryのbox()とPolygon( )以外は必要ありません。

ピクセルごとに、ピクセル境界を各ポ​​リゴンの最小バウンディング ボックスと比較して、ほぼオーバーラップするポリゴンを見つけます。

from shapely.geometry import box, Polygon

for pixel in pixels:
    # say the pixel has llx, lly, urx, ury values.
    pixel_shape = box(llx, lly, urx, ury)

    for polygon in approximately_overlapping:
        # say the polygon has a ``value`` and a 2-D array of coordinates 
        # [[x0,y0],...] named ``xy``.
        polygon_shape = Polygon(xy)
        pixel_value += polygon_shape.intersection(pixel_shape).area * value

ピクセルとポリゴンが交差しない場合、それらの交差領域は 0 になり、そのピクセルに対するそのポリゴンの寄与は消失します。

于 2012-12-18T17:29:01.853 に答える
1

最初の質問にいくつか追加しましたが、これはこれまでのところ有効な解決策です。物事をスピードアップするためのアイデアはありますか? まだかなり遅いです。入力として、100000 を超えるポリゴンがあり、meshgrid には 720*1440 のグリッド セルがあります。交差するポリゴンのないグリッド セルがたくさんあるため、順序を変更したのもそのためです。さらに、グリッド セルと交差するポリゴンが 1 つだけの場合、グリッド セルはポリゴンのデータ値全体を受け取ります。また、「後処理」部分の面積分率とデータ値を格納する必要があるため、交差可能数を10に設定しました。

from shapely.geometry import box, Polygon
import h5py
import numpy as np

f = h5py.File('data.he5','r')
geo = f['geo'][:] #10 columns: 4xlat, lat center, 4xlon, lon center 
product = f['product'][:]
f.close()

#prepare the regular meshgrid
delta = 0.25
darea = delta**-2
llx, lly = np.meshgrid( np.arange(-180, 180, delta), np.arange(-90, 90, delta) )
urx, ury = np.meshgrid( np.arange(-179.75, 180.25, delta), np.arange(-89.75, 90.25, delta) )
lly = np.flipud(lly)
ury = np.flipud(ury)
llx = llx.flatten()
lly = lly.flatten()
urx = urx.flatten()
ury = ury.flatten()

#initialize the data structures
data = np.zeros(len(llx),'f2')+np.nan
counter = np.zeros(len(llx),'f2')
fraction = np.zeros( (len(llx),10),'f2')
value = np.zeros( (len(llx),10),'f2')

#go through all polygons
for ii in np.arange(1000):#len(hcho)):

    percent = (float(ii)/float(len(hcho)))*100
    print("Polygon: %i (%0.3f %%)" % (ii, percent))

    xy = [ [geo[ii,5],geo[ii,0]], [geo[ii,7],geo[ii,2]], [geo[ii,8],geo[ii,3]], [geo[ii,6],geo[ii,1]] ]
    polygon_shape = Polygon(xy)

    # only go through grid cells which might intersect with the polygon    
    minx = np.min( geo[ii,5:9] )
    miny = np.min( geo[ii,:3] )
    maxx = np.max( geo[ii,5:9] )
    maxy = np.max( geo[ii,:3] )
    mask = np.argwhere( (lly>=miny) & (lly<=maxy) & (llx>=minx) & (llx<=maxx) )
    if mask.size:
        cc = 0
        for mm in mask:
            cc = int(counter[mm])
            pixel_shape = box(llx[mm], lly[mm], urx[mm], ury[mm])
            fraction[mm,cc] = polygon_shape.intersection(pixel_shape).area * darea
            value[mm,cc] = hcho[ii]
            counter[mm] += 1

print("post-processing")
mask = np.argwhere(counter>0)
for mm in mask:
    for cc in np.arange(counter[mm]):
        maxfraction = np.sum(fraction[mm,:])
        value[mm,cc] = (fraction[mm,cc]/maxfraction) * value[mm,cc]
    data[mm] = np.mean(value[mm,:int(counter[mm])])

data = data.reshape( 720, 1440 )
于 2012-12-19T15:52:25.783 に答える