Биннинг данных: нерегулярные многоугольники для регулярной сетки

У меня есть тысячи полигонов, хранящихся в формате таблицы (с учетом их четырех угловых координат), которые представляют собой небольшие области земли. Кроме того, каждый многоугольник имеет значение данных. Файл выглядит следующим образом:

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 - это ужасное усилие, и я предполагаю, что на моем пути будет немало препятствий.

Ответ 1

Есть много способов сделать это, но да, Shapely может помочь. Кажется, что ваши многоугольники четырехугольные, но подход, который я нарисую, не учитывает этого. Вам не понадобится ничего, кроме box() и Polygon() из shapely.geometry.

Для каждого пикселя найдите многоугольники, которые приблизительно перекрываются с ним, сравнивая пиксели с минимальным ограничивающим прямоугольником каждого многоугольника.

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, а вклад этого многоугольника в этот пиксель исчезнет.

Ответ 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 )