Creating regular polygon grid in PostGIS

Question

How to create, on the shape of a polygon, a regular grid of polygons/squares of an given size, in PostGIS?

I have thought about a function like Creating regular point grid inside polygon in PostGIS only for squares, so that squares can be 5m x 5m or even 10m x 10m. But have no idea to change this in the correct way.

Answer 1

Here is a set returning function ST_CreateFishnet that creates a 2D grid of polygon geometries:

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 * $3 + $5, i * $4 + $6) AS geom
FROM generate_series(0, $1 - 1) AS i,
     generate_series(0, $2 - 1) AS j,
(
SELECT ('POLYGON((0 0, 0 '||$4||', '||$3||' '||$4||', '||$3||' 0,0 0))')::geometry AS cell
) AS foo;
$$ LANGUAGE sql IMMUTABLE STRICT;

where nrow and ncol are the number of rows and columns, xsize and ysize are the lengths of the cell size, and optional x0 and y0 are coordinates for the bottom-left corner.

The result is row and col numbers, starting from 1 at the bottom-left corner, and geom rectangular polygons for each cell. So for example:

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)

Or to make a single geometry collection for the full grid:

SELECT ST_Collect(cells.geom)
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;

You can add the x0 / y0 origin offsets (these defaulted to zero).

Answer 2

Here's a specific variant of generating, for a situation when you need to create a grid for a geographic map with a constant metric step (cells might be used for grouping values, e.g. lightning density in a region).

Function is not very elegant, but I didn't find any better solution for that task (incl. the Mike Toews' function above). So you have a bound polygon (e.g. arrived from a Google Maps interface), have a step value in meters:

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($1, $3); --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+$2)||' '||Y||', '||(X+$2)||' '||(Y+$2)||', '||X||' '||(Y+$2)||', '||X||' '||Y||'))', $3);

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

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

How to use it:

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

So you can see that lines formatted by generated polygons lie along the geographic parallels and meridians - that's very convenient.

Advice: If you compute something like density (e.g. map of lightning strokes by cells), and the grid is generated dynamically To increase the performance I would suggest using temporary tables for storing cells as geometry polygons, with a spatial index on the the column represents the cell.

Answer 3

You can create a regular grid simply by vectorising an empty raster:

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

Answer 4

I've created a variant of @Alexander's function that doesn't require that we transform to another SRID. This avoids the issue of having to find a projection that uses meters as units for a particular region. It uses ST_Project to properly stride using the given projection. I've also added a width_step and height_step to allow for rectangular tiles instead of requiring they be squares.

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, $2, radians(90))::geometry);
      NextY := ST_Y(ST_Project(CPoint, $3, 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, $3, radians(0))::geometry);
    Y := NextY;
  END LOOP yloop;

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

You can use it like so:

SELECT ST_AsGeoJSON(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
        ),
         10000, -- width step in meters
         10000  -- height step in meters
       ) 
    )
  ) .geom AS cell
)q;

Answer 5

Here is an optimized and efficient algorithm to create fishnet, regular grid, polygon grid, rectangular grid inside any envelope, polygon, or multipolygons. almost handle any SRID;

GitHub Repo: https://github.com/imran-5/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))',
                                            $3, $2, $3, $2), 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 * $2 + x_min, i * $3 + 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;

Use it with a simple query; input should be a valid polygon, Multipolygon, or envelope.

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

Answer 6

I created the following function for creating a grid in PostGIS from a table of polygons. The grid is rectangular and it can be inside a given set of polygons. The code first creates a raster then dump its pixels as polygons. For creating 100,000 grid cells, it takes 0.5s.

DROP FUNCTION IF EXISTS  CreateGrid;

CREATE OR REPLACE FUNCTION CreateGrid(schem TEXT, in_table_name 
TEXT,geom_column TEXT,
                                  out_table_name TEXT, grid_size INTEGER , 
srid INTEGER) 
returns void
AS $$
BEGIN
EXECUTE  'CREATE TABLE ' || schem || '.' || out_table_name ||' as(  SELECT 
(ST_PixelAsPolygons(rast,1)).geom AS geom
               FROM (SELECT ST_AsRaster(ST_Transform( ' || ' ST_Union( ' || 
geom_column  || ')' || ' , ' || srid|| 
               '), ' || 'CAST(' || grid_size || ' AS DOUBLE PRECISION) ' || 
', '
               ||  'CAST(' || grid_size || ' AS DOUBLE PRECISION) ' ||', 
''4BUI'', 1, 0) AS rast
               FROM ' || schem || '.' || in_table_name ||
               ') as foo   )';
END; $$
LANGUAGE plpgsql;

The input arguments are as follows:

1. schem: the schema that the input and output tables are in it.

2. in_table_name: The input table that contains the area that we want to have the grid. this table can have many rows the output grid will be created for all rows of this table.

3. geom_column: The name of the column in input table ("in_table_name") that contains the geometry. out_table_name: The output table that will be created in the same schema and will have thew output grid rows.

4. out_table_name: The name of the output table that will be created to have the grid. this will be created in the same schema.

5. grid_size: The size of the rectangles that create the grid. It should be in the unit of measurement of the srid of output(the srid that you give to this function as last input).

6. srid: The desired srid of the output grid. In fact, the srid of the input table does not matter. You do not need to worry about it.

Therefore, store the preferred boundary of your grid in a table, lets say "grid_boundary" and lets say that the geometry column is call "geom". In addition, lets say that both input and output tables will be in public schema. If you want a grid of 1000m*1000m in SRID:3979, then call the function as:

SELECT CreateGrid('public', 'grid_boundary', 'geom', 
'out_grid', 1000 , 3979);

Note: 1- your input srid can be anything. 2- your grid size should be the same as the measurement unit of your SRID that you give to function as 6th input. 3- you can have many rows of polygon in input table as the boundaries of your output grid

It works like this:

Answer 7

I improved one of the previous answers which was limited by Postgres max array length. This solution generates a grid based on metric extent, metric step and with a chosen CRS. At the same time it will work efficiently for bigger grids. I think it's quite short and clean. :

CREATE OR replace FUNCTION generate_grid(bound_polygon geometry, grid_step integer, srid integer default 2180)
 RETURNS table(id bigint, geom geometry)
 LANGUAGE plpgsql
AS $function$
DECLARE
  Xmin int;
  Xmax int;
  Ymax int;
  Ymin int;
  query_text text;
BEGIN
  Xmin := floor(ST_XMin(bound_polygon));
  Xmax := ceil(ST_XMax(bound_polygon));
  Ymin := floor(ST_YMin(bound_polygon));
  Ymax := ceil(ST_YMax(bound_polygon));
    
  query_text := 'select row_number() over() id, st_makeenvelope(s1, s2, s1+$5, s2+$5, $6) geom
    from generate_series($1, $2+$5, $5) s1, generate_series ($3, $4+$5, $5) s2';

  RETURN QUERY EXECUTE query_text using Xmin, Xmax, Ymin, Ymax, grid_step, srid;
END;
$function$
;

You can later use it like that:

create table grid as select id, geom from generate_grid(geom, step, crs_id);