数据分箱：不规则多边形到规则网格

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
...

许多多边形相交或重叠。现在我想创建一个 n*m 矩阵，范围从 -90° 到 90° 纬度和 -180° 到 180° 经度，例如 0.25°x0.25° 以存储（面积加权）平均数据落在每个像素内的所有多边形的值。

因此，规则网格中的一个像素应获得一个或多个多边形的平均值（如果没有多边形与像素重叠，则没有）。每个多边形都应根据其在该像素内的面积分数对该平均值做出贡献。

基本上规则网格和多边形看起来像这样：

如果您查看像素 2，您会看到该像素内有两个多边形。因此，考虑到它们的面积分数，我必须取两个多边形的平均数据值。然后将结果存储在常规网格像素中。

我环顾网络，到目前为止还没有找到令人满意的方法。因为我在日常工作中使用 Python/Numpy，所以我想坚持下去。这可能吗？包裹shapely 看起来很有希望，但我不知道从哪里开始...... 将所有内容都移植到 postgis 数据库是一项非常艰巨的工作，我想在我的道路上会有很多障碍。

Answer 1

有很多方法可以做到这一点，但是，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，并且该多边形对该像素的贡献消失。

Answer 2

我在最初的问题中添加了一些内容，但到目前为止这是一个可行的解决方案。你有什么想法可以加快速度吗？它仍然很慢。作为输入，我有超过 100000 个多边形，meshgrid 有 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 )