Repository: Toblerity/Shapely (python)

Issue:
How to use strtree to find an element nearest to it but not its own?
## Expected behavior and actual behavior.


How to use strtree to find an element nearest to it but not its own?My code is as follows, what I expect to find is (1,1), not (0,0).
If I want to start from (0,0), find the nearest point (1,1), and then find the second nearest point (2,2), how do I operate?

## Steps to reproduce the problem.

from shapely.strtree import STRtree
from shapely.geometry import Point

points=[Point(0,0),Point(1,1),Point(2,2)]

tree=STRtree(points)

pos=tree.nearest(Point(0,0))
print(pos)

#POINT (0 0)

## Operating system

windows 10    20H2
python3.7.0 64-bit (conda)

## Shapely version and provenance

Shapely-1.7.1 from pip install


Pull request description:
Optionally exclude input geom from results
Resolves #1106

Patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index a22880a..8b12c78 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -10,6 +10,8 @@ Shapely 1.8 will support only Python versions >= 3.6.
 
 New features:
 
+- The STRtree nearest*() methods now take an optional argument that
+  specifies exclusion of the input geometry from results (#1115).
 - A GeometryTypeError has been added to shapely.errors and is consistently
   raised instead of TypeError or ValueError as in version 1.7. For backwards
   compatibility, the new exception will derive from TypeError and Value error
@@ -23,6 +25,8 @@ New features:
 
 Bug fixes:
 
+- We no longer attempt to load libc to get the free function on Linux, but get
+  it from the global symbol table.
 - GEOS error messages printed when GEOS_getCoordSeq() is passed an empty
   geometry are avoided by never passing an empty geometry (#1134).
 - Python's builtin super() is now used only as described in PEP 3135 (#1109).
diff --git a/_vendor/packaging/specifiers.py b/_vendor/packaging/specifiers.py
index 902f6de..a538a62 100644
--- a/_vendor/packaging/specifiers.py
+++ b/_vendor/packaging/specifiers.py
@@ -503,7 +503,7 @@ class Specifier(_IndividualSpecifier):
                 return False
 
         # Ensure that we do not allow a local version of the version mentioned
-        # in the specifier, which is techincally greater than, to match.
+        # in the specifier, which is technically greater than, to match.
         if prospective.local is not None:
             if Version(prospective.base_version) == Version(spec.base_version):
                 return False
@@ -589,7 +589,7 @@ def _pad_version(left, right):
 class SpecifierSet(BaseSpecifier):
 
     def __init__(self, specifiers="", prereleases=None):
-        # Split on , to break each indidivual specifier into it's own item, and
+        # Split on , to break each individual specifier into it's own item, and
         # strip each item to remove leading/trailing whitespace.
         specifiers = [s.strip() for s in specifiers.split(",") if s.strip()]
 
diff --git a/docs/design.rst b/docs/design.rst
index 0c768f0..9a6ccf4 100644
--- a/docs/design.rst
+++ b/docs/design.rst
@@ -3,7 +3,7 @@ Design Notes
 ============
 
 Shapely provides classes that implement, more or less, the interfaces in the
-OGC's simple feature acess specification [1]_. The classes are defined in
+OGC's simple feature access specification [1]_. The classes are defined in
 similarly named modules under ``shapely.geometry``: ``Point`` is in
 ``shapely.geometry.point``, ``MultiPolygon`` is in
 ``shapely.geometry.multipolygon``. These classes derive from
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index 1b128ca..1ad0d07 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -101,7 +101,7 @@ def deserialize_wkb(data):
 def geos_geom_from_py(ob, create_func=None):
     """Helper function for geos_*_from_py functions in each geom type.
 
-    If a create_func is specified the coodinate sequence is cloned and a new
+    If a create_func is specified the coordinate sequence is cloned and a new
     geometry is created with it, otherwise the geometry is cloned directly.
     This behaviour is useful for converting between LineString and LinearRing
     objects.
diff --git a/shapely/geos.py b/shapely/geos.py
index dea6147..fe73419 100644
--- a/shapely/geos.py
+++ b/shapely/geos.py
@@ -97,11 +97,12 @@ if sys.platform.startswith('linux'):
             'libgeos_c.so',
         ]
         _lgeos = load_dll('libgeos_c', fallbacks=alt_paths)
-    # Necessary for environments with only libc.musl
-    c_alt_paths = [
-        'libc.musl-x86_64.so.1'
-    ]
-    free = load_dll('c', fallbacks=c_alt_paths).free
+
+    # ctypes.CDLL(None) internally calls dlopen(NULL), and as the dlopen
+    # manpage says, "If filename is NULL, then the returned handle is for the
+    # main program". This way we can let the linker do the work to figure out
+    # which libc Python is actually using.
+    free = CDLL(None).free
     free.argtypes = [c_void_p]
     free.restype = None
 
@@ -147,10 +148,6 @@ elif sys.platform == 'darwin':
             ]
         _lgeos = load_dll('geos_c', fallbacks=alt_paths)
 
-    # ctypes.CDLL(None) internally calls dlopen(NULL), and as the dlopen
-    # manpage says, "If filename is NULL, then the returned handle is for the
-    # main program". This way we can let the linker do the work to figure out
-    # which libc Python is actually using.
     free = CDLL(None).free
     free.argtypes = [c_void_p]
     free.restype = None
@@ -187,12 +184,13 @@ elif sys.platform == 'win32':
 
 elif sys.platform == 'sunos5':
     _lgeos = load_dll('geos_c', fallbacks=['libgeos_c.so.1', 'libgeos_c.so'])
-    free = CDLL('libc.so.1').free
+    free.restype = None
     free.argtypes = [c_void_p]
     free.restype = None
+
 else:  # other *nix systems
     _lgeos = load_dll('geos_c', fallbacks=['libgeos_c.so.1', 'libgeos_c.so'])
-    free = load_dll('c', fallbacks=['libc.so.6']).free
+    free = CDLL(None).free
     free.argtypes = [c_void_p]
     free.restype = None
 
diff --git a/shapely/ops.py b/shapely/ops.py
index 720d7c3..ba21ab7 100644
--- a/shapely/ops.py
+++ b/shapely/ops.py
@@ -145,7 +145,7 @@ class CollectionOperator:
         """Returns the union of a sequence of geometries
 
         This method replaces :meth:`cascaded_union` as the
-        prefered method for dissolving many polygons.
+        preferred method for dissolving many polygons.
         """
         try:
             if isinstance(geoms, BaseMultipartGeometry):
diff --git a/shapely/speedups/__init__.py b/shapely/speedups/__init__.py
index 1bbb04c..88ec30e 100644
--- a/shapely/speedups/__init__.py
+++ b/shapely/speedups/__init__.py
@@ -33,7 +33,7 @@ def enable():
     """Enable Cython speedups
 
     The shapely.speedups module contains performance enhancements written in C.
-    They are automaticaly installed when Python has access to a compiler and
+    They are automatically installed when Python has access to a compiler and
     GEOS development headers during installation, and are enabled by default.
 
     You can check if speedups are installed with the `available` attribute, and
diff --git a/shapely/strtree.py b/shapely/strtree.py
index 75b0f2c..ac1008f 100644
--- a/shapely/strtree.py
+++ b/shapely/strtree.py
@@ -21,6 +21,7 @@ References
 import ctypes
 import logging
 from typing import Any, ItemsView, Iterable, Iterator, Sequence, Tuple, Union
+import sys
 from warnings import warn
 
 from shapely.errors import ShapelyDeprecationWarning
@@ -249,7 +250,9 @@ class STRtree:
         """
         return self.query_geoms(geom)
 
-    def nearest_item(self, geom: BaseGeometry) -> Union[Any, None]:
+    def nearest_item(
+        self, geom: BaseGeometry, exclusive: bool = False
+    ) -> Union[Any, None]:
         """Query the tree for the node nearest to geom and get the item
         stored in the node.
 
@@ -259,6 +262,9 @@ class STRtree:
         ----------
         geom : geometry object
             The query geometry.
+        exclusive : bool, optional
+            Whether to exclude the item corresponding to the given geom
+            from results or not.  Default: False.
 
         Returns
         -------
@@ -289,10 +295,14 @@ class STRtree:
 
         def callback(item1, item2, distance, userdata):
             try:
+                callback_userdata = ctypes.cast(userdata, ctypes.py_object).value
                 idx = ctypes.cast(item1, ctypes.py_object).value
                 geom2 = ctypes.cast(item2, ctypes.py_object).value
                 dist = ctypes.cast(distance, ctypes.POINTER(ctypes.c_double))
-                lgeos.GEOSDistance(self._rev[idx]._geom, geom2._geom, dist)
+                if callback_userdata["exclusive"] and self._rev[idx].equals(geom2):
+                    dist[0] = sys.float_info.max
+                else:
+                    lgeos.GEOSDistance(self._rev[idx]._geom, geom2._geom, dist)
                 return 1
             except Exception:
                 log.exception("Caught exception")
@@ -303,12 +313,14 @@ class STRtree:
             ctypes.py_object(geom),
             envelope._geom,
             lgeos.GEOSDistanceCallback(callback),
-            None,
+            ctypes.py_object({"exclusive": exclusive}),
         )
         result = ctypes.cast(item, ctypes.py_object).value
         return result
 
-    def nearest_geom(self, geom: BaseGeometry) -> Union[BaseGeometry, None]:
+    def nearest_geom(
+        self, geom: BaseGeometry, exclusive: bool = False
+    ) -> Union[BaseGeometry, None]:
         """Query the tree for the node nearest to geom and get the
         geometry corresponding to the item stored in the node.
 
@@ -316,6 +328,9 @@ class STRtree:
         ----------
         geom : geometry object
             The query geometry.
+        exclusive : bool, optional
+            Whether to exclude the given geom from results or not.
+            Default: False.
 
         Returns
         -------
@@ -325,13 +340,15 @@ class STRtree:
         version 2.0.
 
         """
-        item = self.nearest_item(geom)
+        item = self.nearest_item(geom, exclusive=exclusive)
         if item is None:
             return None
         else:
             return self._rev[item]
 
-    def nearest(self, geom: BaseGeometry) -> Union[BaseGeometry, None]:
+    def nearest(
+        self, geom: BaseGeometry, exclusive: bool = False
+    ) -> Union[BaseGeometry, None]:
         """Query the tree for the node nearest to geom and get the
         geometry corresponding to the item stored in the node.
 
@@ -342,6 +359,9 @@ class STRtree:
         ----------
         geom : geometry object
             The query geometry.
+        exclusive : bool, optional
+            Whether to exclude the given geom from results or not.
+            Default: False.
 
         Returns
         -------
@@ -351,4 +371,4 @@ class STRtree:
         version 2.0.
 
         """
-        return self.nearest_geom(geom)
+        return self.nearest_geom(geom, exclusive=exclusive)


Test patch:
diff --git a/tests/test_strtree.py b/tests/test_strtree.py
index 04d5742..55b61c4 100644
--- a/tests/test_strtree.py
+++ b/tests/test_strtree.py
@@ -177,3 +177,22 @@ def test_nearest_items(geoms, items):
     with pytest.warns(ShapelyDeprecationWarning):
         tree = STRtree(geoms, items)
     assert tree.nearest_item(None) is None
+
+
+@pytest.mark.skipif(geos_version < (3, 6, 0), reason="GEOS 3.6.0 required")
+@pytest.mark.parametrize(
+    "geoms",
+    [
+        [
+            Point(0, 0.5),
+            Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+            Polygon([(0, 2), (1, 2), (1, 3), (0, 3)]),
+        ]
+    ],
+)
+@pytest.mark.parametrize("items", [list(range(1, 4)), list("abc")])
+@pytest.mark.parametrize("query_geom", [Point(0, 0.5)])
+def test_nearest_item_exclusive(geoms, items, query_geom):
+    with pytest.warns(ShapelyDeprecationWarning):
+        tree = STRtree(geoms, items)
+    assert tree.nearest_item(query_geom, exclusive=True) != items[0]