Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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 Create regular point grid 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.

share|improve this question
1  
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

2 Answers 2

up vote 28 down vote accepted

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

share|improve this answer
    
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 at 5:18

Suggest some 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 liability in a region).

Function is not very elegant, but I didn't found any better solution for that task (inc. Mike Toews' function). 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 optimal for the west part of Russia
)
RETURNS public.geometry AS
$body$
DECLARE
  BoundM public.geometry; --Bound polygon transformed to 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 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 bound for one step. You always can crop the result. But if not you may get not quite correct data for outbound polygons (if you calculate frequency per a sector  e.g.)
        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 formated by generated polygons lines lie along the geographic parallels and meridians - that's very convenient.

An example of a grid with 50 km step

Hope it will be useful for somebody.

Advice: If you compute something like density (e.g. map of lightning strokes by cells), and the grid is generated dynamically, using temporary tables for storing cells as geometry polygons don't forget to create a spatial index for a field of cell. Performance will extremely increase!

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.