Monday, 1 January 2018

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.



Answer



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;


4x6 grid


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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...