Repository: Toblerity/Shapely (python)

Issue:
Expose GEOSNormalize
`GEOSNormalize_r` is the C API for [`Geometry::normalize`](https://geos.osgeo.org/doxygen/classgeos_1_1geom_1_1Geometry.html#a2ea3b58d562bc183e118e2a68f86414d), which reorganizes geometries into 'canonical forms': points in a `LineString` may be deduplicated, members of a collection are sorted by type, points in a polygon's rings or a single `LinearRing` start at the top-left-most point and go clockwise, etc. This can be useful for unit testing or geometry deduplication. This was added [in libgeos 3.0.0](https://github.com/libgeos/geos/commit/9ecb499d71c279f4713f4aa8a2c0e9d4656b0e2c).

I tried to add the method myself but I never really played with `ctypes`, so this was mostly playing around, and I always got `-1` (exception). This method updates the geometry in place and I am not sure if this is a good thing; pygeos [implemented](https://github.com/pygeos/pygeos/blob/7a2798ab8ef43acf59e2023b8cfcb49e6d9384d3/src/ufuncs.c#L257) a wrapper that clones before normalizing.

Note that the `LineString::normalize` implementation has a bug in which the first point of a `LinearRing` *was not* the top-left-most point. It was fixed in [this commit](https://github.com/libgeos/geos/commit/a807833875976a29fbcf509297f4a7a5f837aa09#diff-a6950e44d48810137ee9f1460388168dR323), which is not included in any release yet (probably will be in libgeos 3.9.0). The workaround is to create a polygon using the LinearRing as an exterior ring, normalize the polygon, and get the ring back. Assuming normalize would be added to Shapely, we would get something like this:

```python3
>>> list(Polygon(((8, 8), (8, 5), (8, 5), (5, 5), (5, 8))).normalize().exterior.coords)
[(8, 8), (8, 5), (5, 5), (5, 8), (8, 8)]
>>> list(LinearRing(((8, 8), (8, 5), (8, 5), (5, 5), (5, 8))).normalize().coords)
[(5, 5), (5, 8), (8, 8), (8, 5), (5, 5)]
```

Pull request description:
ENH: Add normalize() method
Closes #994

This exposes `GEOSNormalize` as a method on the geometry classes. In GEOS, this works in place, but to keep the Shapely interface consistent (other operations like simplify also return a new geometry), I made the `normalize()` method to return a new geometry (and thus first clone the geometry before normalizing it).

This seems to work for me for some simple cases. cc @tomplex @Lucidiot  

Especially given the overlay changes in GEOS 3.9.0 (where you often first need to normalize the geometry if you want to be able to compare it to an expected result across different GEOS versions), I think it would be useful to still include this in Shapely 1.8.

Patch:
diff --git a/docs/manual.rst b/docs/manual.rst
index dd6e8fd..6d81e86 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -1182,6 +1182,20 @@ A line's endpoints are part of its `boundary` and are therefore not contained.
   >>> [p.wkt for p in contained]
   ['POINT (0.5000000000000000 0.5000000000000000)']
 
+.. method:: object.covers(other)
+
+  Returns ``True`` if every point of `other` is a point on the interior or
+  boundary of `object`. This is similar to ``object.contains(other)`` except
+  that this does not require any interior points of `other` to lie in the 
+  interior of `object`.
+
+.. method:: object.covered_by(other)
+
+  Returns ``True`` if every point of `object` is a point on the interior or
+  boundary of `other`. This is equivalent to ``other.covers(object)``.
+  
+  `New in version 1.8`.
+
 .. method:: object.crosses(other)
 
   Returns ``True`` if the `interior` of the object intersects the `interior` of
diff --git a/shapely/ctypes_declarations.py b/shapely/ctypes_declarations.py
index 51ee01e..a0a4f00 100644
--- a/shapely/ctypes_declarations.py
+++ b/shapely/ctypes_declarations.py
@@ -388,6 +388,9 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSGeom_getDimensions.restype = c_int
     lgeos.GEOSGeom_getDimensions.argtypes = [c_void_p]
 
+    lgeos.GEOSNormalize.restype = c_int
+    lgeos.GEOSNormalize.argtypes = [c_void_p]
+
     # Misc functions
 
     lgeos.GEOSArea.restype = c_double
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index fc4258f..b97de9a 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -637,6 +637,27 @@ class BaseGeometry(object):
             op = self.impl['simplify']
         return geom_factory(op(self, tolerance))
 
+    def normalize(self):
+        """Converts geometry to normal form (or canonical form).
+
+        This method orders the coordinates, rings of a polygon and parts of
+        multi geometries consistently. Typically useful for testing purposes
+        (for example in combination with `equals_exact`).
+
+        Examples
+        --------
+        >>> from shapely.wkt import loads
+        >>> p = loads("MULTILINESTRING((0 0, 1 1), (3 3, 2 2))")
+        >>> p.normalize().wkt
+        'MULTILINESTRING ((2 2, 3 3), (0 0, 1 1))'
+        """
+        # self.impl['normalize'](self)
+        if self._geom is None:
+            raise ValueError("Null geometry supports no operations")
+        geom_cloned = lgeos.GEOSGeom_clone(self._geom)
+        lgeos.GEOSNormalize(geom_cloned)
+        return geom_factory(geom_cloned)
+
     # Binary operations
     # -----------------
 
diff --git a/shapely/geos.py b/shapely/geos.py
index cb8b232..d7e766f 100644
--- a/shapely/geos.py
+++ b/shapely/geos.py
@@ -699,6 +699,7 @@ class LGEOS330(LGEOSBase):
         self.methods['simplify'] = self.GEOSSimplify
         self.methods['topology_preserve_simplify'] = \
             self.GEOSTopologyPreserveSimplify
+        self.methods['normalize'] = self.GEOSNormalize
         self.methods['cascaded_union'] = self.GEOSUnionCascaded
 
         def parallel_offset(geom, distance, resolution=16, join_style=1,
diff --git a/shapely/impl.py b/shapely/impl.py
index d35a6da..9b206e9 100644
--- a/shapely/impl.py
+++ b/shapely/impl.py
@@ -88,6 +88,7 @@ IMPL330 = {
     'simplify': (UnaryTopologicalOp, 'simplify'),
     'topology_preserve_simplify':
         (UnaryTopologicalOp, 'topology_preserve_simplify'),
+    'normalize': (UnaryTopologicalOp, 'normalize'),
     #
     'difference': (BinaryTopologicalOp, 'difference'),
     'intersection': (BinaryTopologicalOp, 'intersection'),
diff --git a/shapely/strtree.py b/shapely/strtree.py
index c528837..70bee19 100644
--- a/shapely/strtree.py
+++ b/shapely/strtree.py
@@ -82,8 +82,8 @@ class STRtree:
         self._geoms = list(geoms)
 
     def _init_tree_handle(self, geoms):
-        # GEOS STRtree capacity has to be > 1
-        self._tree_handle = lgeos.GEOSSTRtree_create(max(2, len(geoms)))
+        node_capacity = 10
+        self._tree_handle = lgeos.GEOSSTRtree_create(node_capacity)
         for geom in geoms:
             lgeos.GEOSSTRtree_insert(self._tree_handle, geom._geom, ctypes.py_object(geom))
 


Test patch:
diff --git a/tests/test_operations.py b/tests/test_operations.py
index b1d5b88..23dfd08 100644
--- a/tests/test_operations.py
+++ b/tests/test_operations.py
@@ -92,3 +92,13 @@ class OperationsTestCase(unittest.TestCase):
         assert(empty_line.is_empty)
         with pytest.raises(ValueError):
             empty_line.interpolate(.5, normalized=True)
+
+    def test_normalize(self):
+        point = Point(1, 1)
+        result = point.normalize()
+        assert result == point
+
+        line = loads("MULTILINESTRING ((1 1, 0 0), (1 1, 1 2))")
+        result = line.normalize()
+        expected = loads("MULTILINESTRING ((1 1, 1 2), (0 0, 1 1))")
+        assert result == expected