0

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.

New contributor
alxgrh is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
  • What if you use some other format than GeoJSON as input? Non-splited polygons around the dataline are not strictly invalid according to tools.ietf.org/html/rfc7946#page-10 but it is recommended that they should not exist Any geometry that crosses the antimeridian SHOULD be represented by cutting it in two such that neither part's representation crosses the antimeridian.. – user30184 Oct 15 at 21:33
  • Input GeoJSON is preliminary loaded into osgeo.gdal.Dataset, which is unified format. I need to split polygons around antimeridian for displaying. That's why I'm trying to use the VectorTranlate. – alxgrh Oct 17 at 8:56
0

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
| improve this answer | |
New contributor
alxgrh is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
0

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.

| improve this answer | |

Your Answer

alxgrh is a new contributor. Be nice, and check out our Code of Conduct.

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.