Repository: Toblerity/Shapely
Pull request #857 (2020-03-01 02:35:40)

Strtree docs
I've update the docs for the STRtree in manual.rst. This should close these issues (#854, #845, #618) and one other pull request (https://github.com/Toblerity/Shapely/pull/623). Note that I have not done anything with the docstrings. I think I would like to work on that in a separate merge later. Maybe waiting to close those issues until that is done as well is a good idea. 

Also note that I have changed the work around for getting indexes to use an auxiliary dictionary keyed on geometry memory IDs instead of adding an index attribute to the geometries themselves (or a deep copy if you don't want to mutate.) 

Compare: 

``` python 
>>> from copy import deepcopy
>>> 
>>> points = [Point(i, i) for i in range(10)]
>>> tree = STRtree(points)
>>>
>>> points_numbered = []
>>> for i, point in enumerate(points):
>>>     point = deepcopy(point)  # let's not mutate the originals
>>>     point.i = i
>>>     points_numbered.append(point)
>>> tree = STRtree(points_numbered)
>>> result = tree.query(Point(2, 2).buffer(1.0))
>>> [(o.i, o.wkt) for o in result]
[(1, 'POINT (1 1)'), (2, 'POINT (2 2)'), (3, 'POINT (3 3)')]
```  
vs 

``` python 
>>> points = [Point(i, i) for i in range(10)]
>>> tree = STRtree(points)
>>>
>>> index_by_id = dict((id(pt), i) for i, pt in enumerate(points))
>>> [(index_by_id[id(pt)], pt.wkt) for pt in tree.query(Point(2,2).buffer(1.0))]
[(1, 'POINT (1 1)'), (2, 'POINT (2 2)'), (3, 'POINT (3 3)')]
```

This appears to be robust even in the case where the list grows, which you might expect to change memory location (id) of the geometries.  Observe:

``` python
>>> points = [Point(-1, -1)]
>>> original_id = id(points[0])
>>> for i in range(1000000):
>>>     points.append([Point(i, i)])
>>>     if id(points[0]) != original_id:
>>>         print(f"Element id changed from {original_id} to {id(points[0])}")
>>>         break
>>> original_id == id(points[0])
True
```

The RST rendering has been tested. Results are available here: https://jeremiah-england.github.io/shapely/build/manual.html#str-packed-r-tree  

Any review and feedback will be appreciated. 

coveralls: 
[![Coverage Status](https://coveralls.io/builds/29051545/badge)](https://coveralls.io/builds/29051545)

Coverage remained the same at 80.948% when pulling **a1b89b79c5d2dc8759660c88f2bb97659be5e9ff on Jeremiah-England:strtree-docs** into **cabea3a3a3aca6166dc1c2e1d0dfdcc5b6694987 on Toblerity:master**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/29051545/badge)](https://coveralls.io/builds/29051545)

Coverage remained the same at 80.948% when pulling **a1b89b79c5d2dc8759660c88f2bb97659be5e9ff on Jeremiah-England:strtree-docs** into **cabea3a3a3aca6166dc1c2e1d0dfdcc5b6694987 on Toblerity:master**.

Jeremiah-England: @sgillies I've fixed the problems you mentioned. Let me know if you find any others. Thanks. 

Patch:
diff --git a/docs/manual.rst b/docs/manual.rst
index a8d0918..48e6bcb 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -1677,8 +1677,8 @@ linestring feature (right).
   object on its right or its left side.
 
   Distance must be a positive float value. The side parameter may be 'left' or
-  'right'. Left and right is determined by following the direction of given 
-  geometric points of the LineString. The resolution of the offset around 
+  'right'. Left and right is determined by following the direction of given
+  geometric points of the LineString. The resolution of the offset around
   each vertex of the object is parameterized as in the buffer method.
 
   The join style is for outside corners between line segments. Accepted integer
@@ -2300,9 +2300,9 @@ The :func:`~shapely.ops.substring` function in `shapely.ops` returns a line segm
     Negative distance values are taken as measured in the reverse
     direction from the end of the geometry. Out-of-range index
     values are handled by clamping them to the valid range of values.
-    
+
     If the start distance equals the end distance, a point is being returned.
-    
+
     If the normalized arg is True, the distance will be interpreted as a
     fraction of the geometry's length.
 
@@ -2409,37 +2409,61 @@ STR-packed R-tree
 
 Shapely provides an interface to the query-only GEOS R-tree packed using the
 Sort-Tile-Recursive algorithm. Pass a list of geometry objects to the STRtree
-constructor to create an R-tree that you can query with another geometric object.
+constructor to create a spatial index that you can query with another geometric
+object. Query-only means that once created, the `STRtree` is immutable. You
+cannot add or remove geometries.
 
 .. class:: strtree.STRtree(geometries)
 
   The `STRtree` constructor takes a sequence of geometric objects.
 
-  These are copied and stored in the R-tree.
+  References to these geometric objects are kept and stored in the R-tree.
 
   `New in version 1.4.0`.
 
-The `query` method on `STRtree` returns a list of all geometries in the tree that
-intersect the provided geometry argument. If you want to match geometries of a
-more specific spatial relationship (eg. crosses, contains, overlaps), consider
-performing the query on the R-tree, followed by a manual search through the
-returned subset using the desired binary predicate.
+  .. method:: strtree.query(goem)
 
-Query-only means that once created, the R-tree is immutable. You cannot 
-add or remove geometries.
+    Returns a list of all geometries in the `strtree` whose extents intersect the
+    extent of `geom`. This means that a subsequent search through the returned
+    subset using the desired binary predicate (eg. intersects, crosses, contains,
+    overlaps) may be necessary to further filter the results according to their
+    specific spatial relationships.
 
-.. code-block:: pycon
+    .. code-block:: pycon
 
-  >>> from shapely.geometry import Point
-  >>> from shapely.strtree import STRtree
-  >>> points = [Point(i, i) for i in range(10)]
-  >>> tree = STRtree(points)
-  >>> tree.query(Point(2,2).buffer(0.99))
-  >>> [o.wkt for o in tree.query(Point(2,2).buffer(0.99))]
-  ['POINT (2 2)']
-  >>> [o.wkt for o in tree.query(Point(2,2).buffer(1.0))]
-  ['POINT (1 1)', 'POINT (2 2)', 'POINT (3 3)']
+      >>> from shapely.strtree import STRtree
+      >>> points = [Point(i, i) for i in range(10)]
+      >>> tree = STRtree(points)
+      >>> query_geom = Point(2,2).buffer(0.99)
+      >>> [o.wkt for o in tree.query(query_geom)]
+      ['POINT (2 2)']
+      >>> query_geom = Point(2, 2).buffer(1.0)
+      >>> [o.wkt for o in tree.query(query_geom)]
+      ['POINT (1 1)', 'POINT (2 2)', 'POINT (3 3)']
+      >>> [o.wkt for o in tree.query(query_geom) if o.intersects(query_geom)]
+      ['POINT (2 2)']
+
+    .. note::
+      To get the original indexes of the query results, create an auxiliary
+      dictionary. But use the geometry `ids` as keys since the shapely geometries
+      themselves are not hashable.
+
+      .. code-block:: pycon
+
+        >>> index_by_id = dict((id(pt), i) for i, pt in enumerate(points))
+        >>> [(index_by_id[id(pt)], pt.wkt) for pt in tree.query(Point(2,2).buffer(1.0))]
+        [(1, 'POINT (1 1)'), (2, 'POINT (2 2)'), (3, 'POINT (3 3)')]
+
+
+  .. method:: strtree.nearest(geom)
+
+    Returns the nearest geometry in `strtree` to `geom`.
+
+    .. code-block:: pycon
 
+      >>> tree = STRtree([Point(i, i) for i in range(10)])
+      >>> tree.nearest(Point(2.2, 2.2)).wkt
+      'Point (2 2)'
 
 Interoperation
 ==============
diff --git a/shapely/geometry/geo.py b/shapely/geometry/geo.py
index ce2a602..c78d355 100644
--- a/shapely/geometry/geo.py
+++ b/shapely/geometry/geo.py
@@ -1,7 +1,6 @@
 """
 Geometry factories based on the geo interface
 """
-
 from .point import Point, asPoint
 from .linestring import LineString, asLineString
 from .polygon import Polygon, asPolygon
@@ -10,6 +9,47 @@ from .multilinestring import MultiLineString, asMultiLineString
 from .multipolygon import MultiPolygon, MultiPolygonAdapter
 from .collection import GeometryCollection
 
+# numpy is an optional dependency
+try:
+    import numpy as np
+except ImportError:
+    _has_numpy = False
+else:
+    _has_numpy = True
+
+
+def _is_coordinates_empty(coordinates):
+    """Helper to identify if coordinates or subset of coordinates are empty"""
+
+    if coordinates is None:
+        return True
+
+    is_numpy_array = _has_numpy and isinstance(coordinates, np.ndarray)
+    if isinstance(coordinates, (list, tuple)) or is_numpy_array:
+        if len(coordinates) == 0:
+            return True
+        return all(map(_is_coordinates_empty, coordinates))
+    else:
+        return False
+
+
+def _empty_shape_for_no_coordinates(geom_type):
+    """Return empty counterpart for geom_type"""
+    if geom_type == 'point':
+        return Point()
+    elif geom_type == 'multipoint':
+        return MultiPoint()
+    elif geom_type == 'linestring':
+        return LineString()
+    elif geom_type == 'multilinestring':
+        return MultiLineString()
+    elif geom_type == 'polygon':
+        return Polygon()
+    elif geom_type == 'multipolygon':
+        return MultiPolygon()
+    else:
+        raise ValueError("Unknown geometry type: %s" % geom_type)
+
 
 def box(minx, miny, maxx, maxy, ccw=True):
     """Returns a rectangular polygon with configurable normal vector"""
@@ -18,6 +58,7 @@ def box(minx, miny, maxx, maxy, ccw=True):
         coords = coords[::-1]
     return Polygon(coords)
 
+
 def shape(context):
     """Returns a new, independent geometry with coordinates *copied* from the
     context.
@@ -27,15 +68,14 @@ def shape(context):
     else:
         ob = context
     geom_type = ob.get("type").lower()
-    if geom_type == "point":
+    if 'coordinates' in ob and _is_coordinates_empty(ob['coordinates']):
+        return _empty_shape_for_no_coordinates(geom_type)
+    elif geom_type == "point":
         return Point(ob["coordinates"])
     elif geom_type == "linestring":
         return LineString(ob["coordinates"])
     elif geom_type == "polygon":
-        if not ob["coordinates"]:
-            return Polygon()
-        else:
-            return Polygon(ob["coordinates"][0], ob["coordinates"][1:])
+        return Polygon(ob["coordinates"][0], ob["coordinates"][1:])
     elif geom_type == "multipoint":
         return MultiPoint(ob["coordinates"])
     elif geom_type == "multilinestring":
@@ -48,6 +88,7 @@ def shape(context):
     else:
         raise ValueError("Unknown geometry type: %s" % geom_type)
 
+
 def asShape(context):
     """Adapts the context to a geometry interface. The coordinates remain
     stored in the context.
@@ -80,6 +121,7 @@ def asShape(context):
     else:
         raise ValueError("Unknown geometry type: %s" % geom_type)
 
+
 def mapping(ob):
     """Returns a GeoJSON-like mapping"""
     return ob.__geo_interface__


Test patch:
diff --git a/tests/test_collection.py b/tests/test_collection.py
index 1e6f0a6..c87fd7a 100644
--- a/tests/test_collection.py
+++ b/tests/test_collection.py
@@ -63,6 +63,19 @@ class CollectionTestCase(unittest.TestCase):
         self.assertEqual(len(m), 0)
         self.assertEqual(m.geoms, [])
 
+    def test_empty_coordinates(self):
+
+        d = {"type": "GeometryCollection", "geometries": [
+            {"type": "Point", "coordinates": ()},
+            {"type": "LineString", "coordinates": (())}
+        ]}
+
+        # shape
+        m = shape(d)
+        self.assertEqual(m.geom_type, "GeometryCollection")
+        self.assertEqual(len(m), 0)
+        self.assertEqual(m.geoms, [])
+
 
 def test_suite():
     return unittest.TestLoader().loadTestsFromTestCase(CollectionTestCase)
diff --git a/tests/test_shape.py b/tests/test_shape.py
new file mode 100644
index 0000000..cf35395
--- /dev/null
+++ b/tests/test_shape.py
@@ -0,0 +1,47 @@
+from . import unittest
+from shapely.geometry import shape, Polygon
+
+# numpy is an optional dependency
+try:
+    import numpy as np
+except ImportError:
+    _has_numpy = False
+else:
+    _has_numpy = True
+
+
+class ShapeTestCase(unittest.TestCase):
+
+    def test_polygon_no_coords(self):
+        # using None
+        d = {"type": "Polygon", "coordinates": None}
+        p = shape(d)
+        self.assertEqual(p, Polygon())
+
+        # using empty list
+        d = {"type": "Polygon", "coordinates": []}
+        p = shape(d)
+        self.assertEqual(p, Polygon())
+
+        # using empty array
+        d = {"type": "Polygon", "coordinates": np.array([])}
+        p = shape(d)
+        self.assertEqual(p, Polygon())
+
+    def test_polygon_with_coords_list(self):
+        # list
+        d = {"type": "Polygon", "coordinates": [[[5, 10], [10, 10], [10, 5]]]}
+        p = shape(d)
+        self.assertEqual(p, Polygon([(5, 10), (10, 10), (10, 5)]))
+
+        # numpy array
+        if _has_numpy:
+            d = {"type": "Polygon", "coordinates": np.array([[[5, 10],
+                                                              [10, 10],
+                                                              [10, 5]]])}
+            p = shape(d)
+            self.assertEqual(p, Polygon([(5, 10), (10, 10), (10, 5)]))
+
+
+def test_suite():
+    return unittest.TestLoader().loadTestsFromTestCase(ShapeTestCase)