Repository: geopandas/geopandas
Pull request #2963 (2023-07-19 16:23:10)

ENH: add GeoSeries.polygonize
Exposing `polygonize` on the GeoPandasBase level. In the end, I have opted to expose both `polygonize` and `polygonize_full` within a single method as it was easy enough and especially dangles can be quite useful even when you are not debugging the result.

Unlike shapely, we also include the noding step but not through `node` but using `unary_union` which has the same effect but is faster and more robust according to https://github.com/libgeos/geos/issues/877#issuecomment-1521056559.

martinfleis: @jorisvandenbossche good points, all in.

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index cf293cb5..caaebe88 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -10,6 +10,6 @@ repos:
           - id: black
             language_version: python3
     - repo: https://github.com/astral-sh/ruff-pre-commit
-      rev: "v0.1.13"
+      rev: "v0.4.4"
       hooks:
         - id: ruff
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2398eff3..9ccb39b3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -44,6 +44,8 @@ New methods:
 - Added `dwithin` method to check for a "distance within" predicate on
   GeoSeries/GeoDataFrame (#3153).
 - Added `to_geo_dict` method to generate GeoJSON-like dictionary from a GeoDataFrame (#3132).
+- Added `polygonize` method exposing both `polygonize` and `polygonize_full` from
+  shapely to GeoSeries/GeoDataframe (#2963).
 - Added `is_valid_reason` method from shapely to GeoSeries/GeoDataframe (#3176).
 
 New features and improvements:
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 1cb8312d..e09388ad 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -115,6 +115,7 @@ Constructive methods and attributes
    GeoSeries.minimum_clearance
    GeoSeries.minimum_rotated_rectangle
    GeoSeries.normalize
+   GeoSeries.polygonize
    GeoSeries.remove_repeated_points
    GeoSeries.reverse
    GeoSeries.sample_points
diff --git a/geopandas/base.py b/geopandas/base.py
index 949ec203..3b2bcd22 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2139,7 +2139,7 @@ GeometryCollection
 
         An object is said to contain `other` if at least one point of `other` lies in
         the interior and no points of `other` lie in the exterior of the object.
-        (Therefore, any given polygon does not contain its own boundary – there is not
+        (Therefore, any given polygon does not contain its own boundary - there is not
         any point that lies in the interior.)
         If either object is empty, this operation returns ``False``.
 
@@ -4583,7 +4583,7 @@ GeometryCollection
         return _binary_geo("shortest_line", self, other, align)
 
     def snap(self, other, tolerance, align=None):
-        """Snaps an input geometry to reference geometry’s vertices.
+        """Snaps an input geometry to reference geometry's vertices.
 
         Vertices of the first geometry are snapped to vertices of the second. geometry,
         returning a new geometry; the input geometries are not modified. The result
@@ -4915,7 +4915,7 @@ GeometryCollection
         each geometry.
 
         The algorithm (Douglas-Peucker) recursively splits the original line
-        into smaller parts and connects these parts’ endpoints
+        into smaller parts and connects these parts' endpoints
         by a straight line. Then, it removes all points whose distance
         to the straight line is smaller than `tolerance`. It does not
         move any points and it always preserves endpoints of
@@ -5880,6 +5880,93 @@ GeometryCollection
             ignore_index=True
         )
 
+    def polygonize(self, node=True, full=False):
+        """Creates polygons formed from the linework of a GeoSeries.
+
+        Polygonizes the GeoSeries that contain linework which represents the
+        edges of a planar graph. Any geometry type may be provided as input; only the
+        constituent lines and rings will be used to create the output polygons.
+
+        Lines or rings that when combined do not completely close a polygon will be
+        ignored. Duplicate segments are ignored.
+
+        Unless you know that the input GeoSeries represents a planar graph with a clean
+        topology (e.g. there is a node on both lines where they intersect), it is
+        recommended to use ``node=True`` which performs noding prior to polygonization.
+        Using ``node=False`` will provide performance benefits but may result in
+        incorrect polygons if the input is not of the proper topology.
+
+        When ``full=True``, the return value is a 4-tuple containing output polygons,
+        along with lines which could not be converted to polygons. The return value
+        consists of 4 elements or varying lenghts:
+
+        - GeoSeries of the valid polygons (same as with ``full=False``)
+        - GeoSeries of cut edges: edges connected on both ends but not part of
+          polygonal output
+        - GeoSeries of dangles: edges connected on one end but not part of polygonal
+          output
+        - GeoSeries of invalid rings: polygons that are formed but are not valid
+          (bowties, etc)
+
+        Parameters
+        ----------
+        node : bool, default True
+            Perform noding prior to polygonization, by default True.
+        full : bool, default False
+            Return the full output composed of a tuple of GeoSeries, by default False.
+
+        Returns
+        -------
+        GeoSeries | tuple(GeoSeries, GeoSeries, GeoSeries, GeoSeries)
+            GeoSeries with the polygons or a tuple of four GeoSeries as
+            ``(polygons, cuts, dangles, invalid)``
+
+        Examples
+        --------
+        >>> from shapely.geometry import LineString
+        >>> s = geopandas.GeoSeries([
+        ...     LineString([(0, 0), (1, 1)]),
+        ...     LineString([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
+        ...     LineString([(0.5, 0.2), (0.5, 0.8)]),
+        ... ])
+        >>> s.polygonize()
+        0        POLYGON ((0 0, 0.5 0.5, 1 1, 1 0, 0 0))
+        1    POLYGON ((0.5 0.5, 0 0, 0 1, 1 1, 0.5 0.5))
+        Name: polygons, dtype: geometry
+
+        >>> polygons, cuts, dangles, invalid = s.polygonize(full=True)
+
+        """
+        from .geoseries import GeoSeries
+
+        if node:
+            geometry_input = [self.geometry.union_all()]
+        else:
+            geometry_input = self.geometry.values
+
+        if full:
+            polygons, cuts, dangles, invalid = shapely.polygonize_full(geometry_input)
+
+            cuts = GeoSeries(cuts, crs=self.crs, name="cut_edges").explode(
+                ignore_index=True
+            )
+            dangles = GeoSeries(dangles, crs=self.crs, name="dangles").explode(
+                ignore_index=True
+            )
+            invalid = GeoSeries(invalid, crs=self.crs, name="invalid_rings").explode(
+                ignore_index=True
+            )
+            polygons = GeoSeries(polygons, crs=self.crs, name="polygons").explode(
+                ignore_index=True
+            )
+
+            return (polygons, cuts, dangles, invalid)
+
+        polygons = shapely.polygonize(geometry_input)
+        return GeoSeries(polygons, crs=self.crs, name="polygons").explode(
+            ignore_index=True
+        )
+
 
 def _get_index_for_parts(orig_idx, outer_idx, ignore_index, index_parts):
     """Helper to handle index when geometries get exploded to parts.
diff --git a/geopandas/plotting.py b/geopandas/plotting.py
index 678abaa4..eab1c0a3 100644
--- a/geopandas/plotting.py
+++ b/geopandas/plotting.py
@@ -855,7 +855,7 @@ def plot_dataframe(
             if scheme is not None:
                 categories = labels
             patches = []
-            for value, cat in enumerate(categories):
+            for i in range(len(categories)):
                 patches.append(
                     Line2D(
                         [0],
@@ -864,7 +864,7 @@ def plot_dataframe(
                         marker="o",
                         alpha=style_kwds.get("alpha", 1),
                         markersize=10,
-                        markerfacecolor=n_cmap.to_rgba(value),
+                        markerfacecolor=n_cmap.to_rgba(i),
                         markeredgewidth=0,
                     )
                 )
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index b006f61d..7b718545 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -52,7 +52,9 @@ class SpatialIndex:
         """
         return PREDICATES
 
-    def query(self, geometry, predicate=None, sort=False, distance=None):
+    def query(
+        self, geometry, predicate=None, sort=False, distance=None, output_format="tuple"
+    ):
         """
         Return the integer indices of all combinations of each input geometry
         and tree geometries where the bounding box of each input geometry
@@ -208,16 +210,36 @@ class SpatialIndex:
 
         indices = self._tree.query(geometry, predicate=predicate, **kwargs)
 
+        if output_format != "tuple":
+            sort = True
+
         if sort:
             if indices.ndim == 1:
-                return np.sort(indices)
+                indices = np.sort(indices)
             else:
                 # sort by first array (geometry) and then second (tree)
                 geo_idx, tree_idx = indices
                 sort_indexer = np.lexsort((tree_idx, geo_idx))
-                return np.vstack((geo_idx[sort_indexer], tree_idx[sort_indexer]))
+                indices = np.vstack((geo_idx[sort_indexer], tree_idx[sort_indexer]))
+
+        if output_format == "sparse":
+            from scipy.sparse import coo_array
+
+            return coo_array(
+                (np.ones(len(indices[0]), dtype=np.bool_), indices),
+                shape=(len(self.geometries), len(geometry)),
+                dtype=np.bool_,
+            )
+
+        if output_format == "dense":
+            dense = np.zeros((len(self.geometries), len(geometry)), dtype=bool)
+            dense[indices] = True
+            return dense
+
+        if output_format == "tuple":
+            return indices
 
-        return indices
+        raise ValueError("Invalid output_format: {}".format(output_format))
 
     @staticmethod
     def _as_geometry_array(geometry):
diff --git a/pyproject.toml b/pyproject.toml
index 2cd6d620..49603022 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -100,22 +100,23 @@ line-length = 88
 
 [tool.ruff]
 line-length = 88
+extend-exclude = ["doc/*", "benchmarks/*", "versioneer.py", "geopandas/_version.py"]
+
+[tool.ruff.lint]
 select = [
     # pyflakes
     "F",
     # pycodestyle
     "E",
     "W",
-    # flake8-2020
-    "YTT",
+    # pyupgrade
+    # "UP",
     # flake8-bugbear
     "B",
-    # flake8-quotes
-    "Q",
     # flake8-debugger
     "T10",
-    # flake8-gettext
-    "INT",
+    # flake8-simplify
+    # "SIM",
     # pylint
     "PLC",
     "PLE",
@@ -123,54 +124,30 @@ select = [
     "PLW",
     # misc lints
     "PIE",
-    # flake8-pyi
-    "PYI",
-    # tidy imports
-    "TID",
     # implicit string concatenation
     "ISC",
     # type-checking imports
     "TCH",
     # comprehensions
     "C4",
-    # pygrep-hooks
-    "PGH",
     # Ruff-specific rules
     "RUF",
+    # isort
+    # "I",
 ]
-target-version = "py39"
-ignore = [ # space before : (needed for how black formats slicing)
-    # "E203",  # not yet implemented
+
+ignore = [
+    ### Intentionally disabled
     # module level import not at top of file
     "E402",
     # do not assign a lambda expression, use a def
     "E731",
-    # line break before binary operator
-    # "W503",  # not yet implemented
-    # line break after binary operator
-    # "W504",  # not yet implemented
-    # controversial
+    # mutable-argument-default
     "B006",
-    # controversial
+    # unused-loop-control-variable
     "B007",
-    # controversial
-    "B008",
-    # setattr is used to side-step mypy
+    # get-attr-with-constant
     "B009",
-    # getattr is used to side-step mypy
-    "B010",
-    # tests use assert False
-    "B011",
-    # tests use comparisons but not their returned value
-    "B015",
-    # false positives
-    "B019",
-    # Loop control variable overrides iterable it iterates
-    "B020",
-    # Function definition does not bind loop variable
-    "B023",
-    # Functions defined inside a loop must not use variables redefined in the loop
-    # "B301",  # not yet implemented
     # Only works with python >=3.10
     "B905",
     # dict literals
@@ -183,50 +160,27 @@ ignore = [ # space before : (needed for how black formats slicing)
     "PLR0912",
     # Too many statements
     "PLR0915",
+    # Magic number
+    "PLR2004",
     # Redefined loop name
     "PLW2901",
     # Global statements are discouraged
     "PLW0603",
-    # Docstrings should not be included in stubs
-    "PYI021",
-    # Use `typing.NamedTuple` instead of `collections.namedtuple`
-    "PYI024",
-    # No builtin `eval()` allowed
-    "PGH001",
     # compare-to-empty-string
     "PLC1901",
-    # Use typing_extensions.TypeAlias for type aliases
-    # "PYI026",  # not yet implemented
-    # Use "collections.abc.*" instead of "typing.*" (PEP 585 syntax)
-    # "PYI027",  # not yet implemented
-    # while int | float can be shortened to float, the former is more explicit
-    # "PYI041",  # not yet implemented
-
-    # Additional checks that don't pass yet
+
+    ### Additional checks that don't pass yet
     # Useless statement
     "B018",
     # Within an except clause, raise exceptions with ...
     "B904",
-    # Magic number
-    "PLR2004",
     # Consider `elif` instead of `else` then `if` to remove indentation level
     "PLR5501",
-    # ambiguous-unicode-character-string
-    "RUF001",
-    # ambiguous-unicode-character-docstring
-    "RUF002",
-    # ambiguous-unicode-character-comment
-    "RUF003",
     # collection-literal-concatenation
     "RUF005",
-    # pairwise-over-zipped (>=PY310 only)
-    "RUF007",
-    # explicit-f-string-type-conversion
-    "RUF010",
     # Mutable class attributes should be annotated with `typing.ClassVar`,
     "RUF012"
 ]
-extend-exclude = ["doc/*", "benchmarks/*", "versioneer.py", "geopandas/_version.py"]
 
-[tool.ruff.per-file-ignores]
+[tool.ruff.lint.per-file-ignores]
 "geopandas/__init__.py" = ["F401"]


Test patch:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index d5b88bd7..9e551024 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -820,7 +820,7 @@ def test_setitem(item):
 
 def test_equality_ops():
     with pytest.raises(ValueError):
-        P[:5] == P[:7]
+        _ = P[:5] == P[:7]
 
     a1 = from_shapely([points[1], points[2], points[3]])
     a2 = from_shapely([points[1], points[0], points[3]])
diff --git a/geopandas/tests/test_decorator.py b/geopandas/tests/test_decorator.py
index efc3d37e..94b1dbee 100644
--- a/geopandas/tests/test_decorator.py
+++ b/geopandas/tests/test_decorator.py
@@ -10,7 +10,6 @@ def cumsum(whatever):
 
     It computes the cumulative {operation}.
     """
-    ...
 
 
 @doc(
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index a60f7857..79d99ac9 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -1300,7 +1300,7 @@ class TestConstructor:
             geometry="geometry",
         )
         check_geodataframe(gdf)
-        gdf.columns == ["geometry", "a"]
+        assert list(gdf.columns) == ["geometry", "a"]
 
         # with non-default index
         gdf = GeoDataFrame(
@@ -1310,7 +1310,7 @@ class TestConstructor:
             geometry="geometry",
         )
         check_geodataframe(gdf)
-        gdf.columns == ["geometry", "a"]
+        assert list(gdf.columns) == ["geometry", "a"]
 
     def test_preserve_series_name(self):
         geoms = [Point(1, 1), Point(2, 2), Point(3, 3)]
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 82550a17..ce9bf54c 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -174,6 +174,20 @@ class TestGeomMethods:
         self.l6 = LineString([(5, 5), (5, 100), (100, 5)])
         self.g12 = GeoSeries([self.l5])
         self.g13 = GeoSeries([self.l6])
+        self.lines = GeoSeries(
+            [
+                LineString([(0, 0), (1, 1)]),
+                LineString([(0, 0), (0, 1)]),
+                LineString([(0, 1), (1, 1)]),
+                LineString([(1, 1), (1, 0)]),
+                LineString([(1, 0), (0, 0)]),
+                LineString([(5, 5), (6, 6)]),
+                LineString([(0.5, -1), (0.5, 2)]),
+                Point(0, 0),
+            ],
+            crs=4326,
+            index=range(2, 10),
+        )
 
         self.l5 = LineString([(100, 0), (0, 0), (0, 100)])
         self.l6 = LineString([(5, 5), (5, 100), (100, 5)])
@@ -1928,6 +1942,61 @@ class TestGeomMethods:
         ):
             GeoSeries([self.l1]).offset_curve(Series([1], index=[99]))
 
+    def test_polygonize(self):
+        expected = GeoSeries.from_wkt(
+            [
+                "POLYGON ((0 0, 0.5 0.5, 0.5 0, 0 0))",
+                "POLYGON ((0.5 0.5, 0 0, 0 1, 0.5 1, 0.5 0.5))",
+                "POLYGON ((0.5 0.5, 1 1, 1 0, 0.5 0, 0.5 0.5))",
+                "POLYGON ((1 1, 0.5 0.5, 0.5 1, 1 1))",
+            ],
+            name="polygons",
+            crs=4326,
+        )
+
+        result = self.lines.polygonize()
+        assert_geoseries_equal(expected, result)
+        assert_index_equal(self.lines.index, Index(range(2, 10)))
+
+    def test_polygonize_no_node(self):
+        expected = GeoSeries.from_wkt(
+            ["POLYGON ((0 0, 1 1, 1 0, 0 0))", "POLYGON ((1 1, 0 0, 0 1, 1 1))"],
+            name="polygons",
+            crs=4326,
+        )
+        result = self.lines.polygonize(node=False)
+        assert_geoseries_equal(expected, result)
+        assert_index_equal(self.lines.index, Index(range(2, 10)))
+
+    def test_polygonize_full(self):
+        expected_poly = GeoSeries.from_wkt(
+            [
+                "POLYGON ((0 0, 0.5 0.5, 0.5 0, 0 0))",
+                "POLYGON ((0.5 0.5, 0 0, 0 1, 0.5 1, 0.5 0.5))",
+                "POLYGON ((0.5 0.5, 1 1, 1 0, 0.5 0, 0.5 0.5))",
+                "POLYGON ((1 1, 0.5 0.5, 0.5 1, 1 1))",
+            ],
+            name="polygons",
+            crs=4326,
+        )
+        expected_cuts = GeoSeries([], name="cut edges", crs=4326)
+        expected_dangles = GeoSeries.from_wkt(
+            [
+                "LINESTRING (5 5, 6 6)",
+                "LINESTRING (0.5 1, 0.5 2)",
+                "LINESTRING (0.5 -1, 0.5 0)",
+            ],
+            name="dangles",
+            crs=4326,
+        )
+        expected_invalid = GeoSeries([], name="invalid ring lines", crs=4326)
+        result = self.lines.polygonize(full=True)
+        assert_geoseries_equal(expected_poly, result[0])
+        assert_geoseries_equal(expected_cuts, result[1])
+        assert_geoseries_equal(expected_dangles, result[2])
+        assert_geoseries_equal(expected_invalid, result[3])
+        assert_index_equal(self.lines.index, Index(range(2, 10)))
+
     @pytest.mark.skipif(shapely.geos_version < (3, 11, 0), reason="requires GEOS>=3.11")
     @pytest.mark.parametrize(
         "geom,expected",