# Creating regular polygon grid in PostGIS?

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 How to create regular point grid inside a 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.

• The generalization you seek is not clear. Are you saying you start with an (arbitrary) single polygon and wish to tile the plane with congruent copies of it? In general this is not possible, but maybe this polygon has particular properties (perhaps it is known to be a parallelogram, triangle, or hexagon, for instance). – whuber Nov 1 '11 at 1:25

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).

• Thanks! Now I just have to bind the fishnet to the BBox of the polygon. – mk.archaeo Nov 1 '11 at 12:29
• This is very helpful.. I have one query. How can I create grids inside a polygon/bbox? – Mohammed shafeek Feb 12 '14 at 5:18
• Nice Work Mike ,This is very helpful. – Mounaim Jan 20 '15 at 15:16

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.

• I wish I could up-vote this again... this was a perfect solution! and the ability to customize the coordinate system is fantastic~! – DPSSpatial Mar 4 '16 at 16:04
• Just a minor suggestion, instead of using `ST_GeomFromText` when creating a box to add to `sectors`, you can use `ST_MakeEnvelope` and just specify the bottom left and top right coordinates of the box. – Matt Oct 4 '17 at 15:13
• This brings potentials – nickves Jan 16 at 9:24

You can create a regular grid simply be 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
``````
• That is such a simple solution, having done it the vector way so many times. – John Powell Feb 13 '18 at 12:36

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);

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);
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;
``````

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;

``````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);
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
``````