无法生成随机点,直到一个点位于多边形内部将无法正常工作,因为这确实花费了不可预测的时间。
#1 楼
首先将多边形分解为三角形,然后在其内部生成点。 (要获得均匀分布,请按其面积对每个三角形加权。)评论
+1简单有效。值得指出的是,可以在一个三角形内生成均匀的随机点而没有任何排斥,因为在任何三角形与等腰的直角三角形之间都有(容易计算的)保留面积的映射,该三角形为半个正方形,例如一半y坐标超过x坐标。生成两个随机坐标,并对它们进行排序以获得等腰三角形中的随机点,然后将其映射回原始三角形。
– hu
2011-2-22在20:19
+1我真的很喜欢您引用的文章所引用的关于三线性坐标的讨论。我想这将适合于其表面表示为三角形镶嵌的球体。在投影平面上,这不是真正的随机分布,对吗?
– Kirk Kuykendall
2011-2-22在20:42
@whuber-向您+1回复。另一种方法(在链接中,但它们是通过手挥动的)是在共享边上反射来自均匀采样四边形的拒绝点,然后反射回三角形。
– Dan S.
2011-2-22在21:07
@Kirk-引用链接是一种反帮助,因为它在“正确”方式之前列出了一堆错误的(非均匀的)采样方法,包括三线性坐标。看起来没有直接的方法可以获取具有三线性坐标的均匀采样。通过将3d中的随机单位向量转换为等效的经/纬度,我将对整个球面进行均匀采样,但这就是我。 (不确定采样是否限于球形三角形/多边形)(也不确定例如wgs84上的真正均匀采样:我认为只是拾取角度会偏向极点。)
– Dan S.
2011-2-22在21:15
@Dan要均匀采样球体,请使用圆柱等面积投影(坐标是经度和纬度的余弦)。如果要在不使用投影的情况下进行采样,则有一个漂亮的技巧:生成三个独立的标准正态变量(x,y,z)并将它们投影到点(Rx / n,Ry / n,R * z / n) n ^ 2 = x ^ 2 + y ^ 2 + z ^ 2,R是地球半径。如果需要的话,转换为(lat,lon)(在椭球上工作时使用真实纬度)。之所以起作用,是因为该三元正态分布是球对称的。对于三角形采样,请坚持投影。
– hu
2011-2-22在21:26
#2 楼
当您在此问题上放置QGIS标签时:随机点工具可用于边界层。如果您在寻找代码,则基础插件源代码应该有帮助。
评论
即使是5年后,仍然确实有帮助!
–搁浅的孩子
16年7月12日在12:29
#3 楼
有一些不错的库可以为您完成大部分繁重的工作。在python中使用[shapely] [1]的示例。
import random
from shapely.geometry import Polygon, Point
def get_random_point_in_polygon(poly):
minx, miny, maxx, maxy = poly.bounds
while True:
p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if poly.contains(p):
return p
p = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
point_in_poly = get_random_point_in_polygon(mypoly)
或使用
.representative_point()
在对象内获取一个点(如戴恩提到的那样):返回廉价计算的对象保证在几何对象内的点。
poly.representative_point().wkt
'POINT (-1.5000000000000000 0.0000000000000000)'
[1]: https://shapely.readthedocs.io
评论
它不应该来自shapely.geometry导入...吗?
– PyMapr
17年5月17日在10:48
您也可以使用Representative_point方法:shapely.readthedocs.io/en/latest/…
–戴恩
19年7月4日在15:27
#4 楼
您可以确定多边形的范围,然后将X和Y值的随机数生成限制在这些范围内。基本过程:
1)确定多边形顶点的maxx,maxy,minx,miny,
2)使用这些值作为边界生成随机点
3)测试与多边形相交的每个点,
4)当有足够的点满足相交测试时,停止生成
这里是相交测试的算法(C#):
bool PointIsInGeometry(PointCollection points, MapPoint point)
{
int i;
int j = points.Count - 1;
bool output = false;
for (i = 0; i < points.Count; i++)
{
if (points[i].X < point.X && points[j].X >= point.X || points[j].X < point.X && points[i].X >= point.X)
{
if (points[i].Y + (point.X - points[i].X) / (points[j].X - points[i].X) * (points[j].Y - points[i].Y) < point.Y)
{
output = !output;
}
}
j = i;
}
return output;
}
#5 楼
如果选择R,请参阅?spsample
封装中的sp
。可以从rgdal软件包中内置的任何GDAL支持的格式中读取多边形,然后spsample
可以使用多种采样选项直接在导入的对象上工作。 评论
+1-由于R是开放源代码,如果要复制,您可以随时访问源代码以了解它们是如何完成的。对于点模式,可能也对spatstat软件包中的仿真工具感兴趣。
– Andy W
2011-2-24在14:05
#6 楼
我想提供一种几乎不需要GIS分析的解决方案。特别是,它不需要对任何多边形进行三角剖分。以下算法以伪代码给出,除了基本的列表处理功能(创建,查找长度,附加,排序,提取之外,还引用了一些简单的操作)。子列表和串联),并在间隔[0,1)中生成随机浮点数:
Area: Return the area of a polygon (0 for an empty polygon).
BoundingBox: Return the bounding box (extent) of a polygon.
Width: Return the width of a rectangle.
Height: Return the height of a rectangle.
Left: Split a rectangle into two halves and return the left half.
Right: ... returning the right half.
Top: ... returning the top half.
Bottom: ... returning the bottom half.
Clip: Clip a polygon to a rectangle.
RandomPoint: Return a random point in a rectangle.
Search: Search a sorted list for a target value. Return the index
of the last element less than the target.
In: Test whether a point is inside a polygon.
这些几乎都可以在任何GIS或图形编程环境中使用(并且如果没有,则易于编码)。
Clip
一定不能返回退化的多边形(即面积为零的多边形)。过程
SimpleRandomSample
有效地获取了随机分布在多边形内的点的列表。它是SRS
的包装,可将多边形分成较小的块,直到每个块足够紧凑以进行有效采样为止。为此,它使用预先计算的随机数列表来确定要分配给每个片段的点数。可以通过更改参数
t
来“调整” SRS。这是可以容忍的最大边界框:多边形面积比。使其变小(但大于1)会导致大多数多边形被分成许多部分;增大它可能会导致某些多边形(弯曲,带条或充满孔的多边形)的许多试验点被拒绝。这样可以保证对原始多边形进行采样的最大时间是可以预测的。Procedure SimpleRandomSample(P:Polygon, N:Integer) {
U = Sorted list of N independent uniform values between 0 and 1
Return SRS(P, BoundingBox(P), U)
}
如有必要,下一个过程将递归调用自身。神秘的表达式
t*N + 5*Sqrt(t*N)
保守地估计了需要多少点的上限,这考虑了机会的可变性。失败的可能性仅为每百万个过程调用0.3。如果您愿意,可以将5增大到6甚至7,以减少这种可能性。Procedure SRS(P:Polygon, B:Rectangle, U:List) {
N = Length(U)
If (N == 0) {Return empty list}
aP = Area(P)
If (aP <= 0) {
Error("Cannot sample degenerate polygons.")
Return empty list
}
t = 2
If (aP*t < Area(B)) {
# Cut P into pieces
If (Width(B) > Height(B)) {
B1 = Left(B); B2 = Right(B)
} Else {
B1 = Bottom(B); B2 = Top(B)
}
P1 = Clip(P, B1); P2 = Clip(P, B2)
K = Search(U, Area(P1) / aP)
V = Concatenate( SRS(P1, B1, U[1::K]), SRS(P2, B2, U[K+1::N]) )
} Else {
# Sample P
V = empty list
maxIter = t*N + 5*Sqrt(t*N)
While(Length(V) < N and maxIter > 0) {
Decrement maxIter
Q = RandomPoint(B)
If (Q In P) {Append Q to V}
}
If (Length(V) < N) {
Error("Too many iterations.")
}
}
Return V
}
#7 楼
如果多边形是凸面的并且您知道所有顶点,则可能需要考虑对顶点进行“随机”凸面加权,以采样一个新点,该点必须位于凸面壳内(在这种情况下为多边形)。例如,假设您有一个带有顶点的N边凸多边形
V_i, i={1,..,N}
然后随机生成N个凸权重
w_1,w_2,..,w_N such that ∑ w_i = 1; w_i>=0
然后由
Y= ∑ w_i*V_i
给出随机采样的点。
可以采用不同的方法采样N个凸权值
>
均匀挑选一个范围内的N-1个数字(不替换),对它们进行排序并归一化它们之间的N个间隔以得到权重。
您也可以从经常使用的Dirichlet分布中取样作为多项式分布的共轭先验,这与您的情况下的凸权重相似。
当多边形不是非常严重地非凸时,您可以考虑先将其转换为凸包。这至少应在很大程度上限制位于多边形外部的点数。
#8 楼
使用v.random在GRASS GIS(一个命令)中解决该任务非常容易。下面的示例介绍了如何向选定的多边形中添加3个随机点(此处为罗利市的邮政编码区域) ,NC)。通过修改SQL“ where”语句,可以选择多边形。
评论
请注意,邮政编码是线,而不是多边形。
–理查德
18 Mar 8 '18 at 17:25
你能详细说明吗?在我看来,它也是指以下领域:en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes
– markusN
18年3月12日在22:49
当然可以:邮政编码参考特定的邮局及其邮件传递路线。因此,邮政编码是线,而不是多边形。它们可以彼此重叠,包含漏洞,并且不一定覆盖整个美国或任何给定的州。因此,使用它们划分区域很危险。人口普查单位(如块组)是更好的选择。另请参阅:此和此。
–理查德
18年3月12日在23:14
谢谢!它也可能取决于国家/地区,例如,请参见en.wikipedia.org/wiki/Postal_codes_in_Germany-但是,邮政编码不是我的核心主题,只是想说明和回答原始问题“生成位于多边形内的点”而不是在这里讨论OT的邮政编码定义:-)
– markusN
18年3月13日在13:35
两项均注明。我可能应该写一些博客文章,以便下次可以更简洁地说:-)
–理查德
18年3月13日在16:19
#9 楼
答案链接https://gis.stackexchange.com/a/307204/103524
三种使用不同方法的算法。
Git Repo:https://github.com/I1mran/Postgis-自定义
这是一种简单且最佳的方法,它使用从x和y方向的实际坐标距离。
内部算法使用WGS 1984(4326)和结果变换
插入的SRID。
功能===================================== ==============================
CREATE OR REPLACE FUNCTION public.I_Grid_Point_Distance(geom public.geometry, x_side decimal, y_side decimal)
RETURNS public.geometry AS $BODY$
DECLARE
x_min decimal;
x_max decimal;
y_max decimal;
x decimal;
y decimal;
returnGeom public.geometry[];
i integer := -1;
srid integer := 4326;
input_srid integer;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
geom := ST_SetSRID(geom, srid);
----RAISE NOTICE 'No SRID Found.';
ELSE
----RAISE NOTICE 'SRID Found.';
END CASE;
input_srid:=st_srid(geom);
geom := st_transform(geom, srid);
x_min := ST_XMin(geom);
x_max := ST_XMax(geom);
y_max := ST_YMax(geom);
y := ST_YMin(geom);
x := x_min;
i := i + 1;
returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
<<yloop>>
LOOP
IF (y > y_max) THEN
EXIT;
END IF;
CASE i WHEN 0 THEN
y := ST_Y(returnGeom[0]);
ELSE
y := ST_Y(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), y_side, radians(0))::geometry);
END CASE;
x := x_min;
<<xloop>>
LOOP
IF (x > x_max) THEN
EXIT;
END IF;
i := i + 1;
returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
x := ST_X(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), x_side, radians(90))::geometry);
END LOOP xloop;
END LOOP yloop;
RETURN
ST_CollectionExtract(st_transform(ST_Intersection(st_collect(returnGeom), geom), input_srid), 1);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE;
通过简单查询使用该函数,几何必须有效和多边形,多面或信封
SELECT I_Grid_Point_Distance(geom, 50, 61) from polygons limit 1;
结果=============================== =======================================
基于NicklasAvén算法的第二个函数。尝试处理任何SRID。
我对算法进行了以下更改。
x和y方向的单独变量用于像素大小,
用于计算距离的新变量
输入任何SRID,对Geom进行功能变换到
椭球或基准椭球基准面的工作环境中,然后将距离应用于每一侧,获取结果,然后转换为输入SRID。
功能================================= ===================================
CREATE OR REPLACE FUNCTION I_Grid_Point(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
geom := ST_SetSRID(geom, srid);
RAISE NOTICE 'SRID Not Found.';
ELSE
RAISE NOTICE 'SRID Found.';
END CASE;
CASE spheroid WHEN false THEN
RAISE NOTICE 'Spheroid False';
srid := 4326;
x_side := x_side / 100000;
y_side := y_side / 100000;
else
srid := 900913;
RAISE NOTICE 'Spheroid True';
END CASE;
input_srid:=st_srid(geom);
geom := st_transform(geom, srid);
x_max := ST_XMax(geom);
y_max := ST_YMax(geom);
x_min := ST_XMin(geom);
y_min := ST_YMin(geom);
RETURN QUERY
WITH res as (SELECT ST_SetSRID(ST_MakePoint(x, y), srid) point FROM
generate_series(x_min, x_max, x_side) as x,
generate_series(y_min, y_max, y_side) as y
WHERE st_intersects(geom, ST_SetSRID(ST_MakePoint(x, y), srid))
) select ST_TRANSFORM(ST_COLLECT(point), input_srid) from res;
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;
将其与简单查询一起使用。
SELECT I_Grid_Point(geom, 22, 15, false) from polygons;
结果======================================= ===========================
基于串联发生器的函数。
功能================================================= =================
CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
x_series DECIMAL;
y_series DECIMAL;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
geom := ST_SetSRID(geom, srid);
RAISE NOTICE 'SRID Not Found.';
ELSE
RAISE NOTICE 'SRID Found.';
END CASE;
CASE spheroid WHEN false THEN
RAISE NOTICE 'Spheroid False';
else
srid := 900913;
RAISE NOTICE 'Spheroid True';
END CASE;
input_srid:=st_srid(geom);
geom := st_transform(geom, srid);
x_max := ST_XMax(geom);
y_max := ST_YMax(geom);
x_min := ST_XMin(geom);
y_min := ST_YMin(geom);
x_series := CEIL ( @( x_max - x_min ) / x_side);
y_series := CEIL ( @( y_max - y_min ) / y_side );
RETURN QUERY
SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM
generate_series(0, x_series) as x,
generate_series(0, y_series) as y
WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;
通过简单查询使用它。
SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons;
结果======= ================================================== =================
评论
相反,时间是可以预测的。它与多边形的范围面积除以多边形的面积之比成正比,再乘以生成和测试单个点所需的时间。时间稍有变化,但是变化与点数的平方根成正比。对于较大的数字,这变得无关紧要。如有必要,可以将弯曲的多边形分成更紧凑的碎片,以使该面积比降低到较低的值。只有分形多边形会给您带来麻烦,但我怀疑您是否拥有它们!使用Quantum Gis可能与功能有关的分布点重复
还有:gis.stackexchange.com/questions/4663/…
@Pablo:很好的发现。但是,这两个问题都是特定于软件的,并且都涉及在多边形内放置点的规则阵列,而不是随机点
一致的观点是多边形内的随机点与规则点生成。