How to create grid in 4326 just like ST_CreateFishnet PostGIS

Question

So I've been projecting my geometries to 102008, then using this function:

            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;

and then making a grid with grid cell size 100m x 100m to overlay my geometries.

The problem is when I project back to 4326 the size is distorted to 95m x 105m.(some examples are distorted more, some less)

I think I should avoid projecting in the first place to avoid this. How can I modify the st_createfishnet to make a 4326 grid, and still be able to pass in the width and height in meters?

So far my attempt has been to pass in the point, and then buffer it to create the cell, I'll always be using squares for cells. but this attempt so far is not working:

        CREATE OR REPLACE FUNCTION ST_CreateFishnet(
            nrow integer, ncol integer,
            xsize float8, ysize float8,
            x0 geometry,
            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 * (st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 1),
                        ST_PointN(ST_ExteriorRing(cell), 2)))/2), i * (st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 1),
                        ST_PointN(ST_ExteriorRing(cell), 2)))/2)) 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
    SELECT st_envelope(st_buffer($5, $4))::geometry AS cell
    ) AS foo
    $$ LANGUAGE sql IMMUTABLE STRICT;
    
    select (ST_CreateFishnet(2, 3, 100, 100, st_geometryfromtext('POINT(-117.20051364779485 53.33216281083014)'))).geom;

Answer 1

Did you try the ST_SquareGrid function?

it's available since PostGIS 3.1

Answer 2

I was able to modify the function like this and it worked:

        CREATE OR REPLACE FUNCTION ST_CreateFishnet(
            nrow integer, ncol integer,
            xsize float8, ysize float8,
            x0 geometry,
            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 * (st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 2),
                        ST_PointN(ST_ExteriorRing(cell), 3))))+(st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 2),
                        ST_PointN(ST_ExteriorRing(cell), 3)))/2), i * (st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 1),
                        ST_PointN(ST_ExteriorRing(cell), 2))))+(st_length(ST_MakeLine(ST_PointN(ST_ExteriorRing(cell), 1),
                        ST_PointN(ST_ExteriorRing(cell), 2)))/2)) 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
    SELECT st_envelope((st_buffer($5::geography, $4/2))::geometry)::geometry AS cell
    ) AS foo
    $$ LANGUAGE sql IMMUTABLE STRICT;