在多边形内生成坐标

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个多边形，所以这个过程必须非常快（几秒钟内完成）。这也是我选择这种方法的原因，因为它不考虑重叠的区域……（就像我之前的问题数据分箱：不规则多边形到规则网格，那种方法非常慢）

Answer 1

我做了一个简单粗暴的解决方案，就是通过计算角落像素之间的坐标来实现。看看这个：

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

Answer 2

你需要测试一下下面的方法，看它是否足够快。首先，你应该把所有的经纬度转换成网格的索引，可能会是小数：

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来构成第二个三角形。

要判断一个点是否在三角形内部，我会使用重心坐标的方法，详细介绍可以参考这里。这个函数首先检查一组点是否在三角形内：

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

所以处理你的80,000个多边形大约需要超过20秒。