在PostGIS中创建常规多边形网格

如何在PostGIS中以多边形的形状创建给定大小的多边形/正方形的规则网格？
我考虑过类似在PostGIS中仅在正方形中在多边形内部创建规则点网格的功能，所以正方形可以是5m x 5m甚至10m x 10m。但是不知道以正确的方式更改它。

您寻求的概括尚不清楚。您是说要从一个（任意）单一多边形开始，并希望用其全等副本对平面进行平铺？通常这是不可能的，但是此多边形可能具有特殊的属性（例如，已知它是平行四边形，三角形或六边形）。

#1 楼

这是一个返回函数ST_CreateFishnet的集合，该函数创建一个二维多边形几何网格：单元格大小的长度以及可选的nrow和ncol是左下角的坐标。

结果是xsize和ysize的数字，从左下角的1开始，以及x0矩形多边形每个单元格因此，例如：

CREATE OR REPLACE FUNCTION ST_CreateFishnet(
        nrow integer, ncol integer,
        xsize float8, ysize float8,
        x0 float8 DEFAULT 0, y0 float8 DEFAULT 0,
        OUT "row" integer, OUT col integer,
        OUT geom geometry)
    RETURNS SETOF record AS
$$
SELECT i + 1 AS row, j + 1 AS col, ST_Translate(cell, j *  + , i *  + ) AS geom
FROM generate_series(0,  - 1) AS i,
     generate_series(0,  - 1) AS j,
(
SELECT ('POLYGON((0 0, 0 '||||', '||||' '||||', '||||' 0,0 0))')::geometry AS cell
) AS foo;
$$ LANGUAGE sql IMMUTABLE STRICT;

或为整个网格制作单个几何图形集：

SELECT *
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;
 row | col |         geom
-----+-----+--------------------------------
   1 |   1 | 0103000000010000000500000000...
   2 |   1 | 0103000000010000000500000000...
   3 |   1 | 0103000000010000000500000000...
   4 |   1 | 0103000000010000000500000000...
   1 |   2 | 0103000000010000000500000000...
   2 |   2 | 0103000000010000000500000000...
   ...
   3 |   6 | 0103000000010000000500000000...
   4 |   6 | 0103000000010000000500000000...
(24 rows)

您可以添加y0 / row原点偏移（默认为零）。

谢谢！现在，我只需要将渔网绑定到多边形的BBox。

–mk.archaeo
2011年11月1日12:29

这非常有帮助。我有一个查询。如何在多边形/ bbox中创建网格？

–穆罕默德·谢菲克（Mohammed shafeek）
2014年2月12日5:18

尼斯工作麦克，这非常有帮助。

– Mounaim
2015年1月20日15:16

#2 楼

这是生成的一种特定变体，适用于需要为地理地图创建具有恒定度量步长的网格的情况（可能使用像元对值进行分组，例如区域中的闪电密度）。

功能不是很好，但是我没有找到更好的解决方案（包括上面的Mike Toews函数）。
所以您有一个绑定的多边形（例如从Google Maps界面到达）以米为单位的阶跃值：

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  grid_step integer,
  metric_srid integer = 28408 --metric SRID (this particular is optimal for the Western Russia)
)
RETURNS public.geometry AS
$body$
DECLARE
  BoundM public.geometry; --Bound polygon transformed to the metric projection (with metric_srid SRID)
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  sectors public.geometry[];
  i INTEGER;
BEGIN
  BoundM := ST_Transform(, ); --From WGS84 (SRID 4326) to the metric projection, to operate with step in meters
  Xmin := ST_XMin(BoundM);
  Xmax := ST_XMax(BoundM);
  Ymax := ST_YMax(BoundM);

  Y := ST_YMin(BoundM); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  --Better if generating polygons exceeds the bound for one step. You always can crop the result. But if not you may get not quite correct data for outbound polygons (e.g. if you calculate frequency per sector)
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      i := i + 1;
      sectors[i] := ST_GeomFromText('POLYGON(('||X||' '||Y||', '||(X+)||' '||Y||', '||(X+)||' '||(Y+)||', '||X||' '||(Y+)||', '||X||' '||Y||'))', );

      X := X + ;
    END LOOP xloop;
    Y := Y + ;
  END LOOP yloop;

  RETURN ST_Transform(ST_Collect(sectors), ST_SRID());
END;
$body$
LANGUAGE 'plpgsql';

使用方法：

SELECT cell FROM 
(SELECT (
ST_Dump(makegrid_2d(ST_GeomFromText('Polygon((35.099577 45.183417,47.283415 45.183417,47.283415 49.640445,35.099577 49.640445,35.099577 45.183417))',
 4326), -- WGS84 SRID
 10000) -- cell step in meters
)).geom AS cell) AS q_grid

所以您可以看到由生成的多边形格式化的线沿地理平行线和子午线分布-非常方便。并动态生成网格。为提高性能，我建议使用临时表将单元格存储为几何多边形，该列上的空间索引表示ce ll。

我希望我能再次投票赞成……这是一个完美的解决方案！而且自定义坐标系的能力也很棒〜！

– DPSSpatial
16-3-4在16:04

这只是一个小建议，您可以使用ST_MakeEnvelope并指定框的左下角和右上角坐标，而不用在创建要添加到扇区的框时使用ST_GeomFromText。

–马特
17-10-4在15:13

这带来了潜力

–尼克斯
19年1月16日在9:24

#3 楼

您只需将空栅格矢量化即可创建常规网格：

SELECT (ST_PixelAsPolygons(ST_AddBand(ST_MakeEmptyRaster(100, 100, 1.1, 1.1, 1.0), '8BSI'::text, 1, 0), 1, false)).geom

这是一个非常简单的解决方案，它多次执行矢量操作。

–约翰·鲍威尔（John Powell）
18年2月13日在12:36

在这种情况下，upperleftx和upperlefty值有什么影响？我需要创建一个700000m x 1400000m的网格，每个“像素”为1000x1000m，但是我不确定为这些值设置什么。

– DeadPassive
6月11日10:52

#4 楼

我创建了@Alexander函数的变体，不需要我们转换为另一个SRID。这样避免了必须找到一个以米为特定区域单位的投影的问题。它使用ST_Project使用给定的投影正确跨步。我还添加了width_step和height_step，以允许使用矩形瓷砖，而不是要求它们是正方形。 br />

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  width_step integer,
  height_step integer
)
RETURNS public.geometry AS
$body$
DECLARE
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  NextX DOUBLE PRECISION;
  NextY DOUBLE PRECISION;
  CPoint public.geometry;
  sectors public.geometry[];
  i INTEGER;
  SRID INTEGER;
BEGIN
  Xmin := ST_XMin(bound_polygon);
  Xmax := ST_XMax(bound_polygon);
  Ymax := ST_YMax(bound_polygon);
  SRID := ST_SRID(bound_polygon);

  Y := ST_YMin(bound_polygon); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
      NextX := ST_X(ST_Project(CPoint, , radians(90))::geometry);
      NextY := ST_Y(ST_Project(CPoint, , radians(0))::geometry);

      i := i + 1;
      sectors[i] := ST_MakeEnvelope(X, Y, NextX, NextY, SRID);

      X := NextX;
    END LOOP xloop;
    CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
    NextY := ST_Y(ST_Project(CPoint, , radians(0))::geometry);
    Y := NextY;
  END LOOP yloop;

  RETURN ST_Collect(sectors);
END;
$body$
LANGUAGE 'plpgsql';

#5 楼

这是一种优化高效的算法，可在任何信封，多边形或Multipolygons内创建渔网，规则网格，多边形网格，矩形网格。
几乎可以处理任何SRID;
GitHub Repo：https：// github。 com / I1mran / Postgis-Custom

     DROP FUNCTION IF EXISTS PUBLIC.I_Grid_Regular(geometry, float8, float8);
    CREATE OR REPLACE FUNCTION PUBLIC.I_Grid_Regular
    ( geom geometry, x_side float8, y_side float8, OUT geometry )
    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;
    geom_cell geometry := ST_GeomFromText(FORMAT('POLYGON((0 0, 0 %s, %s %s, %s 0,0 0))',
                                            , , , ), srid);
    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;
    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 With foo AS (
        SELECT
        ST_Translate( geom_cell, j *  + x_min, i *  + y_min ) AS cell
        FROM
            generate_series ( 0, x_series ) AS j,
            generate_series ( 0, y_series ) AS i
        ) SELECT ST_CollectionExtract(ST_Collect(ST_Transform ( ST_Intersection(cell, geom), input_srid)), 3)
        FROM foo where ST_intersects (cell, geom);
    END;
    $BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

通过简单查询使用它；
输入应该是有效的多边形，多面体或信封。

select I_Grid_Regular(st_setsrid(g.geom, 4326), .0001, .0001 ), geom from polygons limit 1

#6 楼

我创建了以下函数，用于根据多边形表在PostGIS中创建网格。网格是矩形，可以在一组给定的多边形内。代码首先创建一个栅格，然后将其像素转储为多边形。
要创建100,000个网格像元，需要0.5s。
如下：

schem：输入和输出表所在的架构。
in_table_name：包含我们要拥有的区域
的输入表网格。该表可以具有
许多行，将为此表的所有
行创建输出网格。。
out_table_name：将在相同架构中创建的输出表，并将具有w个输出网格行。
out_table_name：将创建为具有网格的输出表的名称。
grid_size：创建网格的矩形的大小。它应该以输出的斜率（您将此功能作为最后一个输入赋予的斜率）的度量单位。
srid：输出网格的所需斜率。实际上，输入表的分隔并不重要。您不必担心。

因此，将网格的首选边界存储在一个表中，让我们说“ grid_boundary”，然后说几何列称为“ geom”。另外，可以说输入和输出表都将处于公共模式中。如果要在SRID：3979中使用1000m * 1000m的网格，则将函数调用为：任何东西。 2-您的网格大小应与用作第6个输入的SRID的度量单位相同。 3-您可以在输入表中有许多行多边形作为输出网格的边界

它的工作原理如下：

编程黑洞网