Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Generating a Grid (fishnet) of points or polygons for PostGIS

Tuesday September 11 2012 at 01:55

Keywordsgrid fishnet points polygons squares postgis
Summary

An article showing how to generate a grid of points or polygons for a geometry object in PostGIS.

I wrote an article on Gridding a PostGIS Geometry object a while back.

I had cause, recently, to use a variation on this for the generation of an array of point objects for a geometry polygon, so I wrote a version that also allows for the generation of a grid of polygons or points depending on a boolean parameter.

  1. -- Create required type
  2. DROP   TYPE IF EXISTS T_Grid CASCADE;
  3. CREATE TYPE T_Grid AS (
  4.   gcol  int4,
  5.   grow  int4,
  6.   geom geometry
  7. );
  8. -- Drop function is exists
  9. DROP FUNCTION IF EXISTS ST_RegularGrid(geometry, NUMERIC, NUMERIC, BOOLEAN);
  10. -- Now create the function
  11. CREATE OR REPLACE FUNCTION ST_RegularGrid(p_geometry   geometry,
  12.                                           p_TileSizeX  NUMERIC,
  13.                                           p_TileSizeY  NUMERIC,
  14.                                           p_point      BOOLEAN DEFAULT TRUE)
  15.   RETURNS SETOF T_Grid AS
  16. $BODY$
  17. DECLARE
  18.    v_mbr   geometry;
  19.    v_srid  int4;
  20.    v_halfX NUMERIC := p_TileSizeX / 2.0;
  21.    v_halfY NUMERIC := p_TileSizeY / 2.0;
  22.    v_loCol int4;
  23.    v_hiCol int4;
  24.    v_loRow int4;
  25.    v_hiRow int4;
  26.    v_grid  T_Grid;
  27. BEGIN
  28.    IF ( p_geometry IS NULL ) THEN
  29.       RETURN;
  30.    END IF;
  31.    v_srid  := ST_SRID(p_geometry);
  32.    v_mbr   := ST_Envelope(p_geometry);
  33.    v_loCol := trunc((ST_XMIN(v_mbr) / p_TileSizeX)::NUMERIC );
  34.    v_hiCol := CEIL( (ST_XMAX(v_mbr) / p_TileSizeX)::NUMERIC ) - 1;
  35.    v_loRow := trunc((ST_YMIN(v_mbr) / p_TileSizeY)::NUMERIC );
  36.    v_hiRow := CEIL( (ST_YMAX(v_mbr) / p_TileSizeY)::NUMERIC ) - 1;
  37.    FOR v_col IN v_loCol..v_hiCol Loop
  38.      FOR v_row IN v_loRow..v_hiRow Loop
  39.          v_grid.gcol := v_col;
  40.          v_grid.grow := v_row;
  41.          IF ( p_point ) THEN
  42.            v_grid.geom := ST_SetSRID(
  43.                              ST_MakePoint((v_col * p_TileSizeX) + v_halfX,
  44.                                           (v_row * p_TileSizeY) + V_HalfY),
  45.                              v_srid);
  46.          ELSE
  47.            v_grid.geom := ST_SetSRID(
  48.                              ST_MakeEnvelope((v_col * p_TileSizeX),
  49.                                              (v_row * p_TileSizeY),
  50.                                              (v_col * p_TileSizeX) + p_TileSizeX,
  51.                                              (v_row * p_TileSizeY) + p_TileSizeY),
  52.                              v_srid);
  53.          END IF;
  54.          RETURN NEXT v_grid;
  55.      END Loop;
  56.    END Loop;
  57. END;
  58. $BODY$
  59.   LANGUAGE plpgsql IMMUTABLE
  60.   COST 100
  61.   ROWS 1000;
  62. -- Assign ownership
  63. ALTER FUNCTION ss.st_regulargrid(geometry, NUMERIC, NUMERIC, BOOLEAN)
  64.   OWNER TO postgres;

Now, let’s do some testing:

First, generate grid of points over a geometry’s MBR.

  1. SELECT gcol, grow,ST_AsText(geom) AS geomWKT
  2.   FROM ss.ST_RegularGrid(ST_GeomFromText('LINESTRING(0 0, 100 100)',0),20,20);

Result:

gcol grow geomWKT
0 0 POINT (10 10)
0 1 POINT (10 30)
0 2 POINT (10 50)
0 3 POINT (10 70)
0 4 POINT (10 90)
1 0 POINT (30 10)
1 1 POINT (30 30)
1 2 POINT (30 50)
1 3 POINT (30 70)
1 4 POINT (30 90)
2 0 POINT (50 10)
2 1 POINT (50 30)
2 2 POINT (50 50)
2 3 POINT (50 70)
2 4 POINT (50 90)
3 0 POINT (70 10)
3 1 POINT (70 30)
3 2 POINT (70 50)
3 3 POINT (70 70)
3 4 POINT (70 90)
4 0 POINT (90 10)
4 1 POINT (90 30)
4 2 POINT (90 50)
4 3 POINT (90 70)
4 4 POINT (90 90)

Secondly, generate grid of polygons over the MBR of the same object.

  1. SELECT gcol, grow,ST_AsText(geom) AS geomWKT
  2.   FROM ss.ST_RegularGrid(ST_GeomFromText('LINESTRING(0 0, 100 100)',0),20,20,FALSE);

Result:

gcol grow geomWKT
0 0 POLYGON ((0 0,0 20,20 20,20 0,0 0))
0 1 POLYGON ((0 20,0 40,20 40,20 20,0 20))
0 2 POLYGON ((0 40,0 60,20 60,20 40,0 40))
0 3 POLYGON ((0 60,0 80,20 80,20 60,0 60))
0 4 POLYGON ((0 80,0 100,20 100,20 80,0 80))
1 0 POLYGON ((20 0,20 20,40 20,40 0,20 0))
1 1 POLYGON ((20 20,20 40,40 40,40 20,20 20))
1 2 POLYGON ((20 40,20 60,40 60,40 40,20 40))
1 3 POLYGON ((20 60,20 80,40 80,40 60,20 60))
1 4 POLYGON ((20 80,20 100,40 100,40 80,20 80))
2 0 POLYGON ((40 0,40 20,60 20,60 0,40 0))
2 1 POLYGON ((40 20,40 40,60 40,60 20,40 20))
2 2 POLYGON ((40 40,40 60,60 60,60 40,40 40))
2 3 POLYGON ((40 60,40 80,60 80,60 60,40 60))
2 4 POLYGON ((40 80,40 100,60 100,60 80,40 80))
3 0 POLYGON ((60 0,60 20,80 20,80 0,60 0))
3 1 POLYGON ((60 20,60 40,80 40,80 20,60 20))
3 2 POLYGON ((60 40,60 60,80 60,80 40,60 40))
3 3 POLYGON ((60 60,60 80,80 80,80 60,60 60))
3 4 POLYGON ((60 80,60 100,80 100,80 80,60 80))
4 0 POLYGON ((80 0,80 20,100 20,100 0,80 0))
4 1 POLYGON ((80 20,80 40,100 40,100 20,80 20))
4 2 POLYGON ((80 40,80 60,100 60,100 40,80 40))
4 3 POLYGON ((80 60,80 80,100 80,100 60,80 60))
4 4 POLYGON ((80 80,80 100,100 100,100 80,80 80))

Visually.

I hope this is useful for someone.

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment [2]

Hi,
A probably more efficient way of generating a grid (of point or polygons) would be to use the generate_series postgresql function.
You can find some similar use cases here :
http://gis.stackexchange.com/questions/4663/how-to-create-regular-point-grid-inside-a-polygon-in-postgis
or here :
http://www.bostongis.com/postgis_translate.snippet

The «Postgis in Action» book includes some great examples of this kind too.

— vincent · 11 September 2012, 13:42 · #

Vincent,

Thanks for your comments.

I am aware of these techniques and have used them myself beforehand.

Sometimes in the spirit of horses for courses I use one or the other.

This technique is useful because it encapsulates the processing inside a function cleaning up the SQL from extraneous detail that would otherwise complicate a SQL’s declarative power in solving a problem.

regards
Simon

— Simon Greener · 11 September 2012, 23:44 · #

Article Navigation:   Previous   Next