Repository: mattijn/topojson

Issue:
Bug: polygons that entirely fill islands in another polygon are often not dedupped
Because they aren't properly dedupped, this situation is created after simplify:
![image](https://user-images.githubusercontent.com/1565789/187953959-8e45dc49-bc0c-4e80-aa57-873f6c26e667.png)

I'm working on a fix...

Discussion:

Patch:
diff --git a/topojson/core/cut.py b/topojson/core/cut.py
index 66b9b0a..ea901c0 100644
--- a/topojson/core/cut.py
+++ b/topojson/core/cut.py
@@ -117,7 +117,9 @@ class Cut(Join):
             slist = []
 
             # create dict with original geometry type per linestring
-            linestring_object_types = self._get_linestring_types(data)
+            linestring_object_types = self._get_linestring_types(
+                objects=data["objects"], bookkeeping_geoms=data["bookkeeping_geoms"],
+            )
 
             # junctions are only existing in coordinates of linestring
             if self.options.shared_coords:
@@ -173,7 +175,28 @@ class Cut(Join):
 
         return data
 
-    def _get_linestring_types(self, data) -> dict:
+    def _get_linestring_types(
+        self, objects, bookkeeping_geoms, bookkeeping_linestrings=None
+    ) -> dict:
+        """
+        Returns the original geometry type for each linestring.
+
+        Parameters
+        ----------
+        objects : list
+            list of original objects that contains a list of arcs for each geometry.
+        bookkeeping_geoms : list
+            list of arc-linestrings for each arc in the objects list.
+        bookkeeping_linestrings : numpy array, optional
+            array with for each arc-linestring the corresponding linestrings that were
+            constructed after splitting them on junctions. Defaults to None.
+
+        Returns
+        -------
+        dict :
+            dict with for each data["linestrings"] index the geometry type of
+            the object the linestring originated from.
+        """
         # create dict with original geometry type per linestring
         def recurse_geometries(object):
             # If object is not a list, make it a list to be able to loop
@@ -190,14 +213,20 @@ class Cut(Join):
                     # For non-Point geometries, loop over arcs
                     for arc_id in object_child["arcs"]:
                         # Find the linestrings for the arc via bookkeeping_geoms
-                        for linestring_id in data["bookkeeping_geoms"][arc_id]:
-                            linestring_object_types[linestring_id] = (
-                                object_child["type"]
-                            )
+                        for arc_line_id in bookkeeping_geoms[arc_id]:
+                            if bookkeeping_linestrings is None:
+                                arc_lines = [arc_line_id]
+                            else:
+                                arc_lines = bookkeeping_linestrings[arc_line_id]
+                            for linestring_id in arc_lines:
+                                if linestring_id >= 0:
+                                    linestring_object_types[linestring_id] = (
+                                        object_child["type"]
+                                    )
 
         linestring_object_types = {}
-        for object_key in data["objects"]:
-            recurse_geometries(data["objects"][object_key])
+        for object_key in objects:
+            recurse_geometries(objects[object_key])
 
         return linestring_object_types
 
diff --git a/topojson/ops.py b/topojson/ops.py
index 77d613e..85c7815 100644
--- a/topojson/ops.py
+++ b/topojson/ops.py
@@ -903,7 +903,7 @@ def cart(arr):
 def find_duplicates(segments_list, type="array"):
     """
     Function for solely detecting and recording duplicate LineStrings. The function
-    converts sorts the coordinates of each linestring and gets the hash. Using the
+    converts and sorts the coordinates of each linestring and gets the hash. Using the
     hashes it can quickly detect duplicates and return the indices.
 
     Parameters
@@ -915,16 +915,26 @@ def find_duplicates(segments_list, type="array"):
 
     """
 
-    # get hash of sorted paths
+    # get hash of sorted linestring coordinates
     hash_segments = []
 
-    if type == "array":
-        for path in segments_list:
-            hash_segments.append(hash(bytes(np.sort(path, axis=0))))
+    if type != "array":
+        segments_list = [
+            np.array(list(linestring.coords)) for linestring in segments_list
+        ]
+
+    for coordinates in segments_list:
+        # If start and end points are the same, remove end point before sorting
+        # Rmark: check if it was originally a ring is not relevant, because lines with
+        # equal start and end point are no probem to be deduplicated with rings.
+        if np.array_equal(coordinates[0], coordinates[-1]):
+            coordinates = coordinates[0:-1]
+            coordinates = np.sort(coordinates, axis=0)
+            coordinates = np.append(coordinates[0:2], coordinates)
+        else:
+            coordinates = np.sort(coordinates, axis=0)
 
-    else:
-        for path in segments_list:
-            hash_segments.append(hash(tuple(sorted(path.coords))))
+        hash_segments.append(hash(bytes(coordinates)))
 
     hash_segments = np.array(hash_segments, dtype=np.int64)
 


Test patch:
diff --git a/tests/test_cut.py b/tests/test_cut.py
index 1ac902e..ad0fc04 100644
--- a/tests/test_cut.py
+++ b/tests/test_cut.py
@@ -295,8 +295,6 @@ def test_cut_shared_paths_reversed_rings_ABCA_ACBA_no_cuts():
 
 
 # cut rotated duplicate rings BCAB & ABCA have no cuts
-# changed the assertion of bookkeeping_duplicates, since the method rotated polygons are
-# not detected anymore in find_duplicate function when shapely is used for shared_paths.
 def test_cut_shared_paths_rotated_duplicates_rings_BCAB_ABCA_no_cuts():
     data = {
         "abca": {"type": "Polygon", "coordinates": [[[0, 0], [1, 0], [2, 1], [0, 0]]]},
@@ -305,7 +303,7 @@ def test_cut_shared_paths_rotated_duplicates_rings_BCAB_ABCA_no_cuts():
     topo = Cut(data, options={"shared_coords": False}).to_dict()
 
     assert len(topo["junctions"]) == 0
-    assert len(topo["bookkeeping_duplicates"]) == 0
+    assert len(topo["bookkeeping_duplicates"]) == 1
 
 
 # cut ring ABCA & line ABCA have no cuts
@@ -327,8 +325,6 @@ def test_cut_shared_paths_ring_ABCA_line_ABCA_no_cuts():
 
 
 # cut ring BCAB & line ABCA have no cuts
-# changed the assertion of bookkeeping_duplicates, since the method rotated polygons are
-# not detected anymore in find_duplicate function when shapely is used for shared_paths.
 def test_cut_shared_paths_ring_BCAB_line_ABCA_no_cuts():
     data = {
         "abcaLine": {
@@ -343,12 +339,10 @@ def test_cut_shared_paths_ring_BCAB_line_ABCA_no_cuts():
     topo = Cut(data, options={"shared_coords": False}).to_dict()
 
     assert len(topo["junctions"]) == 0
-    assert len(topo["bookkeeping_duplicates"]) == 0
+    assert len(topo["bookkeeping_duplicates"]) == 1
 
 
 # cut ring ABCA & line BCAB have no cuts
-# changed the assertion of bookkeeping_duplicates, since the method rotated polygons are
-# not detected anymore in find_duplicate function when shapely is used for shared_paths.
 def test_cut_shared_paths_ring_ABCA_line_BCAB_no_cuts():
     data = {
         "bcabLine": {
@@ -363,7 +357,7 @@ def test_cut_shared_paths_ring_ABCA_line_BCAB_no_cuts():
     topo = Cut(data, options={"shared_coords": False}).to_dict()
 
     assert len(topo["junctions"]) == 0
-    assert len(topo["bookkeeping_duplicates"]) == 0
+    assert len(topo["bookkeeping_duplicates"]) == 1
 
 
 # topoquantize can have low values, but prequantize cannot. this smells as a bug
diff --git a/tests/test_dedup.py b/tests/test_dedup.py
index 77c9afb..b7d2588 100644
--- a/tests/test_dedup.py
+++ b/tests/test_dedup.py
@@ -280,5 +280,5 @@ def test_dedup_merge_continuous_shared_path():
     )
     topo = Dedup(data, options={"shared_coords": False}).to_dict()
 
-    assert len(topo["linestrings"]) == 10
+    assert len(topo["linestrings"]) == 6
     assert len(topo["junctions"]) == 0
diff --git a/tests/test_topology.py b/tests/test_topology.py
index 44992c1..55a89d5 100644
--- a/tests/test_topology.py
+++ b/tests/test_topology.py
@@ -1,13 +1,15 @@
-import os
 import json
-import pytest
+import os
 
 import fiona
 import geojson
 import geopandas
+import geopandas.datasets
+import pytest
 from shapely import geometry, wkt
 
 import topojson
+import topojson.utils
 
 
 # this test was added since geometries of only linestrings resulted in a topojson
@@ -30,6 +32,62 @@ def test_topology_naturalearth_lowres_defaults():
     assert len(topo["objects"]) == 1
 
 
+# Test created for following issue: 
+#     Creating a topology for data without junctions and shared_coords=False, 
+#     prequantize=False gives error (https://github.com/mattijn/topojson/issues/181)
+# Test updated to also cover this issue:
+#     Polygons that entirely fill islands in another polygon are often not dedupped
+#     (https://github.com/mattijn/topojson/issues/183)
+@pytest.mark.parametrize("prequantize", [(True), (False)])
+def test_topology_polygon_filled_island_no_junctions(prequantize):
+    data = geopandas.GeoDataFrame(
+        {
+            "name": ["abcde_fghij", "jihgf"],
+            "geometry": [
+                geometry.Polygon(
+                    shell=[[0, 0], [3, 0], [3, 3], [0, 3], [0, 0]], 
+                    holes=[[[1, 1], [1, 2], [2, 2], [2, 1], [1, 1]]],
+                ),
+                geometry.Polygon([[2, 1], [2, 2], [1, 2], [1, 1], [2, 1]]),
+            ]
+        }
+    )
+    topo = topojson.Topology(data, prequantize=prequantize, shared_coords=False)
+    topo_gdf = topo.to_gdf()
+
+    assert len(topo.output["arcs"]) == 2
+    for index in range(len(topo_gdf)):
+        assert topo_gdf.geometry[index].equals(data["geometry"][index]), (
+            f"{topo_gdf.geometry[index].wkt} != {data['geometry'][index]}"
+        )
+
+# Test created for following issue: 
+#     Polygons that entirely fill islands in another polygon are often not dedupped
+#     (https://github.com/mattijn/topojson/issues/183)
+def test_topology_polygon_filled_island_with_junctions():
+    data = geopandas.GeoDataFrame(
+        {
+            "name": ["abcda_efghie", "fghief", "b__cb"],
+            "geometry": [
+                geometry.Polygon(
+                    shell=[[0, 0], [3, 0], [3, 3], [0, 3], [0, 0]], 
+                    holes=[[[1, 1], [1, 2], [2, 2], [2, 1], [1, 1]]],
+                ),
+                geometry.Polygon([[2, 1], [2, 2], [1, 2], [1, 1], [2, 1]]),
+                geometry.Polygon([[3, 0], [4, 0], [4, 3], [3, 3], [3, 0]]),
+            ]
+        }
+    )
+    topo = topojson.Topology(data, prequantize=False, shared_coords=False)
+    topo_gdf = topo.to_gdf()
+
+    assert len(topo.output["arcs"]) == 4
+    for index in range(len(topo_gdf)):
+        assert topo_gdf.geometry[index].equals(data["geometry"][index]), (
+            f"{topo_gdf.geometry[index].wkt} != {data['geometry'][index]}"
+        )
+
+
 # test winding order using TopoOptions object
 def test_topology_winding_order_TopoOptions():
     data = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
@@ -636,28 +694,3 @@ def test_topology_write_multiple_object_json_dict():
     topo_dict = topo.to_dict()
 
     assert len(topo_dict["objects"]) == 2
-
-
-@pytest.mark.parametrize(
-    "shared_coords, prequantize", 
-    [(True, True), (True, False), (False, True), (False, False)]
-)
-def test_topology_polygon_filled_island_no_junctions(shared_coords, prequantize):
-    data = geopandas.GeoDataFrame(
-        {
-            "name": ["abcde_fghij", "jihgf"],
-            "geometry": [
-                geometry.Polygon(
-                    shell=[[0, 0], [3, 0], [3, 3], [0, 3], [0, 0]], 
-                    holes=[[[1, 1], [1, 2], [2, 2], [2, 1], [1, 1]]],
-                ),
-                geometry.Polygon([[1, 1], [2, 1], [2, 2], [1, 2], [1, 1]]),
-            ]
-        }
-    )
-    topo = topojson.Topology(data, prequantize=prequantize, shared_coords=shared_coords)
-    topo_gdf = topo.to_gdf()
-
-    assert len(topo.output["arcs"]) == 2
-    assert topo_gdf.geometry[0] == data["geometry"][0]
-    assert topo_gdf.geometry[1] == data["geometry"][1]