binning dei dati: poligoni irregolari su mesh normale

Question

Ho migliaia di poligoni memorizzati in un formato tabella (date le loro 4 coordinate d'angolo) che rappresentano piccole regioni della terra. Inoltre, ciascun poligono ha un valore di dati. Il file è per esempio come questo:

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

Molti dei poligoni si intersecano o si sovrappongono. Ora vorrei creare una matrice * m che vada da -90 ° a 90 ° di latitudine e da -180 ° a 180 ° di longitudine in passi di, ad esempio, 0,25 ° x0,25 ° per memorizzare i dati medi (ponderati per l'area) valore di tutti i poligoni che rientrano in ciascun pixel.

Quindi, un pixel nella mesh normale deve ottenere il valore medio di uno o più poligoni (o nessuno se nessun poligono si sovrappone al pixel). Ogni poligono dovrebbe contribuire a questo valore medio in base alla sua frazione di area all'interno di questo pixel.

Fondamentalmente la maglia regolare ei poligoni simile a questa:

Se si guarda al pixel 2, si vede che due poligoni sono dentro questo pixel. Quindi, devo prendere il valore medio dei dati di entrambi i poligoni considerando le loro frazioni di area. Il risultato dovrebbe quindi essere memorizzato nel normale pixel mesh.

Ho guardato intorno al web e non ho trovato nessun approccio soddisfacente fino ad ora. Dal momento che sto usando Python/Numpy per il lavoro quotidiano mi piacerebbe attenervisi. È possibile? Il pacchetto shapely sembra promettente, ma non so da dove cominciare ... Portare tutto su un database PostGG è un enorme sforzo e immagino ci saranno molti ostacoli sulla mia strada.

Answer 1

Ci sono molti modi per farlo, ma sì, Shapely può aiutare. Sembra che i tuoi poligoni siano quadrilateri, ma l'approccio che abbozzerò non conta su questo. Non avrai bisogno di altro che di box() e di Polygon() da shapely.geometry.

Per ogni pixel, trovare i poligoni che approssimativamente si sovrappongono confrontando i limiti dei pixel con il riquadro di delimitazione minimo di ciascun poligono.

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 

Se il pixel e poligono non si intersecano, l'area della loro intersezione sarà 0 e il contributo del poligono a quel pixel svanisce.

Answer 2

Ho aggiunto un paio di cose alla mia domanda iniziale, ma questa è una soluzione funzionante finora. Hai qualche idea per velocizzare le cose? È ancora piuttosto lento. Come input, ho oltre 100000 poligoni e meshgrid ha celle a griglia 720 * 1440. Questo è anche il motivo per cui ho cambiato l'ordine, perché ci sono molte celle della griglia senza poligoni intersecanti. Inoltre, quando c'è un solo poligono che si interseca con una cella della griglia, la cella della griglia riceve l'intero valore dei dati del poligono. Inoltre, dal momento che ho per memorizzare la frazione di area e il valore di dati per la parte "post-processing", ho impostare il numero possibile di intersezioni a 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)