Adapt set_precision API and fix remaining GEOS 3.10 tests (#410)

.github/workflows/test-pip.yml

@@ -142,10 +142,10 @@ jobs:
         run: python setup.py build_ext --inplace
 
       - name: Run tests
-        continue-on-error: ${{ matrix.geos == 'main' || matrix.geos == '3.10.0' }}
+        continue-on-error: ${{ matrix.geos == 'main' }}
         run: pytest pygeos
 
       # Only run doctests on 1 runner (because of typographic differences in doctest results)
       - name: Run doctests
         run: pytest --doctest-modules pygeos --ignore=pygeos/tests
-        if: ${{ matrix.os == 'ubuntu-latest' && matrix.python == '3.8' }}
+        if: ${{ matrix.os == 'ubuntu-latest' && matrix.python == '3.10.0' }}

CHANGELOG.rst

@@ -17,9 +17,17 @@ Version 0.11 (unreleased)
 
 **API Changes**
 
+* The default behaviour of ``pygeos.set_precision`` is now to always return valid geometries.
+  Before, the default was ``preserve_topology=False`` which caused confusion because
+  it mapped to GEOS_PREC_NO_TOPO (the new 'pointwise').
+  At the same time, GEOS < 3.10 implementation was not entirely correct so that some geometries
+  did and some did not preserve topology with this mode. Now, the new ``mode`` argument controls
+  the behaviour and the ``preserve_topology`` argument is deprecated (#410).
 * When constructing a linearring through ``pygeos.linearrings`` or a polygon through 
   ``pygeos.polygons`` now a ``ValueError`` is raised (instead of a ``GEOSException``)
   if the ring contains less than 4 coordinates including ring closure (#378).
+* Removed deprecated ``normalize`` keyword argument in ``pygeos.line_locate_point`` and
+  ``pygeos.line_interpolate_point`` (#410).
 
 **Added GEOS functions**
 

pygeos/constructive.py

@@ -131,7 +131,7 @@ def buffer(
     >>> buffer(Geometry("POINT (10 10)"), 2, quadsegs=1)
     <pygeos.Geometry POLYGON ((12 10, 10 8, 8 10, 10 12, 12 10))>
     >>> buffer(Geometry("POINT (10 10)"), 2, quadsegs=2)
-    <pygeos.Geometry POLYGON ((12 10, 11.4 8.59, 10 8, 8.59 8.59, 8 10, 8.59 11....>
+    <pygeos.Geometry POLYGON ((12 10, 11.414 8.586, 10 8, 8.586 8.586, 8 10, 8.5...>
     >>> buffer(Geometry("POINT (10 10)"), -2, quadsegs=1)
     <pygeos.Geometry POLYGON EMPTY>
     >>> line = Geometry("LINESTRING (10 10, 20 10)")
@@ -147,7 +147,7 @@ def buffer(
     >>> buffer(line2, 2, cap_style="flat", join_style="mitre")
     <pygeos.Geometry POLYGON ((18 12, 18 20, 22 20, 22 8, 10 8, 10 12, 18 12))>
     >>> buffer(line2, 2, cap_style="flat", join_style="mitre", mitre_limit=1)
-    <pygeos.Geometry POLYGON ((18 12, 18 20, 22 20, 21.8 9, 21 8.17, 10 8, 10 12...>
+    <pygeos.Geometry POLYGON ((18 12, 18 20, 22 20, 21.828 9, 21 8.172, 10 8, 10...>
     >>> square = Geometry("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
     >>> buffer(square, 2, join_style="mitre")
     <pygeos.Geometry POLYGON ((-2 -2, -2 12, 12 12, 12 -2, -2 -2))>
@@ -876,11 +876,11 @@ def minimum_bounding_circle(geometry, **kwargs):
     Examples
     --------
     >>> minimum_bounding_circle(Geometry("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))"))
-    <pygeos.Geometry POLYGON ((12.1 5, 11.9 3.62, 11.5 2.29, 10.9 1.07, 10 4.44e...>
+    <pygeos.Geometry POLYGON ((12.071 5, 11.935 3.621, 11.533 2.294, 10.879 1.07...>
     >>> minimum_bounding_circle(Geometry("LINESTRING (1 1, 10 10)"))
-    <pygeos.Geometry POLYGON ((11.9 5.5, 11.7 4.26, 11.4 3.06, 10.8 1.96, 10 1, ...>
+    <pygeos.Geometry POLYGON ((11.864 5.5, 11.742 4.258, 11.38 3.065, 10.791 1.9...>
     >>> minimum_bounding_circle(Geometry("MULTIPOINT (2 2, 4 2)"))
-    <pygeos.Geometry POLYGON ((4 2, 3.98 1.8, 3.92 1.62, 3.83 1.44, 3.71 1.29, 3...>
+    <pygeos.Geometry POLYGON ((4 2, 3.981 1.805, 3.924 1.617, 3.831 1.444, 3.707...>
     >>> minimum_bounding_circle(Geometry("POINT (0 1)"))
     <pygeos.Geometry POINT (0 1)>
     >>> minimum_bounding_circle(Geometry("GEOMETRYCOLLECTION EMPTY"))

pygeos/geometry.py

@@ -1,10 +1,12 @@
+import warnings
 from enum import IntEnum
 
 import numpy as np
 
 from . import Geometry  # NOQA
-from . import _geometry, lib
+from . import _geometry, geos_version, lib
 from .decorators import multithreading_enabled, requires_geos
+from .enum import ParamEnum
 
 __all__ = [
     "GeometryType",
@@ -686,9 +688,15 @@ def get_precision(geometry, **kwargs):
     return lib.get_precision(geometry, **kwargs)
 
 
+class SetPrecisionMode(ParamEnum):
+    valid_output = 0
+    pointwise = 1
+    keep_collapsed = 2
+
+
 @requires_geos("3.6.0")
 @multithreading_enabled
-def set_precision(geometry, grid_size, preserve_topology=False, **kwargs):
+def set_precision(geometry, grid_size, mode="valid_output", **kwargs):
     """Returns geometry with the precision set to a precision grid size.
 
     By default, geometries use double precision coordinates (grid_size = 0).
@@ -716,9 +724,25 @@ def set_precision(geometry, grid_size, preserve_topology=False, **kwargs):
         geometry if precision grid size was not previously set). If this
         value is more precise than input geometry, the input geometry will
         not be modified.
-    preserve_topology : bool, default False
-        If True, will attempt to preserve the topology of a geometry after
-        rounding coordinates.
+    mode :  {'valid_output', 'pointwise', 'keep_collapsed'}, default 'valid_output'
+        This parameter determines how to handle invalid output geometries. There are three modes:
+
+        1. `'valid_output'` (default):  The output is always valid. Collapsed geometry elements
+           (including both polygons and lines) are removed. Duplicate vertices are removed.
+        2. `'pointwise'`: Precision reduction is performed pointwise. Output geometry
+           may be invalid due to collapse or self-intersection. Duplicate vertices are not
+           removed. In GEOS this option is called NO_TOPO.
+
+           .. note::
+
+             'pointwise' mode requires at least GEOS 3.10. It is accepted in earlier versions,
+             but the results may be unexpected.
+        3. `'keep_collapsed'`: Like the default mode, except that collapsed linear geometry
+           elements are preserved. Collapsed polygonal input elements are removed. Duplicate
+           vertices are removed.
+    preserve_topology : bool, optional
+        .. deprecated:: 0.11
+          This parameter is ignored. Use ``mode`` instead.
     **kwargs
         For other keyword-only arguments, see the
         `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
@@ -735,11 +759,33 @@ def set_precision(geometry, grid_size, preserve_topology=False, **kwargs):
     <pygeos.Geometry POINT Z (1 1 0.9)>
     >>> set_precision(Geometry("LINESTRING (0 0, 0 0.1, 0 1, 1 1)"), 1.0)
     <pygeos.Geometry LINESTRING (0 0, 0 1, 1 1)>
+    >>> set_precision(Geometry("LINESTRING (0 0, 0 0.1, 0.1 0.1)"), 1.0, mode="valid_output")
+    <pygeos.Geometry LINESTRING Z EMPTY>
+    >>> set_precision(Geometry("LINESTRING (0 0, 0 0.1, 0.1 0.1)"), 1.0, mode="pointwise")
+    <pygeos.Geometry LINESTRING (0 0, 0 0, 0 0)>
+    >>> set_precision(Geometry("LINESTRING (0 0, 0 0.1, 0.1 0.1)"), 1.0, mode="keep_collapsed")
+    <pygeos.Geometry LINESTRING (0 0, 0 0)>
     >>> set_precision(None, 1.0) is None
     True
     """
-
-    return lib.set_precision(geometry, grid_size, preserve_topology, **kwargs)
+    if isinstance(mode, str):
+        mode = SetPrecisionMode.get_value(mode)
+    elif not np.isscalar(mode):
+        raise TypeError("mode only accepts scalar values")
+    if "preserve_topology" in kwargs:
+        warnings.warn(
+            "preserve_topology is deprecated (ignored), use 'mode' instead",
+            UserWarning,
+            stacklevel=2,
+        )
+        del kwargs["preserve_topology"]
+    if mode == SetPrecisionMode.pointwise and geos_version < (3, 10, 0):
+        warnings.warn(
+            "'pointwise' is only supported for GEOS 3.10",
+            UserWarning,
+            stacklevel=2,
+        )
+    return lib.set_precision(geometry, grid_size, np.intc(mode), **kwargs)
 
 
 @multithreading_enabled

pygeos/linear.py

@@ -1,5 +1,3 @@
-from warnings import warn
-
 from . import Geometry  # NOQA
 from . import lib
 from .decorators import multithreading_enabled
@@ -47,9 +45,6 @@ def line_interpolate_point(line, distance, normalized=False, **kwargs):
     >>> line_interpolate_point(Geometry("LINESTRING EMPTY"), 1)
     <pygeos.Geometry POINT EMPTY>
     """
-    if "normalize" in kwargs:
-        warn("argument 'normalize' is deprecated; use 'normalized'", DeprecationWarning)
-        normalized = kwargs.pop("normalize")
     if normalized:
         return lib.line_interpolate_point_normalized(line, distance)
     else:
@@ -86,9 +81,6 @@ def line_locate_point(line, other, normalized=False, **kwargs):
     >>> line_locate_point(Geometry("LINESTRING EMPTY"), Geometry("POINT(4 4)"))
     nan
     """
-    if "normalize" in kwargs:
-        warn("argument 'normalize' is deprecated; use 'normalized'", DeprecationWarning)
-        normalized = kwargs.pop("normalize")
     if normalized:
         return lib.line_locate_point_normalized(line, other)
     else:

pygeos/predicates.py

@@ -408,7 +408,7 @@ def is_valid_reason(geometry, **kwargs):
     >>> is_valid_reason(Geometry("LINESTRING(0 0, 1 1)"))
     'Valid Geometry'
     >>> is_valid_reason(Geometry("POLYGON((0 0, 1 1, 1 2, 1 1, 0 0))"))
-    'Self-intersection[0 0]'
+    'Ring Self-intersection[1 1]'
     >>> is_valid_reason(None) is None
     True
     """