生成落在多边形或形状内的点网格的最快方法？

Question

我在 python 中使用 shapely，并尝试在最快的 O(n) 时间内在网格中生成均匀分布的点。形状可以是任何封闭的多边形，而不仅仅是正方形或圆形。我目前的做法是：

1. 找到最小/最大 y 和 x 以构建一个矩形。
2. 在给定间距参数（分辨率）的情况下构建点网格
3. 逐一验证点是否在形状内。

有更快的方法吗？

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape
valid_points.extend([i for i in points if polygon.contains(i)])

Answer 1

我看到您回答了您的问题（并且似乎对使用交集很满意），但也请注意 shapely（以及底层的 geos 库）为更高效的批处理准备了几何图形对某些谓词（contains、contains_properly、covers 和 intersects）的操作。 见Prepared geometry operations。

改编自您问题中的代码，可以这样使用：

from shapely.prepared import prep

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# create prepared polygon
prep_polygon = prep(polygon)

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape using
# the prepared polygon
valid_points.extend(filter(prep_polygon.contains, points))

Answer 2

我能想到的最好的就是这样做：

X,Y = np.meshgrid(np.arange(latmin, latmax, resolution),
              np.arange(lonmin, lonmax, resolution))

#create a iterable with the (x,y) coordinates
points = zip(X.flatten(),Y.flatten())

valid_points.extend([i for i in points if polygon.contains(i)])

Answer 3

哦，为什么是的。使用 shapely 的交集方法。

polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# construct rectangle of points
x, y = np.round(np.meshgrid(np.arange(latmin, latmax, resolution), np.arange(lonmin, lonmax, resolution)),4)
points = MultiPoint(list(zip(x.flatten(),y.flatten())))

# validate each point falls inside shapes
valid_points.extend(list(points.intersection(polygon)))

Answer 4

如果你想在shapely.geometry.Polygon 中生成n 点，有一个简单的迭代函数可以做到。管理tol（容差）参数以加快积分生成。

import numpy as np
from shapely.geometry import Point, Polygon


def gen_n_point_in_polygon(self, n_point, polygon, tol = 0.1):
    """
    -----------
    Description
    -----------
    Generate n regular spaced points within a shapely Polygon geometry
    -----------
    Parameters
    -----------
    - n_point (int) : number of points required
    - polygon (shapely.geometry.polygon.Polygon) : Polygon geometry
    - tol (float) : spacing tolerance (Default is 0.1)
    -----------
    Returns
    -----------
    - points (list) : generated point geometries
    -----------
    Examples
    -----------
    >>> geom_pts = gen_n_point_in_polygon(200, polygon)
    >>> points_gs = gpd.GeoSeries(geom_pts)
    >>> points_gs.plot()
    """
    # Get the bounds of the polygon
    minx, miny, maxx, maxy = polygon.bounds    
    # ---- Initialize spacing and point counter
    spacing = polygon.area / n_point
    point_counter = 0
    # Start while loop to find the better spacing according to tolérance increment
    while point_counter <= n_point:
        # --- Generate grid point coordinates
        x = np.arange(np.floor(minx), int(np.ceil(maxx)), spacing)
        y = np.arange(np.floor(miny), int(np.ceil(maxy)), spacing)
        xx, yy = np.meshgrid(x,y)
        # ----
        pts = [Point(X,Y) for X,Y in zip(xx.ravel(),yy.ravel())]
        # ---- Keep only points in polygons
        points = [pt for pt in pts if pt.within(polygon)]
        # ---- Verify number of point generated
        point_counter = len(points)
        spacing -= tol
    # ---- Return
    return points