Creating regular polygon grid in PostGIS
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).
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.
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