4

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:

enter image description here

Output:

enter image description here

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?

  • 2
    Please Edit the question to embed the images within the question. Links are far less likely to be viewed. – Vince Dec 15 '19 at 18:13
  • Thanks for fixing that @user2856 - I tried a few times and couldn't get the images to work when I was posting it. – Mike Honey Dec 16 '19 at 1:47
  • 1
    Are you able to share the input GeoJSON file? – snowman2 Dec 24 '19 at 2:03
  • In the original GeoJSON, is there any difference between the hexagon which gets incorrectly converted and the ones which are correctly converted? In particular: do they all have the same number of vertexes? Any duplicated, overlapping vertexes? – RafDouglas Dec 27 '19 at 17:22
  • 1
    @RafDouglas - they appear identical, mapshaper doesnt report any issues. Ref above I've shared my files. – Mike Honey Dec 29 '19 at 7:39
4
+100

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

enter image description here enter image description here

  • Thanks so much for coming up with that - it worked perfectly. I'm very happy to change to geopandas. – Mike Honey Jan 5 at 6:12

Your Answer

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.