在PostGIS中的多边形内部创建规则点网格

如何在PostGIS中在多边形内部创建点距x，y的规则网格？
例如示例：

我试图制作将此代码与“ postGIS in action”切割代码合并的剪切多边形，但仅创建了一个多边形...我忘记了什么吗？创建或替换功能makegrid（geometry，integer，integer）返回geometry AS'SELECT st_intersection（g1.geom1，g2.geom2）AS geom FROM（SELECT $ 1 AS geom1）AS g1 INNER JOIN（选择st_setsrid（CAST（ST_MakeBox2d（st_setsrid（ ST_Point（x，y），$ 3），st_setsrid（ST_Point（x + $ 2，y + $ 2），$ 3））作为几何），$ 3）as geom2 FROM generate_series（floor（st_xmin（$ 1））:: int，ceiling（ st_xmax（$ 1））:: int，$ 2）as x，generate_series（floor（st_ymin（$ 1））:: int，ceiling（st_ymax（

看我详细的答案。

#1 楼

您可以使用generate_series做到这一点。
如果您不想手动编写网格开始和停止的位置，最简单的方法就是创建一个函数。
我没有正确测试以下内容，但是我认为它应该工作：

 CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_POINT(x, y)) FROM 
generate_series(floor(ST_XMIN())::int, ceiling(ST_XMAX()-ST_XMIN())::int, ) AS x,
generate_series(floor(ST_YMIN())::int, ceiling(ST_YMAX()-ST_YMIN())::int, ) AS y 
WHERE st_intersects(, ST_POINT(x, y))'
LANGUAGE sql

要使用它，您可以执行以下操作：

 SELECT makegrid(the_geom, 1000) from mytable;

其中第一个参数是您希望网格所在的多边形，第二个参数是网格中点之间的距离。
如果需要一个点每行只使用ST_Dump，例如：

 SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

您可能需要将st_setSRID（）添加到st_point函数，否则st_intersects不起作用。

– JaakL
2012年10月31日在9:25

将我的测试版本添加为单独的答案。

– JaakL
2012年10月31日在16:43

#2 楼

我选择了NicklasAvénmakegrid函数代码，并通过读取和使用多边形几何图形中的斜线使其变得更加通用。否则，使用具有已定义srid的多边形会产生错误。
函数：

 CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x, y), ST_SRID())) FROM
generate_series(floor(ST_XMIN())::int, ceiling(ST_XMAX() - ST_XMIN())::int, ) AS x,
generate_series(floor(ST_YMIN())::int, ceiling(ST_YMAX() - ST_YMIN())::int, ) AS y
WHERE st_intersects(, ST_SetSRID(ST_POINT(x,y), ST_SRID()))'
LANGUAGE sql

要使用函数是完全按照NicklasAvén所写的那样完成的：

 SELECT makegrid(the_geom, 1000) from mytable;

，或者如果您希望每行一个点：

 SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

希望对某人有用。
Alex

由于SRID错误，接受的答案不适用于我的数据。此修改效果更好。

– Vitaly Isaev
14-10-20在19:47

您可能要在多边形经过反子午线时添加一些东西吗？我可以想象会导致xmin / xmax出现问题。

–托马斯
16年7月8日在10:23

它对我不起作用。使用Postgres 9.6和PostGIS 2.3.3。在generate_series调用中，我必须将其作为第二个参数“ ceiling（st_xmax（$ 1））:: int”而不是“ ceiling（st_xmax（$ 1）-st_xmin（$ 1））:: int”和“ ceiling（ st_ymax（$ 1））:: int”而不是“ ceiling（st_ymax（$ 1）-st_ymin（$ 1））:: int”

–Vitor Sapucaia
18年3月27日在15:09

我赞成先前的评论； generate_series的上限应该是最大值的上限，而不是差异的上限（max-min）。

– R. Bourgeon
18年6月11日在16:18

@Alexandre Neto-根据Vitor的上述评论，可能值得更新您的答案

– dmci
20-10-20在21:10

#3 楼

使用wgs84几何图形的人可能会在使用此功能时遇到麻烦，因为

generate_series(floor(st_xmin())::int, ceiling(st_xmax())::int,) as x
,generate_series(floor(st_ymin())::int, ceiling(st_ymax())::int,) as y

仅返回整数。除了非常大的几何图形（例如国家/地区）（位于多个纬度和lng度上）之外，这将导致仅收集1个点，这在大多数情况下甚至不与几何图形本身相交... =>空结果！

我的麻烦是我似乎无法在浮点数（例如WSG84）上使用带有十进制距离的generate_series（），这就是为什么我要对该函数进行调整以使其正常工作的原因：

SELECT ST_Collect(st_setsrid(ST_POINT(x/1000000::float,y/1000000::float),st_srid())) FROM 
  generate_series(floor(st_xmin()*1000000)::int, ceiling(st_xmax()*1000000)::int,) as x ,
  generate_series(floor(st_ymin()*1000000)::int, ceiling(st_ymax()*1000000)::int,) as y 
WHERE st_intersects(,ST_SetSRID(ST_POINT(x/1000000::float,y/1000000::float),ST_SRID()))

基本相同。只需乘以1000000，即可在需要时将小数点加入游戏。

肯定有更好的解决方案。
++

这是一个聪明的解决方法。你检查结果了吗？他们一致吗？

–巴勃罗
13年5月29日在13:57

你好是的，巴勃罗。到目前为止，我对结果感到满意。我需要它来建立一些相对高度高于地面的多边形。（我使用SRTM计算每个网格点所需的高度）。我现在只缺少一种方法来包括放置在多边形外围的点。目前，渲染的形状在边缘处已被截断。

–朱利安·加西亚（Julien Garcia）
2013年6月12日11:34

工作，其他所有人的解决方案都失败了，谢谢！

–乔丹·阿塞诺（Jordan Arseno）
2014年8月7日在17:38

#4 楼

此算法应该很好：

 createGridInPolygon(polygon, resolution) {
    for(x=polygon.xmin; x<polygon.xmax; x+=resolution) {
       for(y=polygon.ymin; y<polygon.ymax; y+=resolution) {
          if(polygon.contains(x,y)) createPoint(x,y);
       }
    }
}

其中“多边形”是多边形和“分辨率” '是必需的网格分辨率。

要在PostGIS中实现它，可能需要以下功能：

ST_XMin，ST_XMax，ST_YMin和ST_YMax获取多边形的最小和最大坐标，

ST_包含测试多边形是否包含该点，
和ST_Point创建一个点。

好好运！

请注意，如果您有较大的复杂多边形区域（例如，我有海岸线缓冲区），则此方法并不是最佳选择。

– JaakL
2012年10月31日在8:34

那你有什么建议呢？

–朱利安
2012年10月31日在8:56

#5 楼

三种使用不同方法的算法。
Github Repo：https://github.com/I1mran/Postgis-Custom

最简单，最好的方法，使用从x和x到坐标的实际地球距离。 y方向。
该算法可与任何SRID配合使用，在内部可与WGS 1984（EPSG：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;
结果== ================================================== =====================

#6 楼

当然，这是另一种更快，更容易理解的方法。例如，对于1000m x 1000m的网格：

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0))).geom FROM the_polygon

SRID被保留。

此代码段将多边形的几何图形转换为空栅格，然后将每个像素转换为一个点。优点：我们不必再次检查原始多边形是否与这些点相交。

可选：

您还可以使用参数gridx和gridy添加网格对齐方式。但是由于我们使用每个像素的质心（而不是拐角），因此我们需要使用模来计算正确的值：

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0,mod(1000/2,100),mod(1000/2,100)))).geom FROM the_polygon

这是postgres函数：

CREATE OR REPLACE FUNCTION st_makegrid(geometry, float, integer) RETURNS SETOF geometry AS 'SELECT (ST_PixelAsCentroids(ST_AsRaster(,::float,::float,mod(::int/2,),mod(::int/2,)))).geom' LANGUAGE sql;

可用于：

SELECT st_makegrid(the_geom,1000.0,100) as geom from the_polygon -- makegrid(the_geom,grid_size,alignement)

评论
从最初的测试来看，这种方法在更复杂的多边形几何上似乎很慢。即使使用1000m的网格间距...

–Theo F
20年4月2日在14:13
@TheoF“相当慢”也很主观。应该将其与其他方法进行比较。在我的中档计算机上，我可以在不到6秒的时间内获得〜171000点。

–obchardon
20-4-7在17:04

#7 楼

所以我的固定版本是：

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer) RETURNS geometry AS 'SELECT ST_Collect(st_setsrid(ST_POINT(x,y),)) FROM generate_series(floor(st_xmin())::int, ceiling(st_xmax())::int,) as x ,generate_series(floor(st_ymin())::int, ceiling(st_ymax())::int,) as y where st_intersects(,st_setsrid(ST_POINT(x,y),))' LANGUAGE sql

用法：

SELECT (ST_Dump(makegrid(the_geom, 1000, 3857))).geom as the_geom from my_polygon_table

评论
嗨，我用makegrid函数得到空结果。使用shp2pgsql将shapefile导入到PostGIS。不知道是什么原因引起的麻烦，将srs设置为wgs84。

– Michal Zimmermann
13年3月15日在12:15

#8 楼

对先前答案的一个较小的潜在更新-第三个参数用作wgs84的比例（或对普通参数使用1），并在代码内部四舍五入，以使多个形状上的比例点对齐。

希望这会有所帮助，马丁

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer) RETURNS geometry AS /*geometry column , integer: distance between points, integer: scale factor for distance (useful for wgs84, e.g. use there 50000 as distance and 1000000 as scale factor*/ ' SELECT ST_Collect(st_setsrid(ST_POINT(x/::float,y/::float),st_srid())) FROM generate_series( (round(floor(st_xmin()*)::int/)*)::int, (round(ceiling(st_xmax()*)::int/)*)::int, ) as x , generate_series( (round(floor(st_ymin()*)::int/)*)::int, (round(ceiling(st_ymax()*)::int/)*)::int, ) as y WHERE st_intersects(,ST_SetSRID(ST_POINT(x/::float,y/::float),ST_SRID())) ' LANGUAGE sql

评论
将几何图形转换为特定的SRID（例如EPSG：3857）会比仅乘以比例因子更好吗？

– Nikolaus Krismer
17年4月25日在10:23

#9 楼

首先创建函数：

CREATE OR REPLACE FUNCTION st_polygrid(geometry, integer) RETURNS geometry AS $$ SELECT ST_SetSRID(ST_Collect(ST_POINT(x,y)), ST_SRID()) FROM generate_series(floor(ST_XMin())::numeric, ceiling(ST_xmax())::numeric, ) as x, generate_series(floor(ST_ymin())::numeric, ceiling(ST_ymax())::numeric,) as y WHERE ST_Intersects(,ST_SetSRID(ST_POINT(x,y), ST_SRID())) $$ LANGUAGE sql VOLATILE;

然后使用该函数创建点网格：

CREATE TABLE some.other_table AS SELECT a.gid, st_polygrid(a.the_geom, 1000) AS the_geom FROM some.table a

更改'1000'到所需的点间距。

#10 楼

我不需要功能，只需查询即可运行。根据这里的一些答案，我做了以下操作：

WITH geom AS ( SELECT st_setsrid('polygon-string'::geometry, 4326) AS geometry ), points AS ( SELECT generate_series(ST_XMIN(ST_Transform(geometry, 3857))::int, ST_XMAX(ST_Transform(geometry, 3857))::int, 1000) AS x, generate_series(ST_YMIN(ST_Transform(geometry, 3857))::int, ST_YMAX(ST_Transform(geometry, 3857))::int, 1000) AS y FROM geom ) SELECT ST_Collect(ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326)) AS points FROM points, geom WHERE st_intersects(geometry, ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326));

geom CTE将字符串转换为几何体并确保其SRID为4326
points CTE将几何体格式化为3857，并创建两个相距1000m的点列表
最后的查询将这些点转换回4326，过滤掉不在多边形的形状中，并将它们收集到几何图形集合中

编程黑洞网