gdal.VectorTranslate fails to wrap the antimeridian on specific input

Question

I'm trying to split a GeoJSON polygon which crosses the antimeridian. I use this solution as base gdal.VectorTranslate returns an empty file

My code is:

geojson_str = <some geojson string>
geojson = json.loads(gjstr)
print('Original GeoJSON:')
print(json.dumps(geojson, indent=2))
with open('in.geojson','w') as f:
    json.dump(geojson, f)
gdal.UseExceptions()
g = gdal.OpenEx('in.geojson')
out = gdal.VectorTranslate('./out.geojson', g, format = 'GeoJSON', layerCreationOptions = ["-wrapdateline", '-lco','RFC7946=YES'])
del out
geojson_transformed = json.load(open('./out.geojson'))
print('Transformed GeoJSON:')
print(json.dumps(geojson_transformed, indent=2))

The issue is that gdal.VectorTranslate() has indeterministic behaviuor on different input.

If I use GeoJSON from example, it splits polygon into multipolygon correctly:

Original GeoJSON:
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              172.56933593749997,
              61.58549218152362
            ],
            [
              186.0205078125,
              61.58549218152362
            ],
            [
              186.0205078125,
              66.00015035652659
            ],
            [
              172.56933593749997,
              66.00015035652659
            ],
            [
              172.56933593749997,
              61.58549218152362
            ]
          ]
        ]
      }
    }
  ]
}
Transformed GeoJSON:
{
  "type": "FeatureCollection",
  "name": "in",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                180.0,
                61.5854922
              ],
              [
                180.0,
                66.0001504
              ],
              [
                172.5693359,
                66.0001504
              ],
              [
                172.5693359,
                61.5854922
              ],
              [
                180.0,
                61.5854922
              ]
            ]
          ],
          [
            [
              [
                -180.0,
                66.0001504
              ],
              [
                -180.0,
                61.5854922
              ],
              [
                -173.9794922,
                61.5854922
              ],
              [
                -173.9794922,
                66.0001504
              ],
              [
                -180.0,
                66.0001504
              ]
            ]
          ]
        ]
      }
    }
  ]
}

The example input contains shifted coordinates, where all negative longitudes incresed by 360. The output contains non-shifted coordinates though.

If I try to convert the example coordinates to be non-shifted, gdal.VectorTranslate() doesn't split the polygon:

Original GeoJSON:
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              172.56933593749997,
              61.58549218152362
            ],
            [
              -173.9794921875,
              61.58549218152362
            ],
            [
              -173.9794921875,
              66.00015035652659
            ],
            [
              172.56933593749997,
              66.00015035652659
            ],
            [
              172.56933593749997,
              61.58549218152362
            ]
          ]
        ]
      }
    }
  ]
}
Transformed GeoJSON:
{
  "type": "FeatureCollection",
  "name": "in",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              172.5693359,
              61.5854922
            ],
            [
              172.5693359,
              66.0001504
            ],
            [
              -173.9794922,
              66.0001504
            ],
            [
              -173.9794922,
              61.5854922
            ],
            [
              172.5693359,
              61.5854922
            ]
          ]
        ]
      }
    }
  ]
}

So I can assume that dateline splitting is being performed only if coordinates are centered at longitude 180 (or -180).

Nevertheless when I try to use it on my custom Geojson with shifted coords, then gdal.VectorTranslate() fails to complete:

Original GeoJSON:
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              182.24833333333333,
              61.336666666666666
            ],
            [
              179,
              60.0
            ],
            [
              169.0,
              54.0
            ],
            [
              158.9961111111111,
              50.10194444444445
            ],
            [
              157.8,
              52.25
            ],
            [
              157.5361111111111,
              53.325833333333335
            ],
            [
              158.18333333333334,
              53.79666666666667
            ],
            [
              156.75,
              55.45
            ],
            [
              153.9897222222222,
              58.01416666666667
            ],
            [
              155.9,
              58.86666666666667
            ],
            [
              160.16666666666666,
              60.583333333333336
            ],
            [
              162.21666666666667,
              61.43333333333333
            ],
            [
              164.25,
              62.083333333333336
            ],
            [
              174.75,
              61.833333333333336
            ],
            [
              182.24833333333333,
              61.336666666666666
            ]
          ]
        ]
      },
      "properties": {
        "designator": "UHPP",
        "color": "gray"
      }
    }
  ]
}

RuntimeError                              Traceback (most recent call last)
<ipython-input-218-e81a49bf4254> in <module>()
     33 gdal.UseExceptions()
     34 g = gdal.OpenEx('in.geojson')
---> 35 out = gdal.VectorTranslate('./out.geojson', g, format = 'GeoJSON', layerCreationOptions = ["-wrapdateline", '-lco','RFC7946=YES'])
     36 del out
     37 geojson_transformed = json.load(open('./out.geojson'))
1 frames
/usr/lib/python3/dist-packages/osgeo/gdal.py in wrapper_GDALVectorTranslateDestName(args)
   3139 def wrapper_GDALVectorTranslateDestName(args):
   3140     """wrapper_GDALVectorTranslateDestName(char const * dest, Dataset srcDS, GDALVectorTranslateOptions options, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 3141     return _gdal.wrapper_GDALVectorTranslateDestName(*args)
   3142 class GDALDEMProcessingOptions(_object):
   3143     """Proxy of C++ GDALDEMProcessingOptions class."""
RuntimeError: Terminating translation prematurely after failed
translation of layer in (use -skipfailures to skip errors)

If I use non-shifted coordintates in GeoJSON input, gdal.VectorTranslate() doesn't fail, but produce untransformed result:

Original GeoJSON:
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -177.75166666666667,
              61.336666666666666
            ],
            [
              179,
              60.0
            ],
            [
              169.0,
              54.0
            ],
            [
              158.9961111111111,
              50.10194444444445
            ],
            [
              157.8,
              52.25
            ],
            [
              157.5361111111111,
              53.325833333333335
            ],
            [
              158.18333333333334,
              53.79666666666667
            ],
            [
              156.75,
              55.45
            ],
            [
              153.9897222222222,
              58.01416666666667
            ],
            [
              155.9,
              58.86666666666667
            ],
            [
              160.16666666666666,
              60.583333333333336
            ],
            [
              162.21666666666667,
              61.43333333333333
            ],
            [
              164.25,
              62.083333333333336
            ],
            [
              174.75,
              61.833333333333336
            ],
            [
              -177.75166666666667,
              61.336666666666666
            ]
          ]
        ]
      },
      "properties": {
        "designator": "UHPP",
        "color": "gray"
      }
    }
  ]
}
Transformed GeoJSON:
{
  "type": "FeatureCollection",
  "name": "in",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "designator": "UHPP",
        "color": "gray"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -177.7516667,
              61.3366667
            ],
            [
              174.75,
              61.8333333
            ],
            [
              164.25,
              62.0833333
            ],
            [
              162.2166667,
              61.4333333
            ],
            [
              160.1666667,
              60.5833333
            ],
            [
              155.9,
              58.8666667
            ],
            [
              153.9897222,
              58.0141667
            ],
            [
              156.75,
              55.45
            ],
            [
              158.1833333,
              53.7966667
            ],
            [
              157.5361111,
              53.3258333
            ],
            [
              157.8,
              52.25
            ],
            [
              158.9961111,
              50.1019444
            ],
            [
              169.0,
              54.0
            ],
            [
              179.0,
              60.0
            ],
            [
              -177.7516667,
              61.3366667
            ]
          ]
        ]
      }
    }
  ]
}

Is it possible to determine why gdal.VectorTranslate() fails on the third input?

Are there other reliable and precise solutions to split polygons around antimeridian? I tried geopandas-based solution from here GDAL VectorTranslate creates shards/fragments along anti-meridian, but it doesn't seem to work.

Answer 1

Okay, I haven't found any suitable solution for splitting custom polygon by the antimeridian. I implemented it by myself using Pyproj and Shapely. Maybe someone will find it useful.

from pyproj import CRS, Transformer
from shapely.geometry import Point, LineString
def intersects_antimeridian_precise(start_point, end_point):
    '''
    Detects if line intersects the Antimeridian. Use custom Mercator projection for convenient result.
    '''
    # Initialize custom Mercator projection
    mercator_am_proj = CRS.from_proj4("+proj=merc lon_0={}".format(start_point.x))
    mercator_am_tf = Transformer.from_crs(mercator_am_proj.geodetic_crs, mercator_am_proj)
    mercator_am_tf_inv = Transformer.from_crs( mercator_am_proj,mercator_am_proj.geodetic_crs)
# Project Antimeridian to the Mercator plane
antimeridian= [(180,89.99),(180,-89.99)]
am_sp_m = Point(mercator_am_tf.transform(antimeridian[0][0],antimeridian[0][1]))
am_ep_m = Point(mercator_am_tf.transform(antimeridian[1][0],antimeridian[1][1]))
antimeridian_merc =  LineString([am_sp_m,am_ep_m])

# Project tested line to the Mercator plane
sp_m = Point(mercator_am_tf.transform(start_point.x, start_point.y))
ep_m = Point(mercator_am_tf.transform(end_point.x, end_point.y))
line_m = LineString([sp_m,ep_m])

return line_m.intersects(antimeridian_merc)


def wrap_antimeridian(polygon):
    ''' 
    Splits polygon which intersects the Antimeridian into two new polygons stitched to Antimeridian
    Params:
    polygon - list of (x,y) tuples
Returns:

List of two lists of (x,y) tuples or original polygon if it doesn't intersect the Antimeridian
'''

# Initialize Mercator transformer
mercator_am_proj = CRS.from_proj4("+proj=merc lon_0=179")
mercator_am_tf = Transformer.from_crs(mercator_am_proj.geodetic_crs, mercator_am_proj)
mercator_am_tf_inv = Transformer.from_crs( mercator_am_proj,mercator_am_proj.geodetic_crs)

antimeridian= [(180,89.99),(180,-89.99)]
am_sp_m = Point(mercator_am_tf.transform(antimeridian[0][0],antimeridian[0][1]))
am_ep_m = Point(mercator_am_tf.transform(antimeridian[1][0],antimeridian[1][1]))
antimeridian_merc =  LineString([am_sp_m,am_ep_m])

forw_intersect_idx = None
forw_intersect_line = None
back_intersect_idx = None
back_intersect_line = None

# Direction in which forward line intersects the Antimeridian. Use 1 for West-East and -1 for East-West
forw_intersect_direction = None

for i,p in enumerate(polygon):
    start_point = Point(p[0],p[1])
    if i == len(polygon)-1:
        end_point = Point(*polygon[0])
    else:
        end_point = Point(*polygon[i+1])

    # Find intersection points
    if intersects_antimeridian_precise(start_point, end_point):
        if forw_intersect_idx is None:
            forw_intersect_idx = i
            forw_intersect_line = [start_point, end_point]
        else:
            back_intersect_idx = i        
            back_intersect_line = [start_point, end_point]

if forw_intersect_idx is not None and back_intersect_idx is not None:
    # Convert lat,lon coordinates to Mercator cartesian plane and find forward and backward intersection points
    sp_f_merc = Point(mercator_am_tf.transform(*forw_intersect_line[0].coords[0]))
    ep_f_merc = Point(mercator_am_tf.transform(*forw_intersect_line[1].coords[0]))
    forw_intersect_line_merc = LineString([sp_f_merc,ep_f_merc])
    forw_intersect_point_merc = forw_intersect_line_merc.intersection(antimeridian_merc)
    forw_intersect_point = mercator_am_tf_inv.transform(forw_intersect_point_merc.x,forw_intersect_point_merc.y)

    sp_b_merc = Point(mercator_am_tf.transform(*back_intersect_line[0].coords[0]))
    ep_b_merc = Point(mercator_am_tf.transform(*back_intersect_line[1].coords[0]))
    back_intersect_line_merc = LineString([sp_b_merc,ep_b_merc])
    back_intersect_point_merc = back_intersect_line_merc.intersection(antimeridian_merc)
    back_intersect_point = mercator_am_tf_inv.transform(back_intersect_point_merc.x,back_intersect_point_merc.y)

    # Determine intersection direction
    if sp_f_merc.x < forw_intersect_point_merc.x:
        forw_intersect_direction = 1
    else:
        forw_intersect_direction = -1

    # Cut polygon before first intersection line and append first intersection point. Set longitude sign according to direction
    poly_1 = polygon[:forw_intersect_idx+1]
    poly_1.append((forw_intersect_point[0]*forw_intersect_direction, forw_intersect_point[1] ))

    # Cut polygon after second intersection and append second intersection point before
    poly_1.append((back_intersect_point[0]*forw_intersect_direction, back_intersect_point[1] ))
    poly_1 += polygon[back_intersect_idx+1:]


    # Cut polygon between intersection points, insert intersection points on its edges
    poly_2 = [(forw_intersect_point[0]*forw_intersect_direction*-1, forw_intersect_point[1])] + polygon[forw_intersect_idx+1:back_intersect_idx+1] + [(back_intersect_point[0]*forw_intersect_direction*-1, back_intersect_point[1])]

    return [poly_1, poly_2]
return polygon

Answer 2

I made a test with the json that fails for you. I mean this one

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              182.24833333333333,
              61.336666666666666
            ],
            [
              179,
              60.0
            ],
            [
              169.0,
              54.0
            ],
            [
              158.9961111111111,
              50.10194444444445
            ],
            [
              157.8,
              52.25
            ],
            [
              157.5361111111111,
              53.325833333333335
            ],
            [
              158.18333333333334,
              53.79666666666667
            ],
            [
              156.75,
              55.45
            ],
            [
              153.9897222222222,
              58.01416666666667
            ],
            [
              155.9,
              58.86666666666667
            ],
            [
              160.16666666666666,
              60.583333333333336
            ],
            [
              162.21666666666667,
              61.43333333333333
            ],
            [
              164.25,
              62.083333333333336
            ],
            [
              174.75,
              61.833333333333336
            ],
            [
              182.24833333333333,
              61.336666666666666
            ]
          ]
        ]
      },
      "properties": {
        "designator": "UHPP",
        "color": "gray"
      }
    }
  ]
}

Test with ogr2ogr into GeoJSON and OpenJUMP JML formats:

ogr2ogr -f jml split.jml split.json -wrapdateline
ogr2ogr -f geojson split2.json split.json -wrapdateline -lco RFC7946=YES

In both cases the result is a multipolygon that is split by the dateline.

MULTIPOLYGON (((180.0 60.411493073371,179 60,169 54,158.996111111111 50.1019444444444,157.8 52.25,157.536111111111 53.3258333333333,158.183333333333 53.7966666666667,156.75 55.45,153.989722222222 58.0141666666667,155.9 58.8666666666667,160.166666666667 60.5833333333333,162.216666666667 61.4333333333333,164.25 62.0833333333333,174.75 61.8333333333333,180.0 61.4855893902349,180.0 60.411493073371)),((-177.751666666667 61.3366666666667,-180 60.411493073371,-180 61.4855893902349,-177.751666666667 61.3366666666667)))

I do not know what is the reason for the failure with your code but GDAL seems to be able to handle your input with ogr2ogr.