Creating regular point grid inside polygon in PostGIS

Question

How to create, inside a polygon, a regular grid of point spaced x,y in PostGIS?

Like the example:

alt text

• aurel_nc: 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(
• user103524: look at my detailed answer.
• Peter Krauss: For an uniform sample-set grid you need an equal-area projection: see this answer below or standard ISO DGGS - Discrete Global Grid Systems.
• Peter Krauss: Best performance, for Big Data grids, see solution at gis.stackexchange.com/a/490156/7505

Answer 1

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;

• JaakL: You may need to add st_setSRID() to the st_point functions, otherwise st_intersects does not work.
• JaakL: Added my tested version as separate answer.
• Peter Krauss: something wrong and not make sense, "from min to interval" not runs. I need to correct series to "from min to max" generate_series(floor(ST_XMIN(1ドル))::int, ceiling(ST_XMAX(1ドル) )::int, 2ドル) AS x
• Peter Krauss: ... no answer here after 6 months, ok, see the solution below, suitable for equal-area projections.

Answer 2

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

• Vitaly Isaev: The accepted answer doesn't work with my data due to SRID Errors. This modification works better.
• Thomas: You might want to add something when polygon is crossing the antimeridian? I can imagine it would cause a problem with xmin/xmax.
• Vitor Sapucaia: It didn't work for me. Using Postgres 9.6 and PostGIS 2.3.3. Inside the generate_series call, I had to put this as the second parameter "ceiling(st_xmax(1ドル))::int" instead of "ceiling(st_xmax(1ドル)-st_xmin(1ドル))::int", and "ceiling(st_ymax(1ドル))::int" instead of "ceiling(st_ymax(1ドル)-st_ymin(1ドル))::int"
• R. Bourgeon: I approve the previous comment; the upper bound of generate_series should be the ceiling of the max, and not the ceiling of the difference (max - min).
• dmci: @Alexandre Neto - it might be worth updating your answer based on the comment above from Vitor

Answer 3

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

• Pablo: That's a smart workaround. Have you checked the results? Are they consistent?
• Julien Garcia: 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.
• Jordan Arsenault: worked, when all others' solutions failed, thanks!

Answer 4

Three algorithms using different methods.

Github Repo: https://github.com/muimsd/Postgis-Grids

1. The simple and best approach, using the actual earth distance of coordinates from the x and y direction. The algorithm works with any SRID, internally it works with WGS 1984(EPSG:4326), and the result transforms back to input 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, the geometry must be valid and polygon, multi-polygons or envelope type

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

Result enter image description here

2. Second function based on Nicklas Avén algorithm. I have implemented a way to handle any SRID.

   have upgraded the following changes in the algorithm.

   1. Separate variable for x and y direction for pixel size,
   2. New variable to calculate 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 enter image description here

3. Function-based on a 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

enter image description here

• Guven Degirmenci: Last function actually did the trick while the others fail or suffer from slow query. Thanks!

Answer 5

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!

• JaakL: Note that if you have large complex polygons areas (e.g. I have coastline buffer), then this approach is not quite optimal.
• julien: Then, what do you propose instead?

Answer 6

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 st_makegrid(the_geom,1000.0,100) as geom from the_polygon 
-- makegrid(the_geom,grid_size,alignement)

• Theo F: from an initial test, this method seems quite slow on more complex polygon geometries. Even when using a 1000m grid spacing...
• obchardon: @TheoF "quite slow" is also quite subjective. It should be compared with the others methods. On my mid range computer I can produced ~171000 points in less than 6 seconds.

Answer 7

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

• Michal Zimmermann: 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.

Answer 8

You can use ST_SquareGrid to create square polygons covering your layer, then calculate their centroids:

SELECT
 row_number() over() as id, 
 st_centroid(grid.geom) as geom
FROM 
 test.polygons poly
INNER JOIN
 ST_SquareGrid(150, poly.geom) as grid
ON 
 st_intersects(poly.geom, st_centroid(grid.geom))

enter image description here

Answer 9

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

• Nikolaus Krismer: Wouldn't transforming the geometry to a specific SRID (like EPSG:3857) be better than just multiplying by a scale factor?

Answer 10

Create the function first:

CREATE OR REPLACE FUNCTION st_polygrid(geometry, integer) RETURNS geometry AS
$$
SELECT
 ST_SetSRID(ST_Collect(ST_POINT(x,y)), ST_SRID(1ドル))
FROM
 generate_series(floor(ST_XMin(1ドル))::numeric, ceiling(ST_xmax(1ドル))::numeric, 2ドル) as x,
 generate_series(floor(ST_ymin(1ドル))::numeric, ceiling(ST_ymax(1ドル))::numeric,2ドル) as y 
WHERE
 ST_Intersects(1,ドルST_SetSRID(ST_POINT(x,y), ST_SRID(1ドル)))
$$
 LANGUAGE sql VOLATILE;

Then create the point grid using that function:

CREATE TABLE some.other_table AS
SELECT
 a.gid,
 st_polygrid(a.the_geom, 1000) AS the_geom
FROM
 some.table a

changing the '1000' to the desired point spacing.

Answer 11

Adding yet another answer with named parameters - for the given geometry, specify the delta between each point (same value for x & y). Uses nested query to separate point generation from filtering points to only those that are inside the given geometry.

CREATE
OR REPLACE FUNCTION makegrid(geom geometry, delta numeric)
 RETURNS geometry AS $$
SELECT ST_Collect(point) -- combine all matching points into a single multi-point geometry
FROM (
 -- generate a list of points within the given geometry's bounding box
 SELECT ST_SetSRID(ST_Point(x, y), ST_SRID(geom)) point
 FROM generate_series(
 floor(ST_XMIN(geom))::numeric,
 ceiling(ST_XMAX(geom))::numeric,
 delta) AS x,
 generate_series(
 floor(ST_YMIN(geom))::numeric,
 ceiling(ST_YMAX(geom))::numeric,
 delta) AS y) t
-- only keep the points that are actually inside the geometry
WHERE st_intersects(geom, point) $$ LANGUAGE sql;

• dr_jts: Minor point: you can now write ST_Point(x, y, ST_SRID(geom))

Answer 12

As commented before, the generate_series(min,interval,step) of Avén's solution is not general, fails for example for equal-area projections. Here a correction, with expected "from min to max", generate_series(min,max,step).

To be friendly with user that need "approximately n points per axis", instead "distance between the points in the grid": here also a new parameter suggestion.

To be useful: optional parameters to user-defined origin of the grid, and option to show all the rectangular grid.

For equal-area projections

Using a function name analog to ST_GeneratePoints.

The second parameter is the distance between points, or, if negative, the number of points per axis. Use for example -100 for 100 points in the minor BBBOX axis.

CREATE FUNCTION ST_GenerateGridPoints(
 g geometry, -- a geometry using SRID of equal-area projection
 p_dist integer, -- positive=distance between points;
 -- negative= number of points per axis (some excluded by intersection)
 p_x0 int default NULL, -- (optional) grid origin X
 p_y0 int default NULL, -- (optional) grid origin Y
 p_show_all boolean default false -- to show all rectangular grid
) RETURNS geometry AS $f$
 SELECT ST_Collect( ST_POINT(x, y, ST_SRID(1ドル)) )
 FROM ( --- Parameter calculation: ---
 SELECT CASE 
 WHEN p_dist>0 THEN p_dist 
 WHEN p_dist<0 THEN LEAST(
 ceiling(ST_XMAX(1ドル)-ST_XMIN(1ドル))::int / -p_dist,
 ceiling(ST_YMAX(1ドル)-ST_YMIN(1ドル))::int / -p_dist
 )
 ELSE NULL
 END,
 CASE WHEN p_x0 is null THEN floor(ST_XMIN(1ドル))::int ELSE p_x0 END,
 CASE WHEN p_y0 is null THEN floor(ST_YMIN(1ドル))::int ELSE p_y0 END
 ) t(selected_dist, x0, y0),
 generate_series( --- X axis scan: ---
 x0,
 ceiling(ST_XMAX(1ドル))::int,
 selected_dist
 ) AS x,
 generate_series( --- Y axis scan: ---
 y0,
 ceiling(ST_YMAX(1ドル))::int,
 selected_dist
 ) AS y 
 WHERE p_show_all OR st_intersects(1,ドル ST_POINT(x, y, ST_SRID(1ドル)))
$f$ language SQL;

The ST_POINT() will cast integers to floats, but "all integers" make sense, because (in general) the projected unit is 1 meter, enough precision.

---

Example using the geometry of Brazil (by conic Albers) and 30 points per minor BBOX axis, ST_GenerateGridPoints(geom,-30). Same as ST_GenerateGridPoints(geom,+86012).

enter image description here

Example using the geometry of Cameroon (by near-equal-area projection), with origin at (408600,164150) and distance of 150 km:

• Cameroon points: ST_GenerateGridPoints_list(geom,150000,408600,164150,false)
• All points: ST_GenerateGridPoints_list(geom,150000,408600,164150,true)

enter image description here

Answer 13

Here is my fancy version that lets you set padding to the edges, and centers the plots in the geometry (or geometry extent if its multi). I pass geojson, but it should be trivial to change that to a geometry.

CREATE OR REPLACE FUNCTION gridded_points_in_bounds(_geo_json jsonb, _m_spacing real, _m_buffer real)
 RETURNS table (
 lon float,
 lat float
 ) AS $$
 DECLARE
 _meters_boundary geometry;
 _buffered_extent geometry;
 _x_range float;
 _y_range float;
 _x_steps integer;
 _y_steps integer;
 _x_padding float;
 _y_padding float;
 BEGIN
 SELECT ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(_geo_json), 4326), 3857) INTO _meters_boundary;
 SELECT ST_Buffer(_meters_boundary, -1 * _m_buffer) INTO _buffered_extent;
 SELECT ST_XMax(_buffered_extent) - ST_XMin(_buffered_extent) INTO _x_range;
 SELECT ST_YMax(_buffered_extent) - ST_YMin(_buffered_extent) INTO _y_range;
 SELECT floor(_x_range / _m_spacing) INTO _x_steps;
 SELECT floor(_y_range / _m_spacing) INTO _y_steps;
 SELECT (_x_range - _x_steps * _m_spacing) / 2 INTO _x_padding;
 SELECT (_y_range - _y_steps * _m_spacing) / 2 INTO _y_padding;
 RETURN QUERY
 SELECT ST_X(ST_Centroid(geom)),
 ST_Y(ST_Centroid(geom))
 FROM (
 SELECT ST_Transform(
 ST_SetSRID(
 ST_POINT(x::float + _x_padding, y::float + _y_padding), ST_SRID(_buffered_extent)
 ),
 4326
 ) as geom
 FROM
 generate_series(floor(st_xmin(_buffered_extent))::int, ceiling(st_xmax(_buffered_extent))::int, _m_spacing::int) AS x,
 generate_series(floor(st_ymin(_buffered_extent))::int, ceiling(st_ymax(_buffered_extent))::int, _m_spacing::int) AS y
 WHERE ST_Intersects(
 _buffered_extent,
 ST_SetSRID(ST_POINT(x::float + _x_padding, y::float + _y_padding), ST_SRID(_buffered_extent))
 )
 ) a;
 END
$$ LANGUAGE PLPGSQL;

Answer 14

I didn't need a function, just a query I could run. Based on some of the answers here, I made the following:

WITH geom AS (
 SELECT st_setsrid('polygon-string'::geometry, 4326) AS geometry
),
points AS (
 SELECT generate_series(ST_XMIN(ST_Transform(geometry, 3857))::int, ST_XMAX(ST_Transform(geometry, 3857))::int, 1000) AS x,
 generate_series(ST_YMIN(ST_Transform(geometry, 3857))::int, ST_YMAX(ST_Transform(geometry, 3857))::int, 1000) AS y
 FROM geom
)
SELECT ST_Collect(ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326)) AS points
FROM points, geom
WHERE st_intersects(geometry, ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326));

1. The geom CTE casts the string to geometry and makes sure its SRID to 4326
2. The points CTE formats the geometry into 3857, and creates two lists of points 1000m apart
3. The final query transforms the points back into 4326, filters out points that aren't inside the polygon's shape, and collects them into a geometry collection

• dr_jts: I don't think this produces a grid, since the generate_series functions aren't evaluated independently.