Shapely polygons crossing the antimeridian

Question

Sample data

Consider the following WKT Polygon, crossing the international dateline (antimeridian):

POLYGON((176 49,-65 49,-65 11,176 11,176 49))

And the following points:

POINT(-140 32) # Inside the polygon
POINT(0 32)    # Outside the polygon

The problem

Shapely considers this polygon to span on the other side of the planet - covering Asia and the Atlantic, rather than the US and the Pacific. Therefore, it fails to calculate its centroid and tell whether points are inside or outside it:

from shapely import wkt

polygon_wkt = 'POLYGON((176 49,-65 49,-65 11,176 11,176 49))'
point_in_polygon_wkt = 'POINT(-140 32)'
point_outside_polygon_wkt = 'POINT(0 32)'

polygon = wkt.loads(polygon_wkt)
point_in_polygon = wkt.loads(point_in_polygon_wkt)
point_outside_polygon = wkt.loads(point_outside_polygon_wkt)


print(polygon.centroid)                         # POINT (55.5 30) - Wrong!
print(polygon.contains(point_in_polygon))       # False - Wrong!
print(polygon.contains(point_outside_polygon))  # True - Wrong!

What have I tried

• Using PostGIS - I get the same erroneous results.
• Playing with Shapely arguments - couldn't manage to "wrap" the polygon to the other side of the planet.
• Reading The International Date Line wrap around. To be frank, there does not seem to be an answer there (Except for splitting the polygon).

My question

How can I calculate the centroid, bounding box, and inside/outside predicate for a WGS84 polygon optionally crossing the international dateline (longitude 180 / -180)?

Answer 1

You can create polygons with longitude that goes beyond -180°, in this case the western bound of your polygon will be at -184°.

from shapely import wkt

polygon_wkt = 'POLYGON((-184 49,-65 49,-65 11,-184 11,-184 49))'
point_in_polygon_wkt = 'POINT(-140 32)'
point_outside_polygon_wkt = 'POINT(0 32)'

polygon = wkt.loads(polygon_wkt)
point_in_polygon = wkt.loads(point_in_polygon_wkt)
point_outside_polygon = wkt.loads(point_outside_polygon_wkt)

print(polygon.centroid)          
POINT (-124.5 30)               
print(polygon.contains(point_in_polygon))
True
print(polygon.contains(point_outside_polygon))  
False

shapely does not know about the coordinate system, and it's x/y plane goes to infinity, so it naturally works like that.

But what about in a GIS? If I put the polygon inside geopandas is does not seem to mind.

import geopandas as gpd

polygon_geodf = gpd.GeoDataFrame({'geometry':gpd.GeoSeries([polygon])})
polygon_geodf = polygon_geodf.set_crs('EPSG:4326')
print(polygon_geodf.is_valid)
True

print(polygon_geodf.contains(point_in_polygon))
0    True
dtype: bool

print(polygon_geodf.contains(point_outside_polygon))
0    False
dtype: bool

Intuitively this does makes sense. From the prime meridian, where longitude=0, going 176° east or -184° west will put you in the same spot. Mathematically they are the same too, and will give the same sin and cosine values, to a high precision, after converting to radians.

from math import radians, sin, cos
sin(radians(176))
0.06975647374412552
sin(radians(-184))
0.06975647374412527

cos(radians(176))
-0.9975640502598242
cos(radians(-184))
-0.9975640502598242

But also geopandas does not actually do spherical geometry. So if there is a third point we know to be inside the polygon on a sphere, but it is outside the polygon on the planer coordinates, we are back to the original problem.

point_in_polygon2 = wkt.loads('POINT(178 32)')
polygon_geodf.contains(point_in_polygon2)
0    False
dtype: bool

The geopandas documentation does not say polygons crossing the dateline will "not work", but (during transformations at least) have "undesirable behavior".

Objects crossing the dateline (or other projection boundary) will have undesirable behavior.

In the end if you know your operations will be confined to a certain area, and that assuming a planer coordinate system will work, then you can extend the polygon boundaries past -180 (or 180 for that matter). But if your have points and polygons all over the world you must work with, then you probably need to implement fancier things like in the tutorial @pawarit28 linked to.

Answer 2

Centroid & Bounding Box:

• You'll first have to shift the entire geometry to avoid flipping signs along the antimeridian.
• You can next compute the centroid & bounding box in the new shifted space then finally add or subtract the offset to get the centroid & bounding box in the original space

Predicates

• Here, I think it's easiest to split the polygon along the antimeridian. Hopefully this tutorial + code snippets are helpful (uses Shapely)

Answer 3

For polygons that falls on the edge, vis around 180 degree longitude, I came with a solution that covers all those cases:

    from shapely.geometry import Polygon, MultiPolygon, LineString
    from shapely.ops import split
    def split_poly(ULLon, ULLat, LLLon, LLLat, LRLon, LRLat, URLon, URLat):
       ULLon = 360 + ULLon if ULLon < 0 else ULLon
       URLon = 360 + URLon if URLon < 0 else URLon
       LRLon = 360 + LRLon if LRLon < 0 else LRLon
       LLLon = 360 + LLLon if LLLon < 0 else LLLon
       pol = Polygon([[ULLon, ULLat], [LLLon, LLLat], [LRLon, LRLat], [URLon, URLat], [ULLon, ULLat]])
       line = LineString([(180, 90), (180, -90)])
       spl_pol = split(pol, line)
       for pol in spl_pol:
           if True in [item>180 for item in pol.exterior.coords.xy[0]]:
                new_lon = [item-360 if item >= 180 else item for item in pol.exterior.coords.xy[0]]
                lat = pol.exterior.coords.xy[1]
                poly1 = Polygon(zip(*(new_lon, lat)))
           else:
               poly2 = pol
       return poly1, poly2