GDAL VectorTranslate creates shards/fragments along anti-meridian

Question

I have challenge to take a GeoJSON file and convert the CRS to WGS84/EPSG:4326 and also split the polygons along the anti-meridian. These are both in preparation for presenting the data via mapbox.com.

I came across the VectorTranslate method of gdal which seems to miraculously achieve both of those requirements in 3 lines of code, e.g.

gdal.UseExceptions()
srcDS = gdal.OpenEx(datadir + 'input.geojson')
ds = gdal.VectorTranslate(datadir + 'output.geojson', srcDS=srcDS, format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

My problem is the results are not quite perfect. While it has split most of the polygons correctly, some are fragmented/shards. Here's screen shots of the input and output (via mapshaper.org):

Input:

Output:

The distortion was expected from the crs conversion, and looks ok on mapbox's background layers. But some of the perfectly-aligned hexagons (from the input) which were split along the anti-meridian have fragments / shards missing, e.g. the highlighted polygon at the bottom of the output.

Does anyone know how to avoid or fix this?

Answer 1

It seems that geopandas is the tool for the job in this scenario.

I attempted using transform_geom from Fiona as suggested here but the results were the same as what you show above.

But, since geopandas uses a different method to transform the geometries, it seems to do the trick:

import geopandas

gdf = geopandas.read_file("input.geojson")
gdf_wgs84 = gdf.to_crs("epsg:4326")
gdf_wgs84_cut = gdf_wgs84.set_geometry(gdf_wgs84.geometry.apply(idl_resolve))
gdf_wgs84_cut.to_file("output.geojson", driver="GeoJSON")

Note idl_resolve from https://github.com/Toblerity/Shapely/pull/95/files:

from shapely import geometry


def shift(geom):
    """
    Reads every point in every component of input geometry, and performs the following change:
        if the longitude coordinate is <0, adds 360 to it.
        if the longitude coordinate is >180, subtracts 360 from it.
    Useful for shifting between 0 and 180 centric map
    """

    if geom.is_empty:
        return geom

    if geom.has_z:
        num_dim = 3
    else:
        num_dim = 2

    def shift_pts(pts):
        """Internal function to perform shift of individual points"""
        if num_dim == 2:
            for x, y in pts:
                if x < 0:
                    x += 360
                elif x > 180:
                    x -= 360
                yield (x, y)
        elif num_dim == 3:
            for x, y, z in pts:
                if x < 0:
                    x += 360
                elif x > 180:
                    x -= 360
                yield (x, y, z)

    # Determine the geometry type to call appropriate handler
    if geom.type in ('Point', 'LineString'):
        return type(geom)(list(shift_pts(geom.coords)))
    elif geom.type == 'Polygon':
        ring = geom.exterior
        shell = type(ring)(list(shift_pts(ring.coords)))
        holes = list(geom.interiors)
        for pos, ring in enumerate(holes):
            holes[pos] = type(ring)(list(shift_pts(ring.coords)))
        return type(geom)(shell, holes)
    elif geom.type.startswith('Multi') or geom.type == 'GeometryCollection':
        # Recursive call to shift all components
        return type(geom)([shift(part)
                           for part in geom.geoms])
    else:
        raise ValueError('Type %r not supported' % geom.type)


def idl_resolve(geom, buffer_width=0.0000001):
    """
    Identifies when an intersection is present with -180/180 international date line and corrects it
    Geometry is shifted to 180 centric map and intersection is checked against a line defined as [(180, -90), (180,90)]
    If intersection is identified then the line is buffered by given amount (decimal degrees) and the difference
    between input geometry and buffer result is returned
    If no intersection is identified the passed in geometry is returned
    """

    intersecting_line = geometry.LineString(((180, -90), (180, 90)))

    shifted_geom = shift(geom)

    if shifted_geom.intersects(intersecting_line):
        buffered_line = intersecting_line.buffer(buffer_width)
        difference_geom = shifted_geom.difference(buffered_line)
        geom = shift(difference_geom)

    return geom

Environment information:

>>> import geopandas; geopandas.show_versions()

SYSTEM INFO
-----------
python     : 3.8.0 | packaged by conda-forge | (default, Nov 22 2019, 19:11:38)  [GCC 7.3.0]
executable : ../miniconda/envs/cart/bin/python
machine    : Linux-4.15.0-72-generic-x86_64-with-glibc2.10

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.8.0
GEOS lib   : ../miniconda/envs/cart/lib/libgeos_c.so
GDAL       : 3.0.2
GDAL data dir: ../miniconda/envs/cart/share/gdal
PROJ       : 6.2.1
PROJ data dir: ../miniconda/envs/cart/share/proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.6.2
pandas     : 0.25.3
fiona      : 1.8.13
numpy      : 1.17.3
shapely    : 1.6.4.post2
rtree      : 0.9.3
pyproj     : 2.4.2.post1
matplotlib : 3.1.2
mapclassify: None
pysal      : None
geopy      : None
psycopg2   : None