Generate coordinates inside Polygon

Question

I want to bin the values of polygons to a fine regular grid. For instance, I have the following coordinates:

data = 2.353
data_lats = np.array([57.81000137,  58.15999985,  58.13000107,  57.77999878])
data_lons = np.array([148.67999268,  148.69999695,  148.47999573,  148.92999268])

My regular grid looks like this:

delta = 0.25
grid_lons = np.arange(-180, 180, delta)
grid_lats = np.arange(90, -90, -delta)
llx, lly = np.meshgrid( grid_lons, grid_lats )
rows = lly.shape[0]
cols = llx.shape[1]
grid = np.zeros((rows,cols))

Now I can find the grid pixel that corresponds to the center of my polygon very easily:

centerx, centery = np.mean(data_lons), np.mean(data_lats)
row = int(np.floor( centery/delta ) + (grid.shape[0]/2))
col = int(np.floor( centerx/delta ) + (grid.shape[1]/2))
grid[row,col] = data

However, there are probably a couple of grid pixels that still intersect with the polygon. Hence, I would like to generate a bunch of coordinates inside my polygon (data_lons, data_lats) and find their corresponding grid pixel as before. Do you a suggestion to generate the coordinates randomly or systematically? I failed, but am still trying.

Note: One data set contains around ~80000 polygons so it has to be really fast (a couple of seconds). That is also why I chose this approach, because it does not account the area of overlap... (like my earlier question Data binning: irregular polygons to regular mesh which is VERY slow)

Answer 1

I worked on a quick and dirty solution by simply calculating the coordinates between corner pixels. Take a look:

dlats = np.zeros((data_lats.shape[0],4))+np.nan
dlons = np.zeros((data_lons.shape[0],4))+np.nan
idx = [0,1,3,2,0] #rearrange the corner pixels

for cc in range(4):
    dlats[:,cc] = np.mean((data_lats[:,idx[cc]],data_lats[:,idx[cc+1]]), axis=0)
    dlons[:,cc] = np.mean((data_lons[:,idx[cc]],data_lons[:,idx[cc+1]]), axis=0)

data_lats = np.column_stack(( data_lats, dlats ))
data_lons = np.column_stack(( data_lons, dlons ))

Thus, the red dots represent the original corners - the blue ones the intermediate pixels between them.

I can do this one more time and include the center pixel (geo[:,[4,9]])

dlats = np.zeros((data.shape[0],8))
dlons = np.zeros((data.shape[0],8))

for cc in range(8):
    dlats[:,cc] = np.mean((data_lats[:,cc], geo[:,4]), axis=0)
    dlons[:,cc] = np.mean((data_lons[:,cc], geo[:,9]), axis=0)

data_lats = np.column_stack(( data_lats, dlats, geo[:,4] ))
data_lons = np.column_stack(( data_lons, dlons, geo[:,9] ))

This works really nice, and I can assign each point directly to its corresponding grid pixel like this:

row = np.floor( data_lats/delta ) + (llx.shape[0]/2)
col = np.floor( data_lons/delta ) + (llx.shape[1]/2)

However the final binning now takes ~7sec!!! How can I speed this code up:

for ii in np.arange(len(data)):
    for cc in np.arange(data_lats.shape[1]):
        final_grid[row[ii,cc],col[ii,cc]] += data[ii]
        final_grid_counts[row[ii,cc],col[ii,cc]] += 1