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

Fastest way to produce a grid of points that fall within a polygon or shape?

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

  1. 找到 min/max 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)])

我能想到的最好的办法就是这样做:

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

哦,为什么是的。使用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)))

我看到你回答了你的问题(并且似乎对使用交集很满意)但也注意到 shapely(和底层 geos 库)有 准备好的几何图形 用于对某些谓词进行更高效的批处理操作(包含、contains_properly、覆盖和相交)。 参见 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))

如果您想在 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