1

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;
4
  • 1
    how are you measuring the reprojected fishnet? you would expect some distortion
    – Ian Turton
    Oct 18, 2022 at 16:02
  • Once its in 4326 I can use google earth, postgis, and qgis and they all show the same distortion. For my use case I cannot accept any distortion (or at least very very minimal), so I will need to find a different way to make the grid.
    – tbob
    Oct 18, 2022 at 16:08
  • can you post pictures of what you get and what you would like?
    – Ian Turton
    Oct 18, 2022 at 19:00
  • I actually got it to work, I just posted an answer. I'll do some more testing to make sure it works for the general case.
    – tbob
    Oct 18, 2022 at 19:27

2 Answers 2

0

Did you try the ST_SquareGrid function?

it's available since PostGIS 3.1

1
  • Thank you for the recommendation, I have not tried this, I'm looking into now, I actually do start with a rectangle geometry and typically I will pass into st_createfishnet the "bottom left" corner of the rectangle and then I will rotate the fishnet to line up with the rectangle. for squaregrid would I just pass in the size (100) and the recangle geometry and the resulting grid will already be rotated to align with the recatngel?
    – tbob
    Oct 18, 2022 at 17:26
0

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;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.