データビニング：不規則なポリゴンから規則的なメッシュ

Question

地球の小さな領域を表すテーブル形式（四隅の座標が与えられている）で何千ものポリゴンが格納されています。さらに、各ポリゴンにはデータ値があります。

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°の経度から、例えば0.25°x0.25°のステップで（面積加重）平均データを格納する* m行列を作成したいと思います各ピクセルに含まれるすべてのポリゴンの値。

したがって、通常のメッシュの1つのピクセルは、1つ以上のポリゴンの平均値を取得する（または、ポリゴンがピクセルと重なっていない場合はなし）。各ポリゴンは、このピクセル内の面積率によってこの平均値に寄与するはずです。

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

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

私はウェブの周りを見渡し、これまでのところ満足のいくアプローチは見つけられませんでした。私はPython/Numpyを毎日の作業に使用しているので、それに固執したいと思います。これは可能ですか？パッケージshapelyは有望そうですが、私はどこから始めるべきか分かりません。 ポストギズのデータ​​ベースへのすべての移植は非常に手間がかかり、途中でかなりの障害があると思います。

Answer 1

これにはさまざまな方法がありますが、確かに助けてください。あなたのポリゴンは四角形であるようですが、私がスケッチするアプローチはそれを考慮しません。あなたは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になり、そのピクセルへのポリゴンの寄与は消滅します。

Answer 2

最初の質問に2つのことを追加しましたが、これはこれまでの解決策です。あなたは物事をスピードアップするためのアイデアはありますか？それはまだかなり遅いです。入力として、私は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)