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
)?
Best Answer
You can create polygons with longitude that goes beyond -180°, in this case the western bound of your polygon will be at -184°.
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.
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.
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.
The geopandas documentation does not say polygons crossing the dateline will "not work", but (during transformations at least) 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.