Podział danych: nieregularne wielokąty na regularną siatkę

Question

Mam tysiące wielokątów zapisanych w formacie tabeli (z uwzględnieniem ich 4 współrzędnych narożnych), które reprezentują małe obszary ziemi. Ponadto każdy wielokąt ma wartość danych. Plik wygląda na przykład tak:

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

Wiele wielokątów przecinają lub pokrywają. Teraz chciałbym stworzyć matrycę * m w zakresie od -90 ° do 90 ° szerokości geograficznej i -180 ° do 180 ° długości geograficznej w krokach, na przykład, 0,25 ° x 0,25 ° do przechowywania (średniej ważonej) średniej danych wartość wszystkich wielokątów mieszczących się w każdym pikselu.

Tak więc, jeden piksel w siatce regularnej powinien otrzymać średnią wartość jednego lub więcej wielokątów (lub brak, jeśli wielokąt nie pokrywa się z pikselem). Każdy wielokąt powinien przyczyniać się do tej średniej wartości w zależności od jej frakcji obszaru w obrębie tego piksela.

Zasadniczo regularna siatka i wielokąty wyglądać następująco:

Jeśli spojrzeć na piksel 2, widać, że dwa wielokąty są wewnątrz tego piksela. Tak więc, muszę wziąć średnią wartość danych obu poligonów biorąc pod uwagę ich ułamki powierzchni. Wynik powinien być następnie zapisany w zwykłym pikselach siatki.

Rozejrzałem się po internecie i nie znalazłem dotychczas zadowalającego rozwiązania. Ponieważ używam Pythona/Numpy do codziennej pracy, chciałbym się do tego przyłączyć. czy to możliwe? Pakiet shapely wygląda obiecująco, ale nie wiem, od czego zacząć ... Przeniesienie wszystkiego do bazy danych Postgis to ogromny wysiłek i sądzę, że na mojej drodze będzie sporo przeszkód.

Answer 1

Istnieje wiele sposobów na zrobienie tego, ale tak, Shapely może pomóc. Wygląda na to, że twoje wielokąty są czworokątne, ale podejście, które naszkicuję, nie liczy się z tym. Nie będziesz potrzebować niczego poza box() i Polygon() z shapely.geometry.

Dla każdego piksela znajdź wielokątów, które pokrywają się w przybliżeniu, porównując granice pikseli do minimalnego obwiedni każdego wielokąta.

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 

Jeśli piksel i wielokąt nie przecinają się, obszar ich przecięcia będzie 0 i wkład tego wielokąta do tego piksela znika.

Answer 2

Dodałem kilka rzeczy do mojego wstępnego pytania, ale jest to działające rozwiązanie do tej pory. Czy masz jakieś pomysły na przyspieszenie? Jest wciąż dość powolny. Jako dane wejściowe mam ponad 100000 wielokątów, a siatka mesh ma 720 * 1440 komórek siatki. Właśnie dlatego zmieniłem kolejność, ponieważ istnieje wiele komórek siatki bez przecinających się wielokątów. Ponadto, gdy istnieje tylko jeden wielokąt przecinający się z komórką siatki, komórka siatki odbiera całą wartość danych wielokąta. Ponadto, ponieważ muszę przechowywać część obszaru i wartość danych dla „post-processingu” część, ustawić możliwą liczbę przecięć do 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)