Invalid GeoJSON

[Plantain]

My GeoJSON wasn't rendering properly in combination with geojson-vt, and after filing a bug with them, they identified that the GeoJSON from geojsoncontour was not valid.

mapbox/tippecanoe#580

The fundamental problem is that many of the outer polygon rings are marked in the GeoJSON as being interior rings.

o For Polygons with more than one of these rings, the first MUST be
the exterior ring, and any others MUST be interior rings. The
exterior ring bounds the surface, and the interior rings (if
present) bound holes within the surface.
Ignoring the GeoJSON spec and detecting inner/outer rings by winding rather than by ring order within the geometry gives the intended result:

Can this be fixed within geojsoncontour?

[bartromgens]

@Plantain Thank you for reporting this issue and sorry for my late reply.

From what I understand, this should be fixed in geojsoncontour.
@ucyo wrote the conversion to polygons. I am not sure how much work it is to fix this. I hope to find some time to have a look at it, but don't count on it.

Here is the related code:

geojsoncontour/geojsoncontour/contour.py

Line 68 in 79e3071

def contourf_to_geojson(contourf, geojson_filepath=None, min_angle_deg=None,

PRs are more than welcome!

Hi, we bumped into this as well. Here is a minimal example that demonstrates the problem using the validation feature from the shapely library:

import numpy as np
import matplotlib.pyplot as plt
import shapely.geometry
import json
import geojsoncontour

arr = np.array([
    [0, 1, 0],
    [1, 1, 1],
    [0, 1, 0],
])

contour = plt.contourf(arr)
geodata = geojsoncontour.contourf_to_geojson(contour, serialize=False)
for feature in geodata["features"]:
    print(shapely.geometry.asShape(feature["geometry"]).is_valid)

This code prints False for each feature except the last one.

In case it might help, my coworker managed to fix the GeoJSON afterwards using this code:

for n  in np.arange(0, len(geodata['features'])): 
    poly = shapely.geometry.asPolygon(geodata['features'][n]["geometry"]["coordinates"][0][0])
    for i in np.arange(0, len(geodata['features'][n]["geometry"]["coordinates"])): 
        for k in np.arange(1, len( geodata['features'][n]["geometry"]["coordinates"][i])):
            n_poly = shapely.geometry.asPolygon(geodata['features'][n]["geometry"]["coordinates"][i][k])
            if poly.contains(n_poly):
                poly = poly.difference(n_poly)
            else: 
                poly = poly.union(n_poly)
    geodata["features"][n]["geometry"] = poly.__geo_interface__ 

[MarySymons]

I have this problem also, it is not an issue when testing on geojson.io , that seems to know when polygons are inside others. But becomes apparent in google maps! Has there been any progress in fixing this? :)

I will try the above code as well, though I imagine with lots of data, this would not be too efficient. Will test it out!

[bartromgens]

@MarySymons Thanks for the update. There is no progress on this issue :( I would like to dive into it but cannot find the time right now.

@amarandon Thanks for the minimal example that demonstrates the problem. I hope to give it a try soon.

[MarySymons]

Update for you... I have used @amarandon's colleagues fix using shapely and it works for my problem too! Combining these things isn't particularly efficient as you can imagine, but it gets my desired outcome :)

[Plantain]

Any chance of taking a look at this @bartromgens ? I am using a very inefficient workaround, and so far I haven't been able to understand the code to come up with a better fix.

[ucyo]

Hi @Plantain,
original author of the polygon converter here. Next month I will have more spare time and will work on this 👍

[bartromgens]

I had a look at this issue and found the following issues with the current multipolygon representation in comparison to the geojson standard:

• A single contour level is currently represented as single Polygon of multiple LineStrings (linear rings) inside a MultiPolygon, instead of multiple Polygons
• geojson only considers the first LineString of a Polygon as exterior ring, all others are interior rings (holes)
• Hole coordinates should be defined clockwise

Things to check/fix:

• Independent areas (exterior rings) of a single level should be individual Polygons inside a MultiPolygon
  • This simply means moving the coordinate list one level up in the coordinates of the MultiPolygon. This works, but ignores holes, which become independent exterior rings.
• Holes should be defined as part of a single Polygon
  • The first linear ring of a Polygon is the exterior ring
  • How to know when a matplotlib Path is a hole?
• Holes should be defined clockwise
  • This seems to pose no real problem in most implementations if the two points above are fixed

Relevant geojson documentation:

A linear ring MUST follow the right-hand rule with respect to the
area it bounds, i.e., exterior rings are counterclockwise, and
holes are clockwise.

For Polygons with more than one of these rings, the first MUST be
the exterior ring, and any others MUST be interior rings. The
exterior ring bounds the surface, and the interior rings (if
present) bound holes within the surface.

• https://www.rfc-editor.org/rfc/rfc7946.html#section-3.1.6

[ucyo]

@bartromgens thank you for the research. It seems that „simply“ the ordering/direction of the Polygons are the problem. Matplotlib seems not to care about the specific order contrary to GeoJSON. But this seems reasonable since matplotlib does not support „holes“ in polygons as GeoJSON interprets different clockwise directions as an indicator for this. Anyway... soon I’ll have spare time and work on this 😺

[Plantain]

It seems the winding is fine, because using the --use-source-polygon-winding flag on tippecanoe makes it work with geojsoncontour just fine mapbox/tippecanoe#580 & mapbox/tippecanoe#581

[bartromgens]

I tried to fix the issues by simply making every matplotlib path its own Polygon inside the MultiPolygon (instead of a Polygon with multiple LineStrings).
In this case, http://geojson.io rendered the holes as filled areas while http://geojson.tools rendered the holes correct.

Likely because the winding is indeed correct (holes clockwise, exterior ccw), but geojson.io requires the holes to be part of a Polygon in which it is a hole, as by the standard.

The real problem is knowing when a path is a hole (check winding?) and to which exterior this belongs. If we have this information, it is fairly easy to create the correct MultiPolygon.

[Plantain]

I tried reordering the polygons manually through sorting by winding so the clockwise winding was always the exterior, but it still doesn't seem valid and I'm not sure why:
def ccw(a):
return (shapely.geometry.LinearRing(a).is_ccw)
polygon.coordinates[0] = sorted(polygon.coordinates[0], key=ccw)

[veloman-yunkan]

I think I have a fix for this issue. Need to debug one last problem observed by OP (@Plantain). Then I will submit a pull request.

[veloman-yunkan]

Pull Request #14 submitted