डेटा बिनिंग: नियमित जाल

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

बहुभुज में से कई एक दूसरे को काटना या ओवरलैप। अब मैं -90 डिग्री से 9 0 डिग्री अक्षांश और -180 डिग्री से 180 डिग्री अक्षांश के चरणों में एक * एम मैट्रिक्स बनाना चाहता हूं, उदाहरण के लिए, 0.25 डिग्री x0.25 डिग्री (क्षेत्र-भारित) डेटा डेटा को स्टोर करने के लिए प्रत्येक पिक्सेल के भीतर आने वाले सभी बहुभुजों का मूल्य।

तो, नियमित जाल में एक पिक्सेल को एक या अधिक बहुभुज का औसत मूल्य प्राप्त होगा (या कोई भी पॉलीगॉन पिक्सेल के साथ ओवरलैप नहीं होता है)। प्रत्येक बहुभुज को इस पिक्सेल के भीतर अपने क्षेत्र के अंश के आधार पर इस औसत मूल्य में योगदान देना चाहिए।

मूल रूप से नियमित रूप से जाल और बहुभुज इस तरह दिखेगा:

आप पिक्सेल 2 को देखें, तो आप देखते हैं कि दो बहुभुज इस पिक्सेल के अंदर हैं। इस प्रकार, मुझे अपने क्षेत्र के अंशों पर विचार करते हुए दोनों बहुभुजों का औसत डेटा मान लेना होगा। परिणाम नियमित जाल पिक्सेल में संग्रहीत किया जाना चाहिए।

मैंने वेब के चारों ओर देखा और अब तक इसके लिए कोई संतोषजनक दृष्टिकोण नहीं मिला है। चूंकि मैं दैनिक काम के लिए पाइथन/न्यूम्पी का उपयोग कर रहा हूं, इसलिए मैं इसके साथ रहना चाहता हूं। क्या यह संभव है? पैकेज shapely वादा करता है लेकिन मुझे नहीं पता कि कहां से शुरू करना है ... पोस्टगिस डेटाबेस में सबकुछ पोर्ट करना एक बहुत ही कठिन प्रयास है और मुझे लगता है कि मेरे रास्ते में कुछ बाधाएं होंगी।

Answer 1

ऐसा करने के कई तरीके हैं, लेकिन हाँ, आकार में मदद कर सकते हैं। ऐसा प्रतीत होता है कि आपके बहुभुज चतुर्भुज हैं, लेकिन जिस दृष्टिकोण को मैं स्केच करूंगा उस पर भरोसा नहीं करता है। आपको 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 से अधिक बहुभुज हैं और मेषग्रीड में 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)