在Polygon中生成坐标

Question

我想将多边形的值组合成一个精细的规则网格。 举例来说，我有以下坐标：

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]) 

我的规则网格是这样的：

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

现在，我可以发现，很容易对应到我的多边形的中心网格像素：

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 

但是，可能有一些网格像素仍然与多边形相交。因此，我想在我的多边形（data_lons，data_lats）中生成一堆坐标，并像以前一样查找它们相应的网格像素。你有没有建议随机或系统地生成坐标？我失败了，但我仍然在努力。

注意：一个数据集包含约80000个多边形，所以它必须非常快（几秒钟）。这也是为什么我选择这种方法，因为它没有考虑重叠区域...（就像我之前的问题Data binning: irregular polygons to regular mesh这非常慢）

Answer 1

您需要测试以下方法以查看它是否足够快。首先，你应该修改所有的背阔肌和离子吸附进去，让他们（可能是部分的）指标到您的网格：

idx_lats = (data_lats - lat_grid_start)/lat_grid step 
idx_lons = (data_lons - lon_grid_start)/lon_grid step 

接下来，我们希望您的多边形分割成三角形。对于任何凸多边形，可以将多边形的中心作为所有三角形的一个顶点，然后将多边形的顶点连续成对。但是如果你的多边形都是四边形的话，将它们分成只有两个三角形的速度会更快，第一个使用顶点0,1,2，第二个使用0,2,3。

要知道某个点是否在三角形内，我将使用here中描述的重心坐标方法。第一个函数检查是否一堆点是三角形内：

def check_in_triangle(x, y, x_tri, y_tri) : 
    A = np.vstack((x_tri[0], y_tri[0])) 
    lhs = np.vstack((x_tri[1:], y_tri[1:])) - A 
    rhs = np.vstack((x, y)) - A 
    uv = np.linalg.solve(lhs, rhs) 
    # Equivalent to (uv[0] >= 0) & (uv[1] >= 0) & (uv[0] + uv[1] <= 1) 
    return np.logical_and(uv >= 0, axis=0) & (np.sum(uv, axis=0) <= 1) 

鉴于其顶点的三角形，你可以得到的格点里面，通过在边界框格点运行上面的功能

def lattice_points_in_triangle(x_tri, y_tri) : 
    x_grid = np.arange(np.ceil(np.min(x_tri)), np.floor(np.max(x_tri)) + 1) 
    y_grid = np.arange(np.ceil(np.min(y_tri)), np.floor(np.max(y_tri)) + 1) 
    x, y = np.meshgrid(x_grid, y_grid) 
    x, y = x.reshape(-1), y.reshape(-1) 
    idx = check_in_triangle(x, y, x_tri, y_tri) 
    return x[idx], y[idx] 

而对于一个四边形，你只需调用这个函数最后两次：三角形的

def lattice_points_in_quadrilateral(x_quad, y_quad) : 
    return map(np.concatenate, 
       zip(lattice_points_in_triangle(x_quad[:3], y_quad[:3]), 
        lattice_points_in_triangle(x_quad[[0, 2, 3]], 
               y_quad[[0, 2, 3]]))) 

如果您在示例数据运行此代码，哟你会得到两个返回的空数组：这是因为四边形点的顺序令人吃惊：索引0和1定义了一个对角线，而另一个则定义为2和3。我上面的函数期望在多边形周围排列顶点。如果您确实是以其他方式进行操作，则需要将lattice_points_in_quadrilateral中的第二个呼叫更改为lattice_points_in_triangle，以便使用的索引是[0, 1, 3]而不是[0, 2, 3]。

而现在，这种变化：

>>> idx_lats = (data_lats - (-180))/0.25 
>>> idx_lons = (data_lons - (-90))/0.25 
>>> lattice_points_in_quadrilateral(idx_lats, idx_lons) 
[array([952]), array([955])] 

如果您改变了网格的分辨率为0。1：

>>> idx_lats = (data_lats - (-180))/0.1 
>>> idx_lons = (data_lons - (-90))/0.1 
>>> lattice_points_in_quadrilateral(idx_lats, idx_lons) 
[array([2381, 2380, 2381, 2379, 2380, 2381, 2378, 2379, 2378]), 
array([2385, 2386, 2386, 2387, 2387, 2387, 2388, 2388, 2389])] 

定时明智的这一做法将是，在我的系统，约10倍的速度太慢您的需求：

In [8]: %timeit lattice_points_in_quadrilateral(idx_lats, idx_lons) 
1000 loops, best of 3: 269 us per loop 

那么您正在寻找超过20秒。处理您的80,000个多边形。

Answer 2

我通过简单计算角点像素之间的坐标来处理快速而肮脏的解决方案。看看：

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

因此，红点代表原来的角 - 蓝色的代表它们之间的中间像素。

我能做到这一次，包括中心像素（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])) 

这工作真的很好，并且我可以将每个点直接指定给其对应的网格像素，如下所示：

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

但是最终的分档现在需要〜7秒！我怎样才能加速这一代码了：

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