我具有多边形特征,并希望能够在其中生成点。我需要这项功能来完成一项分类任务。

无法生成随机点,直到一个点位于多边形内部将无法正常工作,因为这确实花费了不可预测的时间。

评论

相反,时间是可以预测的。它与多边形的范围面积除以多边形的面积之比成正比,再乘以生成和测试单个点所需的时间。时间稍有变化,但是变化与点数的平方根成正比。对于较大的数字,这变得无关紧要。如有必要,可以将弯曲的多边形分成更紧凑的碎片,以使该面积比降低到较低的值。只有分形多边形会给您带来麻烦,但我怀疑您是否拥有它们!

使用Quantum Gis可能与功能有关的分布点重复

还有:gis.stackexchange.com/questions/4663/…

@Pablo:很好的发现。但是,这两个问题都是特定于软件的,并且都涉及在多边形内放置点的规则阵列,而不是随机点

一致的观点是多边形内的随机点与规则点生成。

#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;
结果======= ================================================== =================