您需要测试以下方法以查看它是否足够快。首先,您应该将所有 lats 和 lons 修改为,以使它们(可能是分数)索引到您的网格中:
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_triangle
inside lattice_points_in_quadrilateral
,以便使用的索引[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 个多边形。