Vjylku

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

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:

• 
  ST_XMin, ST_XMax, ST_YMin and ST_YMax to get the minimum and maximum coordinates of the polygon,


• 
  ST_Contains to test if the polygon contains the point,


• and ST_Point to create a point.

Good luck!

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

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

A small potential update to the previous answers - third argument as scale for wgs84 (or use 1 for normal ones), and also rounding inside the code so that the scaled points on multiple shapes are aligned.

Hope this helps, Martin

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)

RETURNS geometry AS







/*geometry column , integer: distance between points, integer: scale factor for distance (useful for wgs84, e.g. use there 50000 as distance and 1000000 as scale factor*/



'

SELECT ST_Collect(st_setsrid(ST_POINT(x/$3::float,y/$3::float),st_srid($1))) FROM 

  generate_series(

                (round(floor(st_xmin($1)*$3)::int/$2)*$2)::int, 

                (round(ceiling(st_xmax($1)*$3)::int/$2)*$2)::int,

                $2) as x ,

  generate_series(

                (round(floor(st_ymin($1)*$3)::int/$2)*$2)::int, 

                (round(ceiling(st_ymax($1)*$3)::int/$2)*$2)::int,

                $2) as y 

WHERE st_intersects($1,ST_SetSRID(ST_POINT(x/$3::float,y/$3::float),ST_SRID($1)))

'



LANGUAGE sql

Here is another approach which is certainly faster and easier to understand.

For example for a 1000m by 1000m grid:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0))).geom 

FROM the_polygon

Also the original SRID is preserved.

This snippet convert the polygon's geometry to an empty raster, then convert each pixel into a point. Advantage: we do not have to check again if the original polygon intersects the points.

Optional:

You can also add the grid alignement with the parameter gridx and gridy. But since we use the centroid of each pixel (and not a corner) we need to use a modulo to compute the right value:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0,mod(1000/2,100),mod(1000/2,100)))).geom 

FROM the_polygon

With mod(grid_size::integer/2,grid_precision)

Here is the postgres function:

CREATE OR REPLACE FUNCTION st_makegrid(geometry, float, integer)

RETURNS SETOF geometry AS

'SELECT (ST_PixelAsCentroids(ST_AsRaster($1,$2::float,$2::float,mod($2::int/2,$3),mod($2::int/2,$3)))).geom'

LANGUAGE sql;

Canbe used with:

SELECT makegrid(the_geom,1000.0,100) as geom from the_polygon  

-- makegrid(the_geom,grid_size,alignement)

Three algorithms using different approaches.

Git Repo Link

1. Here is a simple and best approach, using the actual distance of coordinates from x and y direction.
   The internal algorithm use the WGS 1984 (4326) and result transform
   to inserted SRID.

Function ===================================================================

CREATE OR REPLACE FUNCTION public.I_Grid_Point_Distance(geom public.geometry, x_side decimal, y_side decimal)

RETURNS public.geometry AS $BODY$

DECLARE

x_min decimal;

x_max decimal;

y_max decimal;

x decimal;

y decimal;

returnGeom public.geometry[];

i integer := -1;

srid integer := 4326;

input_srid integer;

BEGIN

CASE st_srid(geom) WHEN 0 THEN

    geom := ST_SetSRID(geom, srid);

        ----RAISE NOTICE 'No SRID Found.';

    ELSE

        ----RAISE NOTICE 'SRID Found.';

END CASE;

    input_srid:=st_srid(geom);

    geom := st_transform(geom, srid);

    x_min := ST_XMin(geom);

    x_max := ST_XMax(geom);

    y_max := ST_YMax(geom);

    y := ST_YMin(geom);

    x := x_min;

    i := i + 1;

    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);

<<yloop>>

LOOP

IF (y > y_max) THEN

    EXIT;

END IF;



CASE i WHEN 0 THEN 

    y := ST_Y(returnGeom[0]);

ELSE 

    y := ST_Y(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), y_side, radians(0))::geometry);

END CASE;



x := x_min;

<<xloop>>

LOOP

  IF (x > x_max) THEN

      EXIT;

  END IF;

    i := i + 1;

    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);

    x := ST_X(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), x_side, radians(90))::geometry);

END LOOP xloop;

END LOOP yloop;

RETURN

ST_CollectionExtract(st_transform(ST_Intersection(st_collect(returnGeom), geom), input_srid), 1);

END;

$BODY$ LANGUAGE plpgsql IMMUTABLE;

Use the function with a simple query, geometry must be valid and polygon, multi-polygons or envelope

SELECT I_Grid_Point_Distance(geom, 50, 61) from polygons limit 1;

Result======================================================================

2. 
   

   Second function based on Nicklas Avén algorithm. Have tried to handle any SRID.

   
   
   

   I have applied following changes in the algorithm.

   
   
   
   
   
   1. Separate variable for x and y direction for pixel size,
   
   
   2. New variable for calculates the distance in spheroid or ellipsoid.
   
   
   3. Input any SRID, function transform Geom to the working environment of
      Spheroid or Ellipsoid Datum, then apply the distance to each side, get
      the result and transform to input SRID.
   
   
   
   

Function ===================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)

RETURNS SETOF geometry AS $BODY$ 

DECLARE

x_max decimal; 

y_max decimal;

x_min decimal;

y_min decimal;

srid integer := 4326;

input_srid integer; 

BEGIN

CASE st_srid(geom) WHEN 0 THEN

  geom := ST_SetSRID(geom, srid);

  RAISE NOTICE 'SRID Not Found.';

    ELSE

        RAISE NOTICE 'SRID Found.';

    END CASE;



    CASE spheroid WHEN false THEN

        RAISE NOTICE 'Spheroid False';

        srid := 4326;

        x_side := x_side / 100000;

        y_side := y_side / 100000;

    else

        srid := 900913;

        RAISE NOTICE 'Spheroid True';

    END CASE;

    input_srid:=st_srid(geom);

    geom := st_transform(geom, srid);

    x_max := ST_XMax(geom);

    y_max := ST_YMax(geom);

    x_min := ST_XMin(geom);

    y_min := ST_YMin(geom);

RETURN QUERY

WITH res as (SELECT ST_SetSRID(ST_MakePoint(x, y), srid) point FROM

generate_series(x_min, x_max, x_side) as x,

generate_series(y_min, y_max, y_side) as y

WHERE st_intersects(geom, ST_SetSRID(ST_MakePoint(x, y), srid))

) select ST_TRANSFORM(ST_COLLECT(point), input_srid) from res;

END;

$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Use it with a simple query.

SELECT I_Grid_Point(geom, 22, 15, false) from polygons;

Result===================================================================

3. Function based on series generator.

Function==================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)

RETURNS SETOF geometry AS $BODY$

DECLARE

x_max decimal;

y_max decimal;

x_min decimal;

y_min decimal;

srid integer := 4326;

input_srid integer;

x_series DECIMAL;

y_series DECIMAL;

BEGIN

CASE st_srid(geom) WHEN 0 THEN

  geom := ST_SetSRID(geom, srid);

  RAISE NOTICE 'SRID Not Found.';

    ELSE

        RAISE NOTICE 'SRID Found.';

    END CASE;



    CASE spheroid WHEN false THEN

        RAISE NOTICE 'Spheroid False';

    else

        srid := 900913;

        RAISE NOTICE 'Spheroid True';

    END CASE;

    input_srid:=st_srid(geom);

    geom := st_transform(geom, srid);

    x_max := ST_XMax(geom);

    y_max := ST_YMax(geom);

    x_min := ST_XMin(geom);

    y_min := ST_YMin(geom);



    x_series := CEIL ( @( x_max - x_min ) / x_side);

    y_series := CEIL ( @( y_max - y_min ) / y_side );

RETURN QUERY

SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM

generate_series(0, x_series) as x,

generate_series(0, y_series) as y

WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);

END;

$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Use it with a simple query.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons;
Result==========================================================================