Python OGR – Handling Polygons Across the International Dateline with Python and OGR

antimeridianogrpythonshapelywgs84

I'm trying to generate polygons for satellite orbital swaths.
So far I have a method to generate two lines which represent the edge of each swath in [lat,long]. Some of the swath's cross the international dateline and so wrap round:

I was able to solve this with ogr2ogr -wrapdateline:

ogr2ogr -wrapdateline  -f "ESRI Shapefile" test.shp orbits.shp

Which does split the lines probably

I now want to be able to generate polygons on the interior of both lines. So for example in the case where one edge of the swath crosses the dateline a polygon fills in when it emerges on the other side, like:

I need a method that is automated as I need to repeat the task a lot. Preferably in python as that's how i have generated the lines. Here are the two shapefiles containing the lines: wraparound ; datelinefixed

Best Answer

You can build a custom mercator projection centered approximately on the center of the swath. For example, use for swath 25:

+proj=merc +lon_0=-140 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

In this projection, the swath is not broken by the dateline. You can create the polygon from the line.

Then create a cut polygon between -179.95°E and 179.95°E in EPSG:4326:

Nr;WKT
1;POLYGON ((-179.95 89, 179.95 89, 179.95 -89, -179.95 -89, -179.95 89))

Reproject it to your custom CRS too, and subtract it from the swath polygon.

After reprojecting back to EPSG:4326, the swath is correctly divided by the dateline:

Continue with all swaths that cross the dateline.

Related Solutions

Shapefile – Changing Longitude Range from -180-180 to 0-360

Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite:

ogr2ogr world_shifted.shp world.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,180,0) FROM world"

or:

ogr2ogr -s_srs EPSG:4326 -t_srs "+proj=longlat +ellps=WGS84 +pm=-180 +datum=WGS84 +no_defs" world_shifted.shp world.shp

Both commands produce a longitude offset of 180°, i.e. a prime meridian of -180° is considered. In fact:

>ogrinfo world_shifted.shp world_shifted | grep Extent
Extent: (0.000000, -90.000000) - (360.000000, 83.623596)

The difference between the two commands is that with a longitude offset (2nd try) data are simply reprojected using -180° as prime meridian, while shifting the coordinates geometries (1st try) are altered, even if the result is apparently the same.

EDIT
If there are parts in 0-180 that should not move, it's possible to adapt this working solution: https://gis.stackexchange.com/a/73164/22405

Clip the two parts:

ogr2ogr world_part1.shp world.shp -clipsrc -180 -90 0 90
ogr2ogr world_part2.shp world.shp -clipsrc 0 -90 180 90

Shift only the first part:

ogr2ogr world_part1_shifted.shp world_part1.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,360,0), CNTRY_NAME FROM world_part1"

Then, merge the second part and the first shifted:

ogr2ogr world_0_360_raw.shp world_part2.shp
ogr2ogr -update -append world_0_360_raw.shp world_part1_shifted.shp -nln world_0_360_raw

Finally, dissolve countries boundaries of world_0_360_raw.shp obtaining world_0_360.shp by country names. For instance:

ogr2ogr world_0_360.shp world_0_360_raw.shp -dialect sqlite -sql "SELECT ST_Union(Geometry), CNTRY_NAME FROM world_0_360_raw GROUP BY CNTRY_NAME"

[GIS] Shapely polygons crossing the antimeridian

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.

Best Answer

Related Solutions

Shapefile – Changing Longitude Range from -180-180 to 0-360

[GIS] Shapely polygons crossing the antimeridian

Related Question