Converts a GeoJSON polygon to its antimeridian-compatible constituent polygon(s)

example_input.py

geojson = {
  "type": "Polygon",
  "coordinates": [
    [
      [179.0, 0.0], 
      [-179.0, 0.0], 
      [-179.0, 1.0],
      [179.0, 1.0], 
      [179.0, 0.0]
    ]
  ]
}

geopolygon_utils.py

from typing import Union, List, Generator
from shapely.geometry import mapping, Polygon, GeometryCollection
from shapely import affinity

def check_crossing(lon1: float, lon2: float, validate: bool = False, dlon_threshold: float = 180.0):
    """
    Assuming a minimum travel distance between two provided longitude coordinates,
    checks if the 180th meridian (antimeridian) is crossed.
    """

    if validate and any([abs(x) > 180.0 for x in [lon1, lon2]]):
        raise ValueError("longitudes must be in degrees [-180.0, 180.0]")   
    return abs(lon2 - lon1) > dlon_threshold

def translate_polygons(geometry_collection: GeometryCollection, 
                       output_format: str = "geojson") -> Generator[
                          Union[List[dict], List[Polygon]], None, None
                       ]:
    
  for polygon in geometry_collection:
      (minx, _, maxx, _) = polygon.bounds
      if minx < -180: geo_polygon = affinity.translate(polygon, xoff = 360)
      elif maxx > 180: geo_polygon = affinity.translate(polygon, xoff = -360)
      else: geo_polygon = polygon

      yield_geojson = output_format == "geojson"
      yield json.dumps(mapping(geo_polygon)) if (yield_geojson) else geo_polygon

polygon_splitter.py

import math
import copy
import json
from typing import Union, List
from shapely.geometry import Polygon, LineString, GeometryCollection
from shapely.ops import split

# https://gist.github.com/PawaritL/ec7136c0b718ca65db6df1c33fd1bb11
from geopolygon_utils import check_crossing, translate_polygons

def split_polygon(geojson: dict, output_format: str = "geojson", validate: bool = False) -> Union[
    List[dict], List[Polygon], GeometryCollection
    ]:
    """
    Given a GeoJSON representation of a Polygon, returns a collection of
    'antimeridian-safe' constituent polygons split at the 180th meridian, 
    ensuring compliance with GeoJSON standards (https://tools.ietf.org/html/rfc7946#section-3.1.9)
    Assumptions:
      - Any two consecutive points with over 180 degrees difference in
        longitude are assumed to cross the antimeridian
      - The polygon spans less than 360 degrees in longitude (i.e. does not wrap around the globe)
      - However, the polygon may cross the antimeridian on multiple occasions
    Parameters:
        geojson (dict): GeoJSON of input polygon to be split. For example:
                        {
                        "type": "Polygon",
                        "coordinates": [
                          [
                            [179.0, 0.0], [-179.0, 0.0], [-179.0, 1.0],
                            [179.0, 1.0], [179.0, 0.0]
                          ]
                        ]
                        }
        output_format (str): Available options: "geojson", "polygons", "geometrycollection"
                             If "geometrycollection" returns a Shapely GeometryCollection.
                             Otherwise, returns a list of either GeoJSONs or Shapely Polygons
        validate (bool): Checks if all longitudes are within [-180.0, 180.0]
      
    Returns:
        List[dict]/List[Polygon]/GeometryCollection: antimeridian-safe polygon(s)
    """
    output_format = output_format.replace("-", "").strip().lower()
    coords_shift = copy.deepcopy(geojson["coordinates"])
    shell_minx = shell_maxx = None
    split_meridians = set()
    splitter = None

    for ring_index, ring in enumerate(coords_shift):
        if len(ring) < 1: 
            continue
        else:
            ring_minx = ring_maxx = ring[0][0]
            crossings = 0

        for coord_index, (lon, _) in enumerate(ring[1:], start=1):
            lon_prev = ring[coord_index - 1][0] # [0] corresponds to longitude coordinate
            if check_crossing(lon, lon_prev, validate=validate):
                direction = math.copysign(1, lon - lon_prev)
                coords_shift[ring_index][coord_index][0] = lon - (direction * 360.0)
                crossings += 1

            x_shift = coords_shift[ring_index][coord_index][0]
            if x_shift < ring_minx: ring_minx = x_shift
            if x_shift > ring_maxx: ring_maxx = x_shift

        # Ensure that any holes remain contained within the (translated) outer shell
        if (ring_index == 0): # by GeoJSON definition, first ring is the outer shell
            shell_minx, shell_maxx = (ring_minx, ring_maxx)
        elif (ring_minx < shell_minx):
            ring_shift = [[x + 360, y] for (x, y) in coords_shift[ring_index]]
            coords_shift[ring_index] = ring_shift
            ring_minx, ring_maxx = (x + 360 for x in (ring_minx, ring_maxx))
        elif (ring_maxx > shell_maxx):
            ring_shift = [[x - 360, y] for (x, y) in coords_shift[ring_index]]
            coords_shift[ring_index] = ring_shift
            ring_minx, ring_maxx = (x - 360 for x in (ring_minx, ring_maxx))

        if crossings: # keep track of meridians to split on
            if ring_minx < -180: split_meridians.add(-180)
            if ring_maxx > 180: split_meridians.add(180)

    n_splits = len(split_meridians)
    if n_splits > 1:
        raise NotImplementedError(
            """Splitting a Polygon by multiple meridians (MultiLineString) 
               not supported by Shapely"""
        )
    elif n_splits == 1:
        split_lon = next(iter(split_meridians))
        meridian = [[split_lon, -90.0], [split_lon, 90.0]]
        splitter = LineString(meridian)

    shell, *holes = coords_shift if splitter else geojson["coordinates"]
    if splitter: split_polygons = split(Polygon(shell, holes), splitter)
    else: split_polygons = GeometryCollection([Polygon(shell, holes)])
        
    geo_polygons = list(translate_polygons(split_polygons, output_format))  
    if output_format == "geometrycollection": return GeometryCollection(geo_polygons)
    else: return geo_polygons

AttributeError: 'str' object has no attribute 'bounds' is an error I got while using split_polygons().

I really need this code to work for me, I have a GeoDataFrame and can easily obtain a GeoJSON of it. I took the GeoJSON of a single polygon and passed it in, that is when I got the error above.

With a GeoDataFrame you need to use gdf.geo_interface['features'][index]['geometry'] . The GeoJSON like object that GeoPandas returns by default is a feature collection which is not what this function wants. You need to pass the geometry associated with a single polygon at a time. So you'll have to loop over all the features in the dictionary, grab their geometries independently, and then pass them to this method. By the way, you might also need to cast the coordinates in the returned geometry dictionary to a list instead of a tuple.

ah sorry, i completely missed @tony-zeidan comments! (haven't checked my github in a while)

+1 to @jgodwin's answer. the functions i wrote were intended for the GeoJSON standard (using Python dictionaries and lists), which might be different to the Python objects that GeoPandas produces. however, conversions should be quite doable. feel free to share an example of the problem you're facing on Google Colab and i'll be happy to take a look as well!

Thanks @PawaritL this is very useful code. I am just learning about antimeridian crossings myself. It's impressive how you dealt with holes etc.

I captured the above gists and started writing some python packaging, and some pytests. Also I created 1 other issue I noticed.
https://github.com/guidorice/antimeridian_splitter

Pawarit, if you want to take ownership of the repository, or add a License file, please let me know, glad to hand it to you.

Hi @guidorice, thanks so much for catching that!
It was a typo when I tried to make a quick edit directly on the gist. sorry about that
I've fixed it now and would be happy to discuss the issue you saw on geojson.io