# PostGIS - Polygon crossing latitude 0 and longitude 180

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: 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: 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: Is there a way to get teh result of the first picture?

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): 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: 