winding_order does not work for enclosed polygons

[mattijn]

winding_order does not work for enclosed polygons

[mattijn]

Thanks for raising the issue. This is technically not a bug in this package, but might be a feature request.

TLDR: it turns out that there is no winding-order enforcement in Fiona/OGR (used by geopandas) for polygons, currently used in the .to_gdf() function.

Your input data is following the left-hand rule (clockwise exterior, counterclockwise interior polygons), where this should be following the right-hand rule (counterclockwise exterior, clockwise interior polygons) to be conform GeoJSON standard.

import topojson as tp
from shapely import geometry
import json
import geopandas

gdf = geopandas.GeoDataFrame({
    "name": ["P1"],
    "geometry": [
        geometry.Polygon([(0, 0), (0, 3), (3, 3), (3, 0), (0, 0)], [[(1, 1), (2, 1), (2, 2), (1, 2), (1, 1)]])
    ]
})
print(json.dumps(gdf.__geo_interface__))

{"type": "FeatureCollection", "features": [{"id": "0", "type": "Feature", "properties": {"name": "P1"}, "geometry": {"type": "Polygon", "coordinates": [[[0.0, 0.0], [0.0, 3.0], [3.0, 3.0], [3.0, 0.0], [0.0, 0.0]], [[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]}, "bbox": [0.0, 0.0, 3.0, 3.0]}], "bbox": [0.0, 0.0, 3.0, 3.0]}

Copy+paste output in https://geojsonlint.com/ gives:

Invalid GeoJSON
Line 1: Polygons and MultiPolygons should follow the right-hand rule

TopoJSON standard requires left-hand rule (clockwise exterior, counterclockwise interior polygons), so all input polygons are
enforced (by default) to have a winding order clockwise for exterior and counterclockwise for interior polygons

topo = tp.Topology(gdf, prequantize=False)  # , winding_order='CW_CCW')
topo

Topology(
{'arcs': [[[0.0, 0.0], [0.0, 3.0], [3.0, 3.0], [3.0, 0.0], [0.0, 0.0]],
          [[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]],
 'bbox': (0.0, 0.0, 3.0, 3.0),
 'coordinates': [],
 'objects': {'data': {'geometries': [{'arcs': [[0], [1]],
                                      'properties': {'name': 'P1'},
                                      'type': 'Polygon'}],
                      'type': 'GeometryCollection'}},
 'type': 'Topology'}
)

The topo.to_gdf() function activates the serialize_as_geodataframe function (source code) where the TopoJSON object is directly read as virtual file by Fiona and mapped to a GeoDataFrame. And since it seems there is no winding order enforcement in the OGR model and thus directly parsed into the GeoDataFrame when reading the data:

topo.to_gdf().__geo_interface__

{'type': 'FeatureCollection',
 'features': [{'id': '0',
   'type': 'Feature',
   'properties': {'id': None, 'name': 'P1'},
   'geometry': {'type': 'Polygon',
    'coordinates': (((0.0, 0.0),
      (0.0, 3.0),
      (3.0, 3.0),
      (3.0, 0.0),
      (0.0, 0.0)),
     ((1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)))},
   'bbox': (0.0, 0.0, 3.0, 3.0)}],
 'bbox': (0.0, 0.0, 3.0, 3.0)}

The topo.to_geojson() activates my own implemented TopoJSON 2 GeoJSON parser (source code) which defaults to the correct winding order (only in master btw! NOT in latest release 1.0rc10) rotation for the GeoJSON standard (counterclockwise exterior, clockwise interior polygons).

print(topo.to_geojson(pretty=True))

{
    "type": "FeatureCollection",
    "features": [
        {
            "id": 0,
            "type": "Feature",
            "properties": {"name": "P1"},
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[0.0, 0.0], [3.0, 0.0], [3.0, 3.0], [0.0, 3.0], [0.0, 0.0]],
                    [[1.0, 1.0], [1.0, 2.0], [2.0, 2.0], [2.0, 1.0], [1.0, 1.0]]
                ]
            }
        }
    ]
}

(this can be tested by copy+pasting in https://geojsonlint.com/)

So this approach also could be used to go back to a GeoDataFrame with winding order of polygons following the GeoJSON standard (again, currently on master only):

geopandas.GeoDataFrame().from_features(json.loads(topo.to_geojson())['features']).__geo_interface__

{'type': 'FeatureCollection',
 'features': [{'id': '0',
   'type': 'Feature',
   'properties': {'name': 'P1'},
   'geometry': {'type': 'Polygon',
    'coordinates': (((0.0, 0.0),
      (3.0, 0.0),
      (3.0, 3.0),
      (0.0, 3.0),
      (0.0, 0.0)),
     ((1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)))},
   'bbox': (0.0, 0.0, 3.0, 3.0)}],
 'bbox': (0.0, 0.0, 3.0, 3.0)}

Question to myself (you and others can give suggestions): is it be better to replace the too loosely implemented approach by Fiona/OGR in the serialize_as_geodataframe function by inclusion of the strict enforced .to_geojson() function?

[mattijn]

One more thought:

If you do not store your Topology object as TopoJSON, but will make the roundtrip to a GeoDataFrame anyway, you also can enforce GeoJSON standard winding order when creating the Topology object. And since OGR don't touch the order, the to_gdf() output is as expected following the GeoJSON right-hand rule:

topo = tp.Topology(gdf, winding_order='CCW_CW')
topo.to_gdf().__geo_interface__

{'type': 'FeatureCollection',
 'features': [{'id': '0',
   'type': 'Feature',
   'properties': {'id': None, 'name': 'P1'},
   'geometry': {'type': 'Polygon',
    'coordinates': (((0.0, 0.0),
      (3.0, 0.0),
      (3.0, 3.0),
      (0.0, 3.0),
      (0.0, 0.0)),
     ((1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)))},
   'bbox': (0.0, 0.0, 3.0, 3.0)}],
 'bbox': (0.0, 0.0, 3.0, 3.0)}

[IvKor]

It works with single polygons. I think the issue appears when there is a small polygon filling the hole of a large polygon with a hole.

gdf = geopandas.GeoDataFrame({
    "name": ["P1", "P2"],
    "geometry": [
        geometry.Polygon([(0, 0), (0, 3), (3, 3), (3, 0), (0, 0)], [[(1, 1), (2, 1), (2, 2), (1, 2), (1, 1)]]),
        geometry.Polygon([(1, 1), (1, 2), (2, 2), (2, 1), (1, 1)])
    ]
})

topo = tp.Topology(gdf, prequantize=False)
print(topo)

Topology(
{'arcs': [[[0.0, 0.0], [0.0, 3.0], [3.0, 3.0], [3.0, 0.0], [0.0, 0.0]],
[[1.0, 1.0], [1.0, 2.0], [2.0, 2.0], [2.0, 1.0], [1.0, 1.0]]],
'bbox': (0.0, 0.0, 3.0, 3.0),
'coordinates': [],
'objects': {'data': {'geometries': [{'arcs': [[0], [1]],
'properties': {'name': 'P1'},
'type': 'Polygon'},
{'arcs': [[1]],
'properties': {'name': 'P2'},
'type': 'Polygon'}],
'type': 'GeometryCollection'}},
'type': 'Topology'}
)

It seems the order of arcs is what determines the polygon structure. The arcs [[1]] has CW direction and it is used for both the inner polygon in P1 and the exterior polygon in P2.

print(topo.to_gdf().__geo_interface__)

{'type': 'FeatureCollection', 'features':[{'id': '0', 'type': 'Feature', 'properties': {'id': None, 'name': 'P1'}, 'geometry': {'type': 'Polygon', 'coordinates': (((0.0, 0.0), (0.0, 3.0), (3.0, 3.0), (3.0, 0.0), (0.0, 0.0)), ((1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)))}, 'bbox': (0.0, 0.0, 3.0, 3.0)}, {'id': '1', 'type': 'Feature', 'properties': {'id': None, 'name': 'P2'}, 'geometry': {'type': 'Polygon', 'coordinates': (((1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)),)}, 'bbox': (1.0, 1.0, 2.0, 2.0)}], 'bbox': (0.0, 0.0, 3.0, 3.0)}

[mattijn]

Ok, than my second thought does not work.

And you really have to use the .to_geojson() function (currently master only) which correctly resolves to clockwise order for the interior polygon (P1) and counterclockwise for single polygon in P2:

print(topo.to_geojson(pretty=True))

{
    "type": "FeatureCollection",
    "features": [
        {
            "id": 0,
            "type": "Feature",
            "properties": {"name": "P1"},
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[0.0, 0.0], [3.0, 0.0], [3.0, 3.0], [0.0, 3.0], [0.0, 0.0]],
                    [[1.0, 1.0], [1.0, 2.0], [2.0, 2.0], [2.0, 1.0], [1.0, 1.0]]
                ]
            }
        },
        {
            "id": 1,
            "type": "Feature",
            "properties": {"name": "P2"},
            "geometry": {
                "type": "Polygon",
                "coordinates": [[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]
            }
        }
    ]
}

Given these type of confusing issues, I will probably just remove the Fiona/OGR approach to parse the Topology object into a GeoDataFrame.

Let me leave this issue open until I've released a new version. Thanks again for raising!

[mattijn]

Since #108 the .to_gdf() function does not use the Fiona/OGR approach anymore, but uses internally the to_geojson() route, which has become more robust as well.