數據裝箱：不規則多邊形到常規網格

Question

我有成千上萬個以表格格式存儲的多邊形（給定它們的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 
... 

許多多邊形的交叉或重疊。現在，我想創建一個* m矩陣，範圍從-90°到90°緯度和-180°到180°經度，例如，以0.25°x0.25°的步長存儲（面積加權）平均數據落在每個像素內的所有多邊形的值。

因此，在常規網格一個像素應獲得一個或多個多邊形的平均值（或根本沒有，如果沒有多邊形像素重疊）。每個多邊形應根據其在此像素內的面積分數貢獻此平均值。

基本上常規的網格和多邊形是這樣的：

如果你看看像素2，你會看到兩個多邊形此像素內部。因此，我必須考慮它們的面積分數的兩個多邊形的平均數據值。結果應該存儲在常規網格像素中。

我看着網頁四周，發現這個沒有令人滿意的方法爲止。由於我在日常工作中使用Python/Numpy，我想堅持下去。這可能嗎？該軟件包shapely看起來很有希望，但我不知道在哪裏有... 一切移植到PostGIS的數據庫的努力一個可怕的量開始，我想還會有我的方式相當多的障礙。

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

我加了幾件事情對我最初的問題，但這是一個工作的解決方案爲止。你有什麼想法來加快速度嗎？它仍然很慢。作爲輸入，我有超過10萬個多邊形，網格有720 * 1440個網格單元。這也是我改變順序的原因，因爲有很多網格單元沒有相交的多邊形。此外，當只有一個與網格單元相交的多邊形時，網格單元會接收該多邊形的整個數據值。 另外，因爲我有存儲區部分，爲「後處理」部分的數據值，我交叉口的可能數量設置爲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)