3

I need to find all the couples long/lat (listed in a table) that belong to the area defined by a polygon in PostGIS (version 3.0). My area is defined by the following polygon, which crosses the meridian with longitude 180 and the parallel with latitude 0:

POLYGON((1 1, 70 70, 100 -80, -130 -70, 130 40, -100 70, 1 1 ))

The surface is plotted in this image:

enter image description here

The reference table public.tiles has one row for each disjoint square like POLYGON((-180 -90,-180 -87.5,-177.5 -87.5,-177.5 -90,-180 -90)), POLYGON((-180 -87.5,-180 -85,-177.5 -85,-177.5 -87.5,-180 -87.5)), etc. (written in the column tile), in order to cover the entire map.

To get the coordinates that belong to the area, I tried to use the function ST_Intersects, as in the following query:

SELECT
    St_x(St_centroid(t.tile)) AS st_x,
    St_y(St_centroid(t.tile)) AS st_y
FROM 
    (SELECT ST_GeomFromText('POLYGON((1 1,  70 70,  100 -80,  -130 -70,  130 40,  -100 70,  1 1  ))',4326) AS footprint) AS p
    JOIN  public.tiles  AS y
        ON ST_Intersects(ST_Transform(ST_SetSRID(p.footprint, 4326), 4326),t.tile)
;

When I plot the result, it is really unexpected: enter image description here

Guessing that the problem could be about the crossing of the "borders", I added the function ST_ShiftLongitude to the query:

SELECT
    St_x(St_centroid(t.tile)) AS st_x,
    St_y(St_centroid(t.tile)) AS st_y
FROM 
    (SELECT ST_GeomFromText('POLYGON((1 1,  70 70,  100 -80,  -130 -70,  130 40,  -100 70,  1 1  ))',4326) AS footprint) AS p
    JOIN  public.tiles  AS y
        ON ST_Intersects(ST_ShiftLongitude(ST_Transform(ST_SetSRID(p.footprint, 4326), 4326)),ST_ShiftLongitude(t.tile))
;

Eventually, the result is closer to what I am looking for: enter image description here

Is there a way to get teh result of the first picture?

1 Answer 1

3

You're right with the fact that crossing the borders is the problem. And you're crossing 2 borders: the 180° and the north pole.

The thing is, here you work in 4326. When you build the polygon from the points, it will take the shortest path between points, but not by crossing borders. That's why you end up with the "wrong" figure.

Actually, if you think about it, there needs to be a rule like that, otherwise there is no way to distinguish between the interior and exterior: in your exemple, for what reason the interior would be the above of the line and not the below ?

So if you need to work with polygons all around the world and the north pole, I think the most simple solution is to work in a referential system based around the north pole, for exemple 102016 (depending of what kind of precision you want, you can use others).

You will have with you exemple something like that:

select 
st_astext(
    ST_MakePolygon(ST_MakeLine(
        ARRAY[
            ST_Transform(ST_GeomFromText('POINT(1 1)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(-100 70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(130 40)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(-130 -70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(100 -80)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(70 70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(1 1)', 4326), 102016)
        ]
    ))
)

And the result (display in the correct projection of course, otherwise the transformation will cut your poly again):

enter image description here

If you really want to go back to 4326, maybe you can try cutting the polygon by the rectangle -180 180 -90 90 and reproject the result multipolygon.

EDIT:

I made an example with the point generation and selection:

with poly as (
select 
    ST_MakePolygon(ST_MakeLine(
        ARRAY[
            ST_Transform(ST_GeomFromText('POINT(1 1)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(-100 70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(130 40)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(-130 -70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(100 -80)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(70 70)', 4326), 102016),
            ST_Transform(ST_GeomFromText('POINT(1 1)', 4326), 102016)
        ]
    ))
    AS geom
),
t1 as (select generate_series(-179, 180) as lon),
t2 as (select generate_series(-90, 90) as lat),
points as (select * from t1 cross join t2),
points_geom as (select st_setsrid(st_makepoint(lon, lat), 4326) as geom from points)
select 
    st_astext(st_collect(pg.geom))
    from points_geom as pg, poly
    where st_intersects(st_transform(pg.geom, 102016), poly.geom)

Here is the result:

Points selected

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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