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, inside a polygon, a regular grid of point spaced x,y in postgis? Like the example:

alt text

share|improve this question
    
I tried to make clipping polygons merging this code with "postGIS in action" dicing code but only one polygon is created... Did I forget something? CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer) RETURNS geometry AS 'SELECT st_intersection(g1.geom1, g2.geom2) AS geom FROM (SELECT $1 AS geom1) AS g1 INNER JOIN (Select st_setsrid(CAST(ST_MakeBox2d( st_setsrid(ST_Point(x,y),$3), st_setsrid(ST_Point(x + $2 ,y +$2),$3)) as geometry),$3) as geom2 FROM generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x, generate_series(floor(st_ymin($1))::int, ceiling(st_ymax( –  aurel_nc May 2 '13 at 14:30

5 Answers 5

up vote 18 down vote accepted

You do that with generate_series.

If you don't want to manually write where the grid is to start and stop, the easiest is to create a function.

I have not tested the below properly, but I think it should work:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_POINT(x,y)) FROM 
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y 
where st_intersects($1,ST_POINT(x,y))'
LANGUAGE sql

To use it you can do:

SELECT makegrid(the_geom, 1000) from mytable;

where the first argument is the polygon you want the grid in, and the second argument is the distance between the points in the grid.

If you want one point per row you just use ST_Dump like:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

HTH

Nicklas

share|improve this answer
    
You may need to add st_setSRID() to the st_point functions, otherwise st_intersects does not work. –  JaakL Oct 31 '12 at 9:25
    
Added my tested version as separate answer. –  JaakL Oct 31 '12 at 16:43

This algorithm should be fine:

createGridInPolygon(polygon, resolution) {
    for(x=polygon.xmin; x<polygon.xmax; x+=resolution) {
       for(y=polygon.ymin; y<polygon.ymax; y+=resolution) {
          if(polygon.contains(x,y)) createPoint(x,y);
       }
    }
}

where 'polygon' is the polygon and 'resolution' is the required grid resolution.

To implement it in PostGIS, the following functions may be needed:

Good luck!

share|improve this answer
    
Note that if you have large complex polygons areas (e.g. I have coastline buffer), then this approach is not quite optimal. –  JaakL Oct 31 '12 at 8:34
    
Then, what do you propose instead? –  julien Oct 31 '12 at 8:56

I have picked up Nicklas Avén makegrid function code and made it a bit more generic by reading and using the srid from the polygon geometry. Otherwise using a polygon with a defined srid, would give an error.

The function:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x,y),ST_SRID($1))) FROM 
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y 
where st_intersects($1,ST_SetSRID(ST_POINT(x,y),ST_SRID($1)))'
LANGUAGE sql

To use the function is done exactly as Nicklas Avén wrote:

SELECT makegrid(the_geom, 1000) from mytable;

or if you want one point per row:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

Hope this will be usefull for someone.

Alex

share|improve this answer

People using a wgs84 geometry will probably have trouble with this function since

generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y 

only return integers. Except for very big geometries such as countries (that are laying on multiple lat, lng degrees), this will cause to collect only 1 point which is most of the time not even intersecting the geometry itself... => empty result !

My trouble was I can not seem to use generate_series() with a decimal distance on floating numbers such as those WSG84... This is why I tweaked the function to get it working anyway :

SELECT ST_Collect(st_setsrid(ST_POINT(x/1000000::float,y/1000000::float),st_srid($1))) FROM 
  generate_series(floor(st_xmin($1)*1000000)::int, ceiling(st_xmax($1)*1000000)::int,$2) as x ,
  generate_series(floor(st_ymin($1)*1000000)::int, ceiling(st_ymax($1)*1000000)::int,$2) as y 
WHERE st_intersects($1,ST_SetSRID(ST_POINT(x/1000000::float,y/1000000::float),ST_SRID($1)))

Basically exactly the same. Just multiplying and dividing by 1000000 to get the decimals into the game when I need it.

There is surely a better solution to achieve that. ++

share|improve this answer
    
That's a smart workaround. Have you checked the results? Are they consistent? –  Pablo May 29 '13 at 13:57
    
Hi. Yes pablo. I'm happy with the results so far. I needed it to build some polygon with relative altitude above the ground. (I use SRTM to calculate the altitude I want for every grid point). I'm now only missing a way to include the points that are laying on the perimeter of the polygon as well. Currently the rendered shape is somewhat truncated at the edge. –  Julien Garcia Jun 12 '13 at 11:34

So my fixed version:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)
RETURNS geometry AS
'SELECT ST_Collect(st_setsrid(ST_POINT(x,y),$3)) FROM 
  generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
  ,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y 
where st_intersects($1,st_setsrid(ST_POINT(x,y),$3))'
LANGUAGE sql

Usage:

SELECT (ST_Dump(makegrid(the_geom, 1000, 3857))).geom as the_geom from my_polygon_table
share|improve this answer
    
Hi, I get empty result with makegrid function. The shapefile was imported to PostGIS using shp2pgsql. Have no idea what could be causing trouble, the srs is set to wgs84. –  zimmi Mar 15 '13 at 12:15

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.