4

Sample data

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

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

Polygon crossing the dateline

And the following points:

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

Points inside and 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)?

  • 2
    Shapely uses a cartesian plane system for computing geometries (distance = euclidean distance). That means that if you work with a crs.unit = degree (WGS84 for example) all calculations are wrong.You must first reproject you layer (many examples in GIS SE) – gene Feb 1 '17 at 9:38
  • Of course, but I think that projection is not the problem - the center point is wrong in PostGIS too. It boils down to choosing the right polygon direction - whether it spans on the western or eastern hemisphere. – Adam Matan Feb 1 '17 at 19:43
  • The principal problem is that the units are angular (degrees) and not cartesian – gene Feb 1 '17 at 20:48
  • I think you just have to split it into two. @gene I don't see how the units make a bit of difference; he's not computing geometry attributes. – Jon Apr 3 '18 at 14:46
  • Using PostGIS and QGIS I do not observe your problem, the centroid of the polygon using QGIS is defined correctly, with respect – Cyril May 25 '18 at 6:38

Your Answer

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

Browse other questions tagged or ask your own question.