在面内生成坐标 [英] Generate coordinates inside Polygon

问题描述

我想将多边形的值装箱到精细的规则网格中. 例如，我有以下坐标:

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

我的常规网格如下:

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

但是，可能仍然有几个网格像素仍与多边形相交.因此，我想在多边形(data_lons，data_lats)内部生成一堆坐标，并像以前一样找到它们对应的网格像素.您是否建议随机或系统地生成坐标?我失败了，但仍在尝试.

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.

注意:一个数据集包含大约800,000个多边形，因此它必须非常快(几秒钟).这也是为什么我选择这种方法的原因，因为它没有考虑重叠区域...(例如我之前的问题

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)

推荐答案

您需要测试以下方法，以查看其速度是否足够快.首先，您应该将所有经纬度修改为使它们(可能是分数)的索引成为网格:

You'll need to test the following approach to see if it is fast enough. First, you should modify all your lats and lons into, to make them (possibly fractional) indices into your grid:

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)的速度会更快.

Next, we want to split your polygons into triangles. For any convex polygon, you could take the center of the polygon as one vertex of all triangles, and then the vertices of the polygon in consecutive pairs. But if your polygon are all quadrilaterals, it is going to be faster to divide them into only 2 triangles, using vertices 0, 1, 2 for the first, and 0, 2, 3 for the second.

要知道某个点是否在三角形内，我将使用所述的重心坐标方法

To know if a certain point is inside a triangle, I am going to use the barycentric coordinates approach described here. This first function checks whether a bunch of points are inside a triangle:

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)

通过三角形的顶点给出三角形的点，可以通过在三角形的边界框中的晶格点上运行上述函数来获得其晶格点:

Given a triangle by its vertices, you can get the lattice points inside it, by running the above function on the lattice points in the bounding box of the triangle:

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]

对于四边形，只需将最后一个函数调用两次:

And for a quadrilateral, you simply call this last function twice :

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

If you run this code on your example data, you will get two empty arrays returned: that's because the order of the quadrilateral points is a surprising one: indices 0 and 1 define one diagonal, 2 and 3 the other. My function above was expecting the vertices to be ordered around the polygon. If you really are doing things this other way, you need to change the second call to lattice_points_in_triangle inside lattice_points_in_quadrilateral so that the indices used are [0, 1, 3] instead of [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:

If you change the resolution of your grid to 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倍，无法满足您的需求:

Timing wise this approach is going to be, in my system, about 10x too slow for your needs:

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

因此，您正在等待20秒以上.处理您的80,000个多边形.

So you are looking at over 20 sec. to process your 80,000 polygons.

这篇关于在面内生成坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！