生成落在多边形或形状内的点网格的最快方法?
Fastest way to produce a grid of points that fall within a polygon or shape?
我在 python 中使用 shapely 并尝试在最快的 O(n) 时间内在落入形状内的网格中生成均匀间隔的点。形状可以是任何封闭的多边形,而不仅仅是正方形或圆形。我目前的做法是:
- 找到 min/max y & x 来构建一个矩形。
- 在给定间距参数(分辨率)的情况下构建点网格
- 一一验证点是否落在形状内。
有更快的方法吗?
# 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
我在 python 中使用 shapely 并尝试在最快的 O(n) 时间内在落入形状内的网格中生成均匀间隔的点。形状可以是任何封闭的多边形,而不仅仅是正方形或圆形。我目前的做法是:
- 找到 min/max y & x 来构建一个矩形。
- 在给定间距参数(分辨率)的情况下构建点网格
- 一一验证点是否落在形状内。
有更快的方法吗?
# 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