Repository: Toblerity/Shapely
Pull request #1034 (2020-11-19 01:53:22)

ops.split: Split the line when the boundary of a splitter LineString touches the interior of the line
See #701 for background. 

I thought about using `splitter.relate_pattern(line, '...')` several times. But I'm not sure how much if any of that operation is cached, so I'm calculating one relation string and using it directly for the three checks. Fortunately the checks are all straightforward, each involving only one character.

Ping to @snorfalorpagus, just in case you would want to take a look since you have been heavily involved with the split function in the past.

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

Coverage increased (+0.005%) to 83.525% when pulling **494a8bf7678889f6041c3fe1c172cdfde11ce510 on Jeremiah-England:llsplit** into **6aeb8e7932076fc8700a3c36133255957989464c on Toblerity:master**.


Patch:
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index d0ca286..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,125 +0,0 @@
-language: python
-dist: xenial
-
-cache:
-  directories:
-    - $HOME/geosinstall
-    - ~/.cache/pip
-
-matrix:
-  include:
-    - python: "3.5"
-      env:
-        GEOS_VERSION="3.5.2"
-        SPEEDUPS=1
-        NUMPY=1
-    - python: "3.5"
-      env:
-        GEOS_VERSION="3.5.2"
-        SPEEDUPS=0
-        NUMPY=1
-    - python: "3.5"
-      env:
-        GEOS_VERSION="3.5.2"
-        SPEEDUPS=0
-        NUMPY=0
-    - python: "3.6"
-      env:
-        GEOS_VERSION="3.6.4"
-        SPEEDUPS=1
-        NUMPY=1
-    - python: "3.6"
-      env:
-        GEOS_VERSION="3.6.4"
-        SPEEDUPS=0
-        NUMPY=1
-    - python: "3.6"
-      env:
-        GEOS_VERSION="3.6.4"
-        SPEEDUPS=0
-        NUMPY=0
-    - python: "3.7"
-      env:
-        GEOS_VERSION="3.7.3"
-        SPEEDUPS=1
-        NUMPY=1
-    - python: "3.7"
-      env:
-        GEOS_VERSION="3.7.3"
-        SPEEDUPS=0
-        NUMPY=1
-    - python: "3.7"
-      env:
-        GEOS_VERSION="3.7.3"
-        SPEEDUPS=0
-        NUMPY=0
-    - python: "3.8"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=1
-        NUMPY=1
-    - python: "3.8"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=0
-        NUMPY=1
-    - python: "3.8"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=0
-        NUMPY=0
-    - python: "3.9"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=1
-        NUMPY=1
-    - python: "3.9"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=0
-        NUMPY=1
-    - python: "3.9"
-      env:
-        GEOS_VERSION="3.8.1"
-        SPEEDUPS=0
-        NUMPY=0
-    - python: "3.9-dev"
-      env:
-        GEOS_VERSION="master"
-        SPEEDUPS=1
-        NUMPY=1
-  allow_failures:
-    - python: "3.9-dev"
-      env:
-        GEOS_VERSION="master"
-        SPEEDUPS=1
-        NUMPY=1
-
-before_install:
-  - ./ci/travis/install_geos.sh
-  - pip install --disable-pip-version-check --upgrade pip
-  - pip install --upgrade wheel
-  # if building with speedups install cython
-  - if [ "$SPEEDUPS" == "1" ]; then pip install --install-option="--no-cython-compile" cython; fi
-  # if testing without numpy explicitly remove it
-  - if [ "$NUMPY" == "0" ]; then pip uninstall --yes numpy; fi
-  # convert SPEEDUPS to --with-speedups/--without-speedups
-  - if [ "$SPEEDUPS" == "1" ]; then SPEEDUPS_FLAG=--with-speedups; else SPEEDUPS_FLAG=--without-speedups; fi
-  - pip install --upgrade coveralls pytest-cov pytest>=3.8
-
-install:
-  - export GEOS_CONFIG=$HOME/geosinstall/geos-$GEOS_VERSION/bin/geos-config
-  - pip install -v -e .[all]
-
-script:
-  - export LD_LIBRARY_PATH=$HOME/geosinstall/geos-$GEOS_VERSION/lib
-  - export DYLD_LIBRARY_PATH=$HOME/geosinstall/geos-$GEOS_VERSION/lib
-  - python -c "from shapely.geos import geos_version; print(geos_version)"
-  - python -m pytest --cov shapely --cov-report term-missing "${SPEEDUPS_FLAG}"
-  - if [[ "$TRAVIS_PYTHON_VERSION" == "3.8" && "$NUMPY" == "1" && "$SPEEDUPS" == "1" ]]; then python -m pytest shapely --doctest-modules; fi
-
-after_success:
-  - coveralls || echo "!! intermittent coveralls failure"
-
-notifications:
-    email: false
diff --git a/CHANGES.txt b/CHANGES.txt
index 03877d0..5339dce 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,15 +1,30 @@
 Changes
 =======
 
-1.8.0 (TBD)
+1.8a1 (TBD)
 -----------
 
+Deprecations:
+
+- ops.cascaded_union will be removed in version 2.0.0.
+- More ...
+
 New features:
 
+- Add a normalize() method to geometry classes, exposing the GEOSNormalize
+  algorithm (#1090).
+- Initialize STRtree with a capacity of 10 items per node (#1070).
+- Load libraries relocated to shapely/.libs by auditwheel versions < 3.1 or
+  relocated to Shapely.libs by auditwheel versions >= 3.1.
 - shapely.ops.voronoi_diagram() computes the Voronoi Diagram of a geometry or
   geometry collection (#833, #851).
 - shapely.validation.make_valid() fixes invalid geometries (#883)
 
+Bug fixes:
+
+- Several issues with duplication and order of vertices in ops.substring have
+  been fixed (#1008).
+
 1.7.1 (2020-08-20)
 ------------------
 
diff --git a/CITATION.txt b/CITATION.txt
index b3ea848..d3c3946 100644
--- a/CITATION.txt
+++ b/CITATION.txt
@@ -1,7 +1,7 @@
 If you use Shapely for any published work, please cite it using the reference
 below:
 
-@Misc{,
+@Misc{shapely,
   author =    {Sean Gillies and others},
   organization = {toblerity.org},
   title =     {Shapely: manipulation and analysis of geometric objects},
diff --git a/README.rst b/README.rst
index 4adbb71..3005d66 100644
--- a/README.rst
+++ b/README.rst
@@ -2,10 +2,10 @@
 Shapely
 =======
 
-|travis| |appveyor| |coveralls|
+|github-actions| |appveyor| |coveralls|
 
-.. |travis| image:: https://travis-ci.org/Toblerity/Shapely.svg?branch=master
-   :target: https://travis-ci.org/Toblerity/Shapely
+.. |github-actions| image:: https://github.com/Toblerity/Shapely/workflows/Tests%20(Linux)/badge.svg
+   :target: https://github.com/Toblerity/Shapely/actions
 
 .. |appveyor| image:: https://ci.appveyor.com/api/projects/status/github/Toblerity/Shapely?branch=master&svg=true
    :target: https://ci.appveyor.com/project/frsci/shapely?branch=master
diff --git a/appveyor.yml b/appveyor.yml
index 0a09480..77365bf 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -1,8 +1,11 @@
+---
 image: Visual Studio 2019
 platform: x64
 
 environment:
-  COMPILER: msvc2019
+  global:
+    COMPILER: msvc2019
+
   matrix:
     # Interleave PYTHON, ARCH and GEOS_VERSION
     - PYTHON: "C:\\Python35-x64"
@@ -19,7 +22,7 @@ environment:
       GEOS_VERSION: "3.8.1"
     - PYTHON: "C:\\Python39-x64"
       ARCH: x64
-      GEOS_VERSION: "3.8.1"
+      GEOS_VERSION: "3.9.1"
 
 
 install:
@@ -33,7 +36,7 @@ install:
         Where-Object pullRequestId -eq $env:APPVEYOR_PULL_REQUEST_NUMBER)[0].buildNumber) { `
           throw "There are newer queued builds for this pull request, failing early." }
 
-
+build: false  # disable automatic build
 build_script:
   - set GEOS_INSTALL=C:\projects\deps\geos-%GEOS_VERSION%-%COMPILER%-%ARCH%
   - set CYTHONPF=C:\\projects\\deps\\geos-%GEOS_VERSION%-%COMPILER%-%ARCH%
@@ -49,7 +52,7 @@ build_script:
 
   - ps: 'Write-Host "Checking GEOS build $env:GEOS_INSTALL" -ForegroundColor Magenta'
   - if not exist C:\projects\deps mkdir C:\projects\deps
-  - call %APPVEYOR_BUILD_FOLDER%\ci\appveyor\build_geos.cmd
+  - call %APPVEYOR_BUILD_FOLDER%\ci\install_geos.cmd
 
   - ps: 'Write-Host "Checking geos_c.dll" -ForegroundColor Magenta'
   - python -c "import sys; print('Python ' + sys.version)"
@@ -68,4 +71,4 @@ build_script:
 
 
 cache:
-  - C:\projects\deps -> %APPVEYOR_BUILD_FOLDER%\ci\appveyor\build_geos.cmd
+  - C:\projects\deps -> %APPVEYOR_BUILD_FOLDER%\ci\install_geos.cmd
diff --git a/ci/appveyor/build_geos.cmd b/ci/appveyor/build_geos.cmd
deleted file mode 100644
index 4dbaf27..0000000
--- a/ci/appveyor/build_geos.cmd
+++ /dev/null
@@ -1,25 +0,0 @@
-REM This script is called from appveyor.yml
-
-if exist %GEOS_INSTALL% (
-  echo Using cached %GEOS_INSTALL%
-) else (
-  echo Building %GEOS_INSTALL%
-
-  cd C:\projects
-
-  curl -fsSO http://download.osgeo.org/geos/geos-%GEOS_VERSION%.tar.bz2
-  7z x geos-%GEOS_VERSION%.tar.bz2
-  7z x geos-%GEOS_VERSION%.tar
-  cd geos-%GEOS_VERSION% || exit /B 5
-
-  pip install ninja
-  cmake --version
-
-  mkdir build
-  cd build
-  cmake -GNinja -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=%GEOS_INSTALL% .. || exit /B 1
-  cmake --build . --config Release || exit /B 2
-  ctest . --config Release || exit /B 3
-  cmake --install . --config Release || exit /B 4
-  cd ..
-)
diff --git a/ci/install_geos.cmd b/ci/install_geos.cmd
new file mode 100644
index 0000000..24f78dd
--- /dev/null
+++ b/ci/install_geos.cmd
@@ -0,0 +1,28 @@
+:: Build and install GEOS on Windows system, save cache for later
+::
+:: This script requires environment variables to be set
+::  - set GEOS_INSTALL=C:\path\to\cached\prefix -- to build or use as cache
+::  - set GEOS_VERSION=3.7.3 -- to download and compile
+
+if exist %GEOS_INSTALL% (
+  echo Using cached %GEOS_INSTALL%
+) else (
+  echo Building %GEOS_INSTALL%
+
+  curl -fsSO http://download.osgeo.org/geos/geos-%GEOS_VERSION%.tar.bz2
+  7z x geos-%GEOS_VERSION%.tar.bz2
+  7z x geos-%GEOS_VERSION%.tar
+  cd geos-%GEOS_VERSION% || exit /B 1
+
+  pip install ninja
+  cmake --version
+
+  mkdir build
+  cd build
+  cmake -GNinja -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=%GEOS_INSTALL% .. || exit /B 2
+  cmake --build . || exit /B 3
+  ctest . || exit /B 4
+  ctest .
+  cmake --install . || exit /B 5
+  cd ..
+)
diff --git a/ci/install_geos.sh b/ci/install_geos.sh
new file mode 100755
index 0000000..b39e093
--- /dev/null
+++ b/ci/install_geos.sh
@@ -0,0 +1,70 @@
+#!/bin/sh
+
+# Build and install GEOS on a POSIX system, save cache for later
+#
+# This script requires environment variables to be set
+#  - export GEOS_INSTALL=/path/to/cached/prefix -- to build or use as cache
+#  - export GEOS_VERSION=3.7.3 or master -- to download and compile
+
+set -e
+
+if [ -z "$GEOS_INSTALL" ]; then
+    echo "GEOS_INSTALL must be set"
+    exit 1
+elif [ -z "$GEOS_VERSION" ]; then
+    echo "GEOS_VERSION must be set"
+    exit 1
+fi
+
+# Create directories, if they don't exit
+mkdir -p $GEOS_INSTALL
+
+# Download and build GEOS outside other source tree
+GEOS_BUILD=$HOME/geosbuild
+
+prepare_geos_build_dir(){
+  rm -rf $GEOS_BUILD
+  mkdir -p $GEOS_BUILD
+  cd $GEOS_BUILD
+}
+
+build_geos(){
+    echo "Building geos-$GEOS_VERSION"
+    mkdir build
+    cd build
+    cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$GEOS_INSTALL ..
+    make -j 2
+    ctest .
+    make install
+}
+
+if [ "$GEOS_VERSION" = "master" ]; then
+    prepare_geos_build_dir
+    # use GitHub mirror
+    git clone --depth 1 https://github.com/libgeos/geos.git geos-$GEOS_VERSION
+    cd geos-$GEOS_VERSION
+    git rev-parse HEAD > newrev.txt
+    BUILD=no
+    # Only build if nothing cached or if the GEOS revision changed
+    if test ! -f $GEOS_INSTALL/rev.txt; then
+        BUILD=yes
+    elif ! diff newrev.txt $GEOS_INSTALL/rev.txt >/dev/null; then
+        BUILD=yes
+    fi
+    if test "$BUILD" = "no"; then
+        echo "Using cached install $GEOS_INSTALL"
+    else
+        cp newrev.txt $GEOS_INSTALL/rev.txt
+        build_geos
+    fi
+else
+    if [ -d "$GEOS_INSTALL/include/geos" ]; then
+        echo "Using cached install $GEOS_INSTALL"
+    else
+        prepare_geos_build_dir
+        wget -q http://download.osgeo.org/geos/geos-$GEOS_VERSION.tar.bz2
+        tar xfj geos-$GEOS_VERSION.tar.bz2
+        cd geos-$GEOS_VERSION
+        build_geos
+    fi
+fi
diff --git a/ci/travis/install_geos.sh b/ci/travis/install_geos.sh
deleted file mode 100755
index db75afc..0000000
--- a/ci/travis/install_geos.sh
+++ /dev/null
@@ -1,49 +0,0 @@
-#!/bin/bash
-# based on contribution by @rbuffat to Toblerity/Fiona
-set -e
-
-# Cache install directory from .travis.yml
-CACHE_GEOS_INST=$HOME/geosinstall
-
-# Create directories, if they don't exit
-GEOS_INST_VERSION=$CACHE_GEOS_INST/geos-$GEOS_VERSION
-mkdir -p $GEOS_INST_VERSION
-
-function build_geos {
-    echo "Building geos-$GEOS_VERSION"
-    mkdir build
-    cd build
-    cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$GEOS_INST_VERSION ..
-    make -j 2
-    ctest .
-    make install
-}
-
-if [ "$GEOS_VERSION" = "master" ]; then
-    # use GitHub mirror
-    git clone --depth 1 https://github.com/libgeos/geos.git geos-$GEOS_VERSION
-    cd geos-$GEOS_VERSION
-    git rev-parse HEAD > newrev.txt
-    BUILD=no
-    # Only build if nothing cached or if the GEOS revision changed
-    if test ! -f $GEOS_INST_VERSION/rev.txt; then
-        BUILD=yes
-    elif ! diff newrev.txt $GEOS_INST_VERSION/rev.txt >/dev/null; then
-        BUILD=yes
-    fi
-    if test "$BUILD" = "no"; then
-        echo "Using cached install $GEOS_INST_VERSION"
-    else
-        cp newrev.txt $GEOS_INST_VERSION/rev.txt
-        build_geos
-    fi
-else
-    if [ -d "$GEOS_INST_VERSION/include/geos" ]; then
-        echo "Using cached install $GEOS_INST_VERSION"
-    else
-        wget -q http://download.osgeo.org/geos/geos-$GEOS_VERSION.tar.bz2
-        tar xfj geos-$GEOS_VERSION.tar.bz2
-        cd geos-$GEOS_VERSION
-        build_geos
-    fi
-fi
diff --git a/docs/manual.rst b/docs/manual.rst
index 187f1c4..6d81e86 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -775,14 +775,12 @@ Its `x-y` bounding box is a ``(minx, miny, maxx, maxy)`` tuple.
   (-1.0, -1.0, 2.0, 2.0)
 
 Its members are instances of `Polygon` and are accessed via the ``geoms``
-property or via the iterator protocol using ``in`` or ``list()``.
+property.
 
 .. code-block:: pycon
 
   >>> len(polygons.geoms)
   3
-  >>> len(polygons)
-  3
 
 .. _empties:
 
@@ -876,10 +874,6 @@ supports linear referencing based on length or distance, evaluating the
 distance along a geometric object to the projection of a given point, or the
 point at a given distance along the object.
 
-.. note::
-
-  Linear referencing methods require GEOS 3.2.0 or later.
-
 .. method:: object.interpolate(distance[, normalized=False])
 
   Return a point at the specified distance along a linear geometric object.
@@ -1188,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
@@ -1653,11 +1661,11 @@ mileage may vary when cleaning surfaces.
   True
   >>> clean
   <shapely.geometry.multipolygon.MultiPolygon object at ...>
-  >>> len(clean)
+  >>> len(clean.geoms)
   2
-  >>> list(clean[0].exterior.coords)
+  >>> list(clean.geoms[0].exterior.coords)
   [(0.0, 0.0), (0.0, 2.0), (1.0, 1.0), (0.0, 0.0)]
-  >>> list(clean[1].exterior.coords)
+  >>> list(clean.geoms[1].exterior.coords)
   [(1.0, 1.0), (2.0, 2.0), (2.0, 0.0), (1.0, 1.0)]
 
 Buffering splits the polygon in two at the point where they touch.
@@ -2151,7 +2159,7 @@ using functions in the :mod:`shapely.ops` module.
     ...     ((1, 1), (100, 100)),
     ...     ]
     >>> result, dangles, cuts, invalids = polygonize_full(lines)
-    >>> len(result)
+    >>> len(result.geoms)
     2
     >>> list(result.geoms)
     [<shapely.geometry.polygon.Polygon object at ...>, <shapely.geometry.polygon.Polygon object at ...>]
@@ -2356,7 +2364,7 @@ Shared paths
 ------------
 
 The :func:`~shapely.ops.shared_paths` function in `shapely.ops` finds the shared
-paths between two lineal geometries.
+paths between two linear geometries.
 
 .. function:: shapely.ops.shared_paths(geom1, geom2)
 
@@ -2414,27 +2422,53 @@ The :func:`~shapely.ops.split` function in `shapely.ops` splits a geometry by an
 Substring
 ---------
 
-The :func:`~shapely.ops.substring` function in `shapely.ops` returns a line segment between specified distances along a linear geometry.
+The :func:`~shapely.ops.substring` function in :mod:`shapely.ops` returns a line segment
+between specified distances along a `LineString`.
 
 .. function:: shapely.ops.substring(geom, start_dist, end_dist[, normalized=False])
 
-    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.
+  Return the `LineString` between `start_dist` and `end_dist` or a `Point`
+  if they are at the same location
 
-    If the start distance equals the end distance, a point is being returned.
+  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 normalized arg is True, the distance will be interpreted as a
-    fraction of the geometry's length.
+  If the start distance equals the end distance, a point is being returned.
 
-   `New in version 1.7.0`
+  If the start distance is actually past the end distance, then the
+  reversed substring is returned such that the start distance is
+  at the first coordinate.
 
-.. code-block:: pycon
+  If the normalized arg is ``True``, the distance will be interpreted as a
+  fraction of the geometry's length
 
-  >>> line = LineString(([0, 0], [2, 0]))
-  >>> result = substring(line, 0.5, 0.6)
-  >>> result.wkt
-  'LINESTRING (0.5 0, 0.6 0)'
+  `New in version 1.7.0`
+
+  Here are some examples that return `LineString` geometries.
+
+  .. code-block:: pycon
+
+    >>> from shapely.geometry import LineString
+    >>> from shapely.ops import substring
+    >>> ls = LineString((i, 0) for i in range(6))
+    >>> ls.wkt
+    'LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)'
+    >>> substring(ls, start_dist=1, end_dist=3).wkt
+    'LINESTRING (1 0, 2 0, 3 0)'
+    >>> substring(ls, start_dist=3, end_dist=1).wkt
+    'LINESTRING (3 0, 2 0, 1 0)'
+    >>> substring(ls, start_dist=1, end_dist=-3).wkt
+    'LINESTRING (1 0, 2 0)'
+    >>> substring(ls, start_dist=0.2, end_dist=-0.6, normalized=True).wkt
+    'LINESTRING (1 0, 2 0)'
+
+  And here is an example that returns a `Point`.
+
+  .. code-block:: pycon
+
+    >>> substring(ls, start_dist=2.5, end_dist=-2.5)
+    'POINT (2.5 0)'
 
 Prepared Geometry Operations
 ----------------------------
diff --git a/shapely/ctypes_declarations.py b/shapely/ctypes_declarations.py
index adcd90d..a0a4f00 100644
--- a/shapely/ctypes_declarations.py
+++ b/shapely/ctypes_declarations.py
@@ -104,65 +104,51 @@ def prototype(lgeos, geos_version):
 
     # Linear referencing
 
-    if geos_version >= (3, 2, 0):
+    lgeos.GEOSProject.restype = c_double
+    lgeos.GEOSProject.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSProject.restype = c_double
-        lgeos.GEOSProject.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSInterpolate.restype = c_void_p
+    lgeos.GEOSInterpolate.argtypes = [c_void_p, c_double]
 
-        lgeos.GEOSInterpolate.restype = c_void_p
-        lgeos.GEOSInterpolate.argtypes = [c_void_p, c_double]
+    lgeos.GEOSProjectNormalized.restype = c_double
+    lgeos.GEOSProjectNormalized.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSProjectNormalized.restype = c_double
-        lgeos.GEOSProjectNormalized.argtypes = [c_void_p, c_void_p]
-
-        lgeos.GEOSInterpolateNormalized.restype = c_void_p
-        lgeos.GEOSInterpolateNormalized.argtypes = [c_void_p, c_double]
+    lgeos.GEOSInterpolateNormalized.restype = c_void_p
+    lgeos.GEOSInterpolateNormalized.argtypes = [c_void_p, c_double]
 
     # Buffer related
 
     lgeos.GEOSBuffer.restype = c_void_p
     lgeos.GEOSBuffer.argtypes = [c_void_p, c_double, c_int]
 
-    if geos_version >= (3, 2, 0):
-
-        lgeos.GEOSBufferWithStyle.restype = c_void_p
-        lgeos.GEOSBufferWithStyle.argtypes = [
-            c_void_p, c_double, c_int, c_int, c_int, c_double]
-
-        if geos_version >= (3, 3, 0):
-
-            lgeos.GEOSOffsetCurve.restype = c_void_p
-            lgeos.GEOSOffsetCurve.argtypes = [
-                c_void_p, c_double, c_int, c_int, c_double]
+    lgeos.GEOSBufferWithStyle.restype = c_void_p
+    lgeos.GEOSBufferWithStyle.argtypes = [
+        c_void_p, c_double, c_int, c_int, c_int, c_double]
 
-            lgeos.GEOSBufferParams_create.restype = c_void_p
-            lgeos.GEOSBufferParams_create.argtypes = None
+    lgeos.GEOSOffsetCurve.restype = c_void_p
+    lgeos.GEOSOffsetCurve.argtypes = [
+        c_void_p, c_double, c_int, c_int, c_double]
 
-            lgeos.GEOSBufferParams_setEndCapStyle.restype = c_int
-            lgeos.GEOSBufferParams_setEndCapStyle.argtypes = [c_void_p, c_int]
+    lgeos.GEOSBufferParams_create.restype = c_void_p
+    lgeos.GEOSBufferParams_create.argtypes = None
 
-            lgeos.GEOSBufferParams_setJoinStyle.restype = c_int
-            lgeos.GEOSBufferParams_setJoinStyle.argtypes = [c_void_p, c_int]
+    lgeos.GEOSBufferParams_setEndCapStyle.restype = c_int
+    lgeos.GEOSBufferParams_setEndCapStyle.argtypes = [c_void_p, c_int]
 
-            lgeos.GEOSBufferParams_setMitreLimit.restype = c_int
-            lgeos.GEOSBufferParams_setMitreLimit.argtypes = [c_void_p, c_double]
+    lgeos.GEOSBufferParams_setJoinStyle.restype = c_int
+    lgeos.GEOSBufferParams_setJoinStyle.argtypes = [c_void_p, c_int]
 
-            lgeos.GEOSBufferParams_setQuadrantSegments.restype = c_int
-            lgeos.GEOSBufferParams_setQuadrantSegments.argtypes = [c_void_p, c_int]
+    lgeos.GEOSBufferParams_setMitreLimit.restype = c_int
+    lgeos.GEOSBufferParams_setMitreLimit.argtypes = [c_void_p, c_double]
 
-            lgeos.GEOSBufferParams_setSingleSided.restype = c_int
-            lgeos.GEOSBufferParams_setSingleSided.argtypes = [c_void_p, c_int]
+    lgeos.GEOSBufferParams_setQuadrantSegments.restype = c_int
+    lgeos.GEOSBufferParams_setQuadrantSegments.argtypes = [c_void_p, c_int]
 
-            lgeos.GEOSBufferWithParams.restype = c_void_p
-            lgeos.GEOSBufferWithParams.argtypes = [
-                c_void_p, c_void_p, c_double]
+    lgeos.GEOSBufferParams_setSingleSided.restype = c_int
+    lgeos.GEOSBufferParams_setSingleSided.argtypes = [c_void_p, c_int]
 
-        else:
-
-            # deprecated in GEOS 3.3.0 in favour of GEOSOffsetCurve
-            lgeos.GEOSSingleSidedBuffer.restype = c_void_p
-            lgeos.GEOSSingleSidedBuffer.argtypes = [
-                c_void_p, c_double, c_int, c_int, c_double, c_int]
+    lgeos.GEOSBufferWithParams.restype = c_void_p
+    lgeos.GEOSBufferWithParams.argtypes = [c_void_p, c_void_p, c_double]
 
     # Geometry constructors
 
@@ -181,9 +167,8 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSGeom_createCollection.restype = c_void_p
     lgeos.GEOSGeom_createCollection.argtypes = [c_int, c_void_p, c_uint]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSGeom_createEmptyCollection.restype = c_void_p
-        lgeos.GEOSGeom_createEmptyCollection.argtypes = [c_int]
+    lgeos.GEOSGeom_createEmptyCollection.restype = c_void_p
+    lgeos.GEOSGeom_createEmptyCollection.argtypes = [c_int]
 
     lgeos.GEOSGeom_clone.restype = c_void_p
     lgeos.GEOSGeom_clone.argtypes = [c_void_p]
@@ -217,14 +202,12 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSUnion.restype = c_void_p
     lgeos.GEOSUnion.argtypes = [c_void_p, c_void_p]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSUnaryUnion.restype = c_void_p
-        lgeos.GEOSUnaryUnion.argtypes = [c_void_p]
+    lgeos.GEOSUnaryUnion.restype = c_void_p
+    lgeos.GEOSUnaryUnion.argtypes = [c_void_p]
 
-    if geos_version >= (3, 1, 0):
-        # deprecated in 3.3.0: use GEOSUnaryUnion instead
-        lgeos.GEOSUnionCascaded.restype = c_void_p
-        lgeos.GEOSUnionCascaded.argtypes = [c_void_p]
+    # deprecated in 3.3.0: use GEOSUnaryUnion instead
+    lgeos.GEOSUnionCascaded.restype = c_void_p
+    lgeos.GEOSUnionCascaded.argtypes = [c_void_p]
 
     lgeos.GEOSPointOnSurface.restype = c_void_p
     lgeos.GEOSPointOnSurface.argtypes = [c_void_p]
@@ -239,10 +222,9 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSPolygonize.restype = c_void_p
     lgeos.GEOSPolygonize.argtypes = [c_void_p, c_uint]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSPolygonize_full.restype = c_void_p
-        lgeos.GEOSPolygonize_full.argtypes = [
-            c_void_p, c_void_p, c_void_p, c_void_p]
+    lgeos.GEOSPolygonize_full.restype = c_void_p
+    lgeos.GEOSPolygonize_full.argtypes = [
+        c_void_p, c_void_p, c_void_p, c_void_p]
 
     if geos_version >= (3, 4, 0):
         lgeos.GEOSDelaunayTriangulation.restype = c_void_p
@@ -306,9 +288,8 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSisValid.restype = c_byte
     lgeos.GEOSisValid.argtypes = [c_void_p]
 
-    if geos_version >= (3, 1, 0):
-        lgeos.GEOSisValidReason.restype = allocated_c_char_p
-        lgeos.GEOSisValidReason.argtypes = [c_void_p]
+    lgeos.GEOSisValidReason.restype = allocated_c_char_p
+    lgeos.GEOSisValidReason.argtypes = [c_void_p]
 
     lgeos.GEOSisSimple.restype = c_byte
     lgeos.GEOSisSimple.argtypes = [c_void_p]
@@ -316,9 +297,8 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSisRing.restype = c_byte
     lgeos.GEOSisRing.argtypes = [c_void_p]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSisClosed.restype = c_byte
-        lgeos.GEOSisClosed.argtypes = [c_void_p]
+    lgeos.GEOSisClosed.restype = c_byte
+    lgeos.GEOSisClosed.argtypes = [c_void_p]
 
     lgeos.GEOSHasZ.restype = c_byte
     lgeos.GEOSHasZ.argtypes = [c_void_p]
@@ -331,47 +311,44 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSRelatePattern.restype = c_byte
     lgeos.GEOSRelatePattern.argtypes = [c_void_p, c_void_p, c_char_p]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSRelatePatternMatch.restype = c_byte
-        lgeos.GEOSRelatePatternMatch.argtypes = [c_char_p, c_char_p]
+    lgeos.GEOSRelatePatternMatch.restype = c_byte
+    lgeos.GEOSRelatePatternMatch.argtypes = [c_char_p, c_char_p]
 
     # Prepared Geometry Binary predicates
     # Return 2 on exception, 1 on true, 0 on false
 
-    if geos_version >= (3, 1, 0):
-
-        lgeos.GEOSPrepare.restype = c_void_p
-        lgeos.GEOSPrepare.argtypes = [c_void_p]
+    lgeos.GEOSPrepare.restype = c_void_p
+    lgeos.GEOSPrepare.argtypes = [c_void_p]
 
-        lgeos.GEOSPreparedGeom_destroy.restype = None
-        lgeos.GEOSPreparedGeom_destroy.argtypes = [c_void_p]
+    lgeos.GEOSPreparedGeom_destroy.restype = None
+    lgeos.GEOSPreparedGeom_destroy.argtypes = [c_void_p]
 
-        lgeos.GEOSPreparedDisjoint.restype = c_byte
-        lgeos.GEOSPreparedDisjoint.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedDisjoint.restype = c_byte
+    lgeos.GEOSPreparedDisjoint.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedTouches.restype = c_byte
-        lgeos.GEOSPreparedTouches.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedTouches.restype = c_byte
+    lgeos.GEOSPreparedTouches.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedIntersects.restype = c_byte
-        lgeos.GEOSPreparedIntersects.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedIntersects.restype = c_byte
+    lgeos.GEOSPreparedIntersects.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedCrosses.restype = c_byte
-        lgeos.GEOSPreparedCrosses.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedCrosses.restype = c_byte
+    lgeos.GEOSPreparedCrosses.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedWithin.restype = c_byte
-        lgeos.GEOSPreparedWithin.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedWithin.restype = c_byte
+    lgeos.GEOSPreparedWithin.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedContains.restype = c_byte
-        lgeos.GEOSPreparedContains.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedContains.restype = c_byte
+    lgeos.GEOSPreparedContains.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedContainsProperly.restype = c_byte
-        lgeos.GEOSPreparedContainsProperly.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedContainsProperly.restype = c_byte
+    lgeos.GEOSPreparedContainsProperly.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedOverlaps.restype = c_byte
-        lgeos.GEOSPreparedOverlaps.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedOverlaps.restype = c_byte
+    lgeos.GEOSPreparedOverlaps.argtypes = [c_void_p, c_void_p]
 
-        lgeos.GEOSPreparedCovers.restype = c_byte
-        lgeos.GEOSPreparedCovers.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSPreparedCovers.restype = c_byte
+    lgeos.GEOSPreparedCovers.argtypes = [c_void_p, c_void_p]
 
     # Geometry info
 
@@ -411,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
@@ -422,9 +402,8 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSDistance.restype = c_int
     lgeos.GEOSDistance.argtypes = [c_void_p, c_void_p, c_void_p]
 
-    if geos_version >= (3, 2, 0):
-        lgeos.GEOSHausdorffDistance.restype = c_int
-        lgeos.GEOSHausdorffDistance.argtypes = [c_void_p, c_void_p, c_void_p]
+    lgeos.GEOSHausdorffDistance.restype = c_int
+    lgeos.GEOSHausdorffDistance.argtypes = [c_void_p, c_void_p, c_void_p]
 
     # Reader and Writer APIs
 
@@ -446,22 +425,20 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSWKTWriter_write.restype = allocated_c_char_p
     lgeos.GEOSWKTWriter_write.argtypes = [c_void_p, c_void_p]
 
-    if geos_version >= (3, 3, 0):
-
-        lgeos.GEOSWKTWriter_setTrim.restype = None
-        lgeos.GEOSWKTWriter_setTrim.argtypes = [c_void_p, c_int]
+    lgeos.GEOSWKTWriter_setTrim.restype = None
+    lgeos.GEOSWKTWriter_setTrim.argtypes = [c_void_p, c_int]
 
-        lgeos.GEOSWKTWriter_setRoundingPrecision.restype = None
-        lgeos.GEOSWKTWriter_setRoundingPrecision.argtypes = [c_void_p, c_int]
+    lgeos.GEOSWKTWriter_setRoundingPrecision.restype = None
+    lgeos.GEOSWKTWriter_setRoundingPrecision.argtypes = [c_void_p, c_int]
 
-        lgeos.GEOSWKTWriter_setOutputDimension.restype = None
-        lgeos.GEOSWKTWriter_setOutputDimension.argtypes = [c_void_p, c_int]
+    lgeos.GEOSWKTWriter_setOutputDimension.restype = None
+    lgeos.GEOSWKTWriter_setOutputDimension.argtypes = [c_void_p, c_int]
 
-        lgeos.GEOSWKTWriter_getOutputDimension.restype = c_int
-        lgeos.GEOSWKTWriter_getOutputDimension.argtypes = [c_void_p]
+    lgeos.GEOSWKTWriter_getOutputDimension.restype = c_int
+    lgeos.GEOSWKTWriter_getOutputDimension.argtypes = [c_void_p]
 
-        lgeos.GEOSWKTWriter_setOld3D.restype = None
-        lgeos.GEOSWKTWriter_setOld3D.argtypes = [c_void_p, c_int]
+    lgeos.GEOSWKTWriter_setOld3D.restype = None
+    lgeos.GEOSWKTWriter_setOld3D.argtypes = [c_void_p, c_int]
 
     lgeos.GEOSWKBReader_create.restype = c_void_p
     lgeos.GEOSWKBReader_create.argtypes = []
@@ -505,16 +482,14 @@ def prototype(lgeos, geos_version):
     lgeos.GEOSWKBWriter_setIncludeSRID.restype = None
     lgeos.GEOSWKBWriter_setIncludeSRID.argtypes = [c_void_p, c_int]
 
-    if geos_version >= (3, 1, 1):
-        lgeos.GEOSFree.restype = None
-        lgeos.GEOSFree.argtypes = [c_void_p]
+    lgeos.GEOSFree.restype = None
+    lgeos.GEOSFree.argtypes = [c_void_p]
 
-    if geos_version >= (3, 3, 0):
-        lgeos.GEOSSnap.restype = c_void_p
-        lgeos.GEOSSnap.argtypes = [c_void_p, c_void_p, c_double]
+    lgeos.GEOSSnap.restype = c_void_p
+    lgeos.GEOSSnap.argtypes = [c_void_p, c_void_p, c_double]
 
-        lgeos.GEOSSharedPaths.restype = c_void_p
-        lgeos.GEOSSharedPaths.argtypes = [c_void_p, c_void_p]
+    lgeos.GEOSSharedPaths.restype = c_void_p
+    lgeos.GEOSSharedPaths.argtypes = [c_void_p, c_void_p]
 
     if geos_version >= (3, 4, 0):
         lgeos.GEOSNearestPoints.restype = c_void_p
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index ea8ee68..b97de9a 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -616,10 +616,6 @@ class BaseGeometry(object):
         if cap_style == CAP_STYLE.round and join_style == JOIN_STYLE.round:
             return geom_factory(self.impl['buffer'](self, distance, res))
 
-        if 'buffer_with_style' not in self.impl:
-            raise NotImplementedError("Styled buffering not available for "
-                                      "GEOS versions < 3.2.")
-
         return geom_factory(self.impl['buffer_with_style'](self, distance, res,
                                                            cap_style,
                                                            join_style,
@@ -641,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
     # -----------------
 
@@ -874,6 +891,11 @@ class BaseMultipartGeometry(BaseGeometry):
             return iter([])
 
     def __len__(self):
+        warn(
+            "__len__ for multi-part geometries is deprecated and will be removed in "
+            "Shapely 2.0. Check the length of the `geoms` property instead to get the "
+            " number of parts of a multi-part geometry.",
+            ShapelyDeprecationWarning, stacklevel=2)
         if not self.is_empty:
             return len(self.geoms)
         else:
@@ -898,7 +920,7 @@ class BaseMultipartGeometry(BaseGeometry):
     def __eq__(self, other):
         return (
             type(other) == type(self) and
-            len(self) == len(other) and
+            len(self.geoms) == len(other.geoms) and
             all(x == y for x, y in zip(self.geoms, other.geoms))
         )
 
diff --git a/shapely/geometry/multipoint.py b/shapely/geometry/multipoint.py
index bf13f85..6464e9f 100644
--- a/shapely/geometry/multipoint.py
+++ b/shapely/geometry/multipoint.py
@@ -48,7 +48,7 @@ class MultiPoint(BaseMultipartGeometry):
         """
         super(MultiPoint, self).__init__()
 
-        if points is None or len(points) == 0:
+        if points is None or (not isinstance(points, MultiPoint) and len(points) == 0):
             # allow creation of empty multipoints, to support unpickling
             pass
         else:
diff --git a/shapely/geos.py b/shapely/geos.py
index 8cf4b70..d7e766f 100644
--- a/shapely/geos.py
+++ b/shapely/geos.py
@@ -62,18 +62,27 @@ def exists_conda_env():
 
 
 if sys.platform.startswith('linux'):
-    # Test to see if we have a wheel repaired by 'auditwheel' containing its
-    # own libgeos_c
+    # Test to see if we have a wheel repaired by auditwheel which contains its
+    # own libgeos_c. Note: auditwheel 3.1 changed the location of libs.
     geos_whl_so = glob.glob(
+        os.path.abspath(os.path.join(os.path.dirname(__file__), ".libs/libgeos*.so*"))
+    ) or glob.glob(
         os.path.abspath(
             os.path.join(
-                os.path.dirname(__file__), ".libs/libgeos_c-*.so.*"
+                os.path.dirname(__file__), "..", "Shapely.libs", "libgeos*.so*"
             )
         )
     )
-    if len(geos_whl_so) == 1:
-        _lgeos = CDLL(geos_whl_so[0])
+
+    if len(geos_whl_so) > 0:
+        # We have observed problems with CDLL of libgeos_c not automatically
+        # loading the sibling c++ library since the change made by auditwheel
+        # 3.1, so we explicitly load them both.
+        geos_whl_so = sorted(geos_whl_so)
+        CDLL(geos_whl_so[0])
+        _lgeos = CDLL(geos_whl_so[-1])
         LOG.debug("Found GEOS DLL: %r, using it.", _lgeos)
+
     elif hasattr(sys, 'frozen'):
         geos_pyinstaller_so = glob.glob(os.path.join(sys.prefix, 'libgeos_c-*.so.*'))
         if len(geos_pyinstaller_so) == 1:
@@ -87,7 +96,7 @@ if sys.platform.startswith('linux'):
             'libgeos_c.so.1',
             'libgeos_c.so',
         ]
-        _lgeos = load_dll('geos_c', fallbacks=alt_paths)
+        _lgeos = load_dll('libgeos_c', fallbacks=alt_paths)
     # Necessary for environments with only libc.musl
     c_alt_paths = [
         'libc.musl-x86_64.so.1'
@@ -201,40 +210,38 @@ def _geos_version():
 geos_version_string, geos_version, geos_capi_version = _geos_version()
 
 
-# If we have the new interface, then record a baseline so that we know what
-# additional functions are declared in ctypes_declarations.
-if geos_version >= (3, 1, 0):
-    start_set = set(_lgeos.__dict__)
+# Record a baseline so that we know what additional functions are declared
+# in ctypes_declarations.
+start_set = set(_lgeos.__dict__)
 
 # Apply prototypes for the libgeos_c functions
 prototype(_lgeos, geos_version)
 
-# If we have the new interface, automatically detect all function
-# declarations, and declare their re-entrant counterpart.
-if geos_version >= (3, 1, 0):
-    end_set = set(_lgeos.__dict__)
-    new_func_names = end_set - start_set
-
-    for func_name in new_func_names:
-        new_func_name = "%s_r" % func_name
-        if hasattr(_lgeos, new_func_name):
-            new_func = getattr(_lgeos, new_func_name)
-            old_func = getattr(_lgeos, func_name)
-            new_func.restype = old_func.restype
-            if old_func.argtypes is None:
-                # Handle functions that didn't take an argument before,
-                # finishGEOS.
-                new_func.argtypes = [c_void_p]
-            else:
-                new_func.argtypes = [c_void_p] + list(old_func.argtypes)
-            if old_func.errcheck is not None:
-                new_func.errcheck = old_func.errcheck
+# Automatically detect all function declarations, and declare their
+# re-entrant counterpart.
+end_set = set(_lgeos.__dict__)
+new_func_names = end_set - start_set
+
+for func_name in new_func_names:
+    new_func_name = "%s_r" % func_name
+    if hasattr(_lgeos, new_func_name):
+        new_func = getattr(_lgeos, new_func_name)
+        old_func = getattr(_lgeos, func_name)
+        new_func.restype = old_func.restype
+        if old_func.argtypes is None:
+            # Handle functions that didn't take an argument before,
+            # finishGEOS.
+            new_func.argtypes = [c_void_p]
+        else:
+            new_func.argtypes = [c_void_p] + list(old_func.argtypes)
+        if old_func.errcheck is not None:
+            new_func.errcheck = old_func.errcheck
 
-    # Handle special case.
-    _lgeos.initGEOS_r.restype = c_void_p
-    _lgeos.initGEOS_r.argtypes = \
-        [EXCEPTION_HANDLER_FUNCTYPE, EXCEPTION_HANDLER_FUNCTYPE]
-    _lgeos.finishGEOS_r.argtypes = [c_void_p]
+# Handle special case.
+_lgeos.initGEOS_r.restype = c_void_p
+_lgeos.initGEOS_r.argtypes = \
+    [EXCEPTION_HANDLER_FUNCTYPE, EXCEPTION_HANDLER_FUNCTYPE]
+_lgeos.finishGEOS_r.argtypes = [c_void_p]
 
 
 def make_logging_callback(func):
@@ -300,71 +307,63 @@ class WKTWriter(object):
     _writer = None
 
     # Establish default output settings
-    defaults = {}
-
-    if geos_version >= (3, 3, 0):
-
-        defaults['trim'] = True
-        defaults['output_dimension'] = 3
-
-        # GEOS' defaults for methods without "get"
-        _trim = False
-        _rounding_precision = -1
-        _old_3d = False
-
-        @property
-        def trim(self):
-            """Trimming of unnecessary decimals (default: True)"""
-            return getattr(self, '_trim')
-
-        @trim.setter
-        def trim(self, value):
-            self._trim = bool(value)
-            self._lgeos.GEOSWKTWriter_setTrim(self._writer, self._trim)
-
-        @property
-        def rounding_precision(self):
-            """Rounding precision when writing the WKT.
-            A precision of -1 (default) disables it."""
-            return getattr(self, '_rounding_precision')
-
-        @rounding_precision.setter
-        def rounding_precision(self, value):
-            self._rounding_precision = int(value)
-            self._lgeos.GEOSWKTWriter_setRoundingPrecision(
-                self._writer, self._rounding_precision)
-
-        @property
-        def output_dimension(self):
-            """Output dimension, either 2 or 3 (default)"""
-            return self._lgeos.GEOSWKTWriter_getOutputDimension(
-                self._writer)
-
-        @output_dimension.setter
-        def output_dimension(self, value):
-            self._lgeos.GEOSWKTWriter_setOutputDimension(
-                self._writer, int(value))
-
-        @property
-        def old_3d(self):
-            """Show older style for 3D WKT, without 'Z' (default: False)"""
-            return getattr(self, '_old_3d')
-
-        @old_3d.setter
-        def old_3d(self, value):
-            self._old_3d = bool(value)
-            self._lgeos.GEOSWKTWriter_setOld3D(self._writer, self._old_3d)
+    defaults = {
+        'trim': True,
+        'output_dimension': 3,
+    }
 
-    def __init__(self, lgeos, **settings):
-        """Create WKT Writer
+    # GEOS' defaults for methods without "get"
+    _trim = False
+    _rounding_precision = -1
+    _old_3d = False
+
+    @property
+    def trim(self):
+        """Trimming of unnecessary decimals (default: True)"""
+        return getattr(self, '_trim')
+
+    @trim.setter
+    def trim(self, value):
+        self._trim = bool(value)
+        self._lgeos.GEOSWKTWriter_setTrim(self._writer, self._trim)
+
+    @property
+    def rounding_precision(self):
+        """Rounding precision when writing the WKT.
+        A precision of -1 (default) disables it."""
+        return getattr(self, '_rounding_precision')
+
+    @rounding_precision.setter
+    def rounding_precision(self, value):
+        self._rounding_precision = int(value)
+        self._lgeos.GEOSWKTWriter_setRoundingPrecision(
+            self._writer, self._rounding_precision)
+
+    @property
+    def output_dimension(self):
+        """Output dimension, either 2 or 3 (default)"""
+        return self._lgeos.GEOSWKTWriter_getOutputDimension(
+            self._writer)
+
+    @output_dimension.setter
+    def output_dimension(self, value):
+        self._lgeos.GEOSWKTWriter_setOutputDimension(
+            self._writer, int(value))
 
-        Note: writer defaults are set differently for GEOS 3.3.0 and up.
-        For example, with 'POINT Z (1 2 3)':
+    @property
+    def old_3d(self):
+        """Show older style for 3D WKT, without 'Z' (default: False)"""
+        return getattr(self, '_old_3d')
 
-            newer: POINT Z (1 2 3)
-            older: POINT (1.0000000000000000 2.0000000000000000)
+    @old_3d.setter
+    def old_3d(self, value):
+        self._old_3d = bool(value)
+        self._lgeos.GEOSWKTWriter_setOld3D(self._writer, self._old_3d)
 
-        The older formatting can be achieved for GEOS 3.3.0 and up by setting
+    def __init__(self, lgeos, **settings):
+        """Create WKT Writer
+
+        Note: older formatting before GEOS 3.3.0 can be achieved by setting
         the properties:
             trim = False
             output_dimension = 2
@@ -594,83 +593,14 @@ class LGEOSBase(threading.local):
             self.geos_handle = None
 
 
-class LGEOS300(LGEOSBase):
-    """Proxy for GEOS 3.0.0-CAPI-1.4.1
-    """
-    geos_version = (3, 0, 0)
-    geos_capi_version = (1, 4, 0)
-
-    def __init__(self, dll):
-        super(LGEOS300, self).__init__(dll)
-        self.geos_handle = self._lgeos.initGEOS(notice_h, error_h)
-        keys = list(self._lgeos.__dict__.keys())
-        for key in keys:
-            setattr(self, key, getattr(self._lgeos, key))
-        self.GEOSFree = self._lgeos.free
-        # Deprecated
-        self.GEOSGeomToWKB_buf.errcheck = errcheck_wkb
-        self.GEOSGeomToWKT.errcheck = errcheck_just_free
-        self.GEOSRelate.errcheck = errcheck_null_exception
-        for pred in (
-                self.GEOSDisjoint,
-                self.GEOSTouches,
-                self.GEOSIntersects,
-                self.GEOSCrosses,
-                self.GEOSWithin,
-                self.GEOSContains,
-                self.GEOSOverlaps,
-                self.GEOSEquals,
-                self.GEOSEqualsExact,
-                self.GEOSRelatePattern,
-                self.GEOSisEmpty,
-                self.GEOSisValid,
-                self.GEOSisSimple,
-                self.GEOSisRing,
-                self.GEOSHasZ):
-            pred.errcheck = errcheck_predicate
-
-        self.methods['area'] = self.GEOSArea
-        self.methods['boundary'] = self.GEOSBoundary
-        self.methods['buffer'] = self.GEOSBuffer
-        self.methods['centroid'] = self.GEOSGetCentroid
-        self.methods['representative_point'] = self.GEOSPointOnSurface
-        self.methods['convex_hull'] = self.GEOSConvexHull
-        self.methods['distance'] = self.GEOSDistance
-        self.methods['envelope'] = self.GEOSEnvelope
-        self.methods['length'] = self.GEOSLength
-        self.methods['has_z'] = self.GEOSHasZ
-        self.methods['is_empty'] = self.GEOSisEmpty
-        self.methods['is_ring'] = self.GEOSisRing
-        self.methods['is_simple'] = self.GEOSisSimple
-        self.methods['is_valid'] = self.GEOSisValid
-        self.methods['disjoint'] = self.GEOSDisjoint
-        self.methods['touches'] = self.GEOSTouches
-        self.methods['intersects'] = self.GEOSIntersects
-        self.methods['crosses'] = self.GEOSCrosses
-        self.methods['within'] = self.GEOSWithin
-        self.methods['contains'] = self.GEOSContains
-        self.methods['overlaps'] = self.GEOSOverlaps
-        self.methods['equals'] = self.GEOSEquals
-        self.methods['equals_exact'] = self.GEOSEqualsExact
-        self.methods['relate'] = self.GEOSRelate
-        self.methods['difference'] = self.GEOSDifference
-        self.methods['symmetric_difference'] = self.GEOSSymDifference
-        self.methods['union'] = self.GEOSUnion
-        self.methods['intersection'] = self.GEOSIntersection
-        self.methods['relate_pattern'] = self.GEOSRelatePattern
-        self.methods['simplify'] = self.GEOSSimplify
-        self.methods['topology_preserve_simplify'] = \
-            self.GEOSTopologyPreserveSimplify
-
-
-class LGEOS310(LGEOSBase):
-    """Proxy for GEOS 3.1.0-CAPI-1.5.0
+class LGEOS330(LGEOSBase):
+    """Proxy for GEOS 3.3.0-CAPI-1.7.0
     """
-    geos_version = (3, 1, 0)
-    geos_capi_version = (1, 5, 0)
+    geos_version = (3, 3, 0)
+    geos_capi_version = (1, 7, 0)
 
     def __init__(self, dll):
-        super(LGEOS310, self).__init__(dll)
+        super(LGEOS330, self).__init__(dll)
         self.geos_handle = self._lgeos.initGEOS_r(notice_h, error_h)
         keys = list(self._lgeos.__dict__.keys())
         for key in [x for x in keys if not x.endswith('_r')]:
@@ -681,9 +611,15 @@ class LGEOS310(LGEOSBase):
                 setattr(self, key, attr)
             else:
                 setattr(self, key, getattr(self._lgeos, key))
-        if not hasattr(self, 'GEOSFree'):
-            # GEOS < 3.1.1
-            self.GEOSFree = self._lgeos.free
+
+        # GEOS 3.3.8 from homebrew has, but doesn't advertise
+        # GEOSPolygonize_full. We patch it in explicitly here.
+        key = 'GEOSPolygonize_full'
+        func = getattr(self._lgeos, key + '_r')
+        attr = partial(func, self.geos_handle)
+        attr.__name__ = func.__name__
+        setattr(self, key, attr)
+
         # Deprecated
         self.GEOSGeomToWKB_buf.func.errcheck = errcheck_wkb
         self.GEOSGeomToWKT.func.errcheck = errcheck_just_free
@@ -713,7 +649,9 @@ class LGEOS310(LGEOSBase):
                 self.GEOSisValid,
                 self.GEOSisSimple,
                 self.GEOSisRing,
-                self.GEOSHasZ):
+                self.GEOSHasZ,
+                self.GEOSisClosed,
+                self.GEOSCoveredBy):
             pred.func.errcheck = errcheck_predicate
 
         self.GEOSisValidReason.func.errcheck = errcheck_just_free
@@ -761,70 +699,9 @@ class LGEOS310(LGEOSBase):
         self.methods['simplify'] = self.GEOSSimplify
         self.methods['topology_preserve_simplify'] = \
             self.GEOSTopologyPreserveSimplify
+        self.methods['normalize'] = self.GEOSNormalize
         self.methods['cascaded_union'] = self.GEOSUnionCascaded
 
-
-class LGEOS311(LGEOS310):
-    """Proxy for GEOS 3.1.1-CAPI-1.6.0
-    """
-    geos_version = (3, 1, 1)
-    geos_capi_version = (1, 6, 0)
-
-    def __init__(self, dll):
-        super(LGEOS311, self).__init__(dll)
-
-
-class LGEOS320(LGEOS311):
-    """Proxy for GEOS 3.2.0-CAPI-1.6.0
-    """
-    geos_version = (3, 2, 0)
-    geos_capi_version = (1, 6, 0)
-
-    def __init__(self, dll):
-        super(LGEOS320, self).__init__(dll)
-
-        if geos_version >= (3, 2, 0):
-
-            def parallel_offset(geom, distance, resolution=16, join_style=1,
-                                mitre_limit=5.0, side='right'):
-                side = side == 'left'
-                if distance < 0:
-                    distance = abs(distance)
-                    side = not side
-                return self.GEOSSingleSidedBuffer(
-                    geom, distance, resolution, join_style, mitre_limit, side)
-
-            self.methods['parallel_offset'] = parallel_offset
-
-        self.methods['project'] = self.GEOSProject
-        self.methods['project_normalized'] = self.GEOSProjectNormalized
-        self.methods['interpolate'] = self.GEOSInterpolate
-        self.methods['interpolate_normalized'] = \
-            self.GEOSInterpolateNormalized
-        self.methods['buffer_with_style'] = self.GEOSBufferWithStyle
-        self.methods['hausdorff_distance'] = self.GEOSHausdorffDistance
-
-
-class LGEOS330(LGEOS320):
-    """Proxy for GEOS 3.3.0-CAPI-1.7.0
-    """
-    geos_version = (3, 3, 0)
-    geos_capi_version = (1, 7, 0)
-
-    def __init__(self, dll):
-        super(LGEOS330, self).__init__(dll)
-
-        # GEOS 3.3.8 from homebrew has, but doesn't advertise
-        # GEOSPolygonize_full. We patch it in explicitly here.
-        key = 'GEOSPolygonize_full'
-        func = getattr(self._lgeos, key + '_r')
-        attr = partial(func, self.geos_handle)
-        attr.__name__ = func.__name__
-        setattr(self, key, attr)
-
-        for pred in (self.GEOSisClosed, self.GEOSCoveredBy):
-            pred.func.errcheck = errcheck_predicate
-
         def parallel_offset(geom, distance, resolution=16, join_style=1,
                             mitre_limit=5.0, side='right'):
             if side == 'right':
@@ -833,9 +710,16 @@ class LGEOS330(LGEOS320):
                 geom, distance, resolution, join_style, mitre_limit)
 
         self.methods['parallel_offset'] = parallel_offset
+        self.methods['project'] = self.GEOSProject
+        self.methods['project_normalized'] = self.GEOSProjectNormalized
+        self.methods['interpolate'] = self.GEOSInterpolate
+        self.methods['interpolate_normalized'] = \
+            self.GEOSInterpolateNormalized
+        self.methods['buffer_with_style'] = self.GEOSBufferWithStyle
+        self.methods['hausdorff_distance'] = self.GEOSHausdorffDistance
         self.methods['unary_union'] = self.GEOSUnaryUnion
-        self.methods['is_closed'] = self.GEOSisClosed
         self.methods['cascaded_union'] = self.methods['unary_union']
+        self.methods['is_closed'] = self.GEOSisClosed
         self.methods['snap'] = self.GEOSSnap
         self.methods['shared_paths'] = self.GEOSSharedPaths
         self.methods['buffer_with_params'] = self.GEOSBufferWithParams
@@ -898,14 +782,8 @@ elif geos_version >= (3, 4, 0):
     L = LGEOS340
 elif geos_version >= (3, 3, 0):
     L = LGEOS330
-elif geos_version >= (3, 2, 0):
-    L = LGEOS320
-elif geos_version >= (3, 1, 1):
-    L = LGEOS311
-elif geos_version >= (3, 1, 0):
-    L = LGEOS310
 else:
-    L = LGEOS300
+    raise ValueError('unexpected geos_version: ' + str(geos_version))
 
 lgeos = L(_lgeos)
 
diff --git a/shapely/impl.py b/shapely/impl.py
index 1a7481b..9b206e9 100644
--- a/shapely/impl.py
+++ b/shapely/impl.py
@@ -73,7 +73,7 @@ class GEOSImpl(BaseImpl):
             lgeos.geos_capi_version,)
 
 
-IMPL300 = {
+IMPL330 = {
     'area': (UnaryRealProperty, 'area'),
     'distance': (BinaryRealProperty, 'distance'),
     'length': (UnaryRealProperty, 'length'),
@@ -85,6 +85,10 @@ IMPL300 = {
     'envelope': (UnaryTopologicalOp, 'envelope'),
     'convex_hull': (UnaryTopologicalOp, 'convex_hull'),
     'buffer': (UnaryTopologicalOp, 'buffer'),
+    'simplify': (UnaryTopologicalOp, 'simplify'),
+    'topology_preserve_simplify':
+        (UnaryTopologicalOp, 'topology_preserve_simplify'),
+    'normalize': (UnaryTopologicalOp, 'normalize'),
     #
     'difference': (BinaryTopologicalOp, 'difference'),
     'intersection': (BinaryTopologicalOp, 'intersection'),
@@ -96,6 +100,7 @@ IMPL300 = {
     'is_ring': (UnaryPredicate, 'is_ring'),
     'is_simple': (UnaryPredicate, 'is_simple'),
     'is_valid': (UnaryPredicate, 'is_valid'),
+    'is_closed': (UnaryPredicate, 'is_closed'),
     #
     'relate': (BinaryPredicate, 'relate'),
     'contains': (BinaryPredicate, 'contains'),
@@ -107,17 +112,9 @@ IMPL300 = {
     'touches': (BinaryPredicate, 'touches'),
     'within': (BinaryPredicate, 'within'),
     'covers': (BinaryPredicate, 'covers'),
+    'covered_by': (BinaryPredicate, 'covered_by'),
     'equals_exact': (BinaryPredicate, 'equals_exact'),
     'relate_pattern': (BinaryPredicate, 'relate_pattern'),
-
-    # First pure Python implementation
-    'is_ccw': (cga.is_ccw_impl, 'is_ccw'),
-    }
-
-IMPL310 = {
-    'simplify': (UnaryTopologicalOp, 'simplify'),
-    'topology_preserve_simplify':
-        (UnaryTopologicalOp, 'topology_preserve_simplify'),
     'prepared_disjoint': (BinaryPredicate, 'prepared_disjoint'),
     'prepared_touches': (BinaryPredicate, 'prepared_touches'),
     'prepared_crosses': (BinaryPredicate, 'prepared_crosses'),
@@ -128,26 +125,19 @@ IMPL310 = {
     'prepared_contains_properly':
         (BinaryPredicate, 'prepared_contains_properly'),
     'prepared_covers': (BinaryPredicate, 'prepared_covers'),
-    }
-
-IMPL311 = {
-    }
 
-IMPL320 = {
     'parallel_offset': (UnaryTopologicalOp, 'parallel_offset'),
     'project_normalized': (ProjectOp, 'project_normalized'),
     'project': (ProjectOp, 'project'),
     'interpolate_normalized': (InterpolateOp, 'interpolate_normalized'),
     'interpolate': (InterpolateOp, 'interpolate'),
     'buffer_with_style': (UnaryTopologicalOp, 'buffer_with_style'),
+    'buffer_with_params': (UnaryTopologicalOp, 'buffer_with_params'),
     'hausdorff_distance': (BinaryRealProperty, 'hausdorff_distance'),
-    }
 
-IMPL330 = {
-    'is_closed': (UnaryPredicate, 'is_closed'),
-    'buffer_with_params': (UnaryTopologicalOp, 'buffer_with_params'),
-    'covered_by': (BinaryPredicate, 'covered_by')
-}
+    # First pure Python implementation
+    'is_ccw': (cga.is_ccw_impl, 'is_ccw'),
+    }
 
 IMPL360 = {
     'minimum_clearance': (UnaryRealProperty, 'minimum_clearance')
@@ -156,15 +146,7 @@ IMPL360 = {
 def impl_items(defs):
     return [(k, v[0](v[1])) for k, v in list(defs.items())]
 
-imp = GEOSImpl(dict(impl_items(IMPL300)))
-if lgeos.geos_version >= (3, 1, 0):
-    imp.update(impl_items(IMPL310))
-if lgeos.geos_version >= (3, 1, 1):
-    imp.update(impl_items(IMPL311))
-if lgeos.geos_version >= (3, 2, 0):
-    imp.update(impl_items(IMPL320))
-if lgeos.geos_version >= (3, 3, 0):
-    imp.update(impl_items(IMPL330))
+imp = GEOSImpl(dict(impl_items(IMPL330)))
 if lgeos.geos_version >= (3, 6, 0):
     imp.update(impl_items(IMPL360))
 
diff --git a/shapely/iterops.py b/shapely/iterops.py
index a0c7deb..4011e6d 100644
--- a/shapely/iterops.py
+++ b/shapely/iterops.py
@@ -1,7 +1,9 @@
 """
 Iterative forms of operations
 """
+import warnings
 
+from shapely.errors import ShapelyDeprecationWarning
 from shapely.topology import Delegating
 
 
@@ -11,6 +13,10 @@ class IterOp(Delegating):
     """
 
     def __call__(self, context, iterator, value=True):
+        warnings.warn(
+            "The '{0}' function is deprecated and will be removed in "
+            "Shapely 2.0".format(self._name),
+            ShapelyDeprecationWarning, stacklevel=2)
         if context._geom is None:
             raise ValueError("Null geometry supports no operations")
         for item in iterator:
diff --git a/shapely/ops.py b/shapely/ops.py
index b9dcbb0..ccbbc86 100644
--- a/shapely/ops.py
+++ b/shapely/ops.py
@@ -129,9 +129,9 @@ class CollectionOperator(object):
             "Use 'unary_union()' instead.",
             ShapelyDeprecationWarning, stacklevel=2)
         try:
-            L = len(geoms)
             if isinstance(geoms, BaseMultipartGeometry):
                 geoms = geoms.geoms
+            L = len(geoms)
         except TypeError:
             geoms = [geoms]
             L = 1
@@ -148,9 +148,9 @@ class CollectionOperator(object):
         prefered method for dissolving many polygons.
         """
         try:
-            L = len(geoms)
             if isinstance(geoms, BaseMultipartGeometry):
                 geoms = geoms.geoms
+            L = len(geoms)
         except TypeError:
             geoms = [geoms]
             L = 1
@@ -429,14 +429,20 @@ class SplitOp(object):
         assert(isinstance(line, LineString))
         assert(isinstance(splitter, LineString) or isinstance(splitter, MultiLineString))
 
-        if splitter.crosses(line):
-            # The lines cross --> return multilinestring from the split
-            return line.difference(splitter)
-        elif splitter.relate_pattern(line, '1********'):
+        # |    s\l   | Interior | Boundary | Exterior |
+        # |----------|----------|----------|----------|
+        # | Interior |  0 or F  |    *     |    *     |   At least one of these two must be 0
+        # | Boundary |  0 or F  |    *     |    *     |   So either '0********' or '[0F]**0*****'
+        # | Exterior |    *     |    *     |    *     |   No overlapping interiors ('1********')
+        relation = splitter.relate(line)
+        if relation[0] == '1':
             # The lines overlap at some segment (linear intersection of interiors)
             raise ValueError('Input geometry segment overlaps with the splitter.')
+        elif relation[0] == '0' or relation[3] == '0':
+            # The splitter crosses or touches the line's interior --> return multilinestring from the split
+            return line.difference(splitter)
         else:
-            # The lines do not cross --> return collection with identity line
+            # The splitter does not cross or touch the line's interior --> return collection with identity line
             return [line]
 
     @staticmethod
@@ -560,17 +566,66 @@ split = SplitOp.split
 
 
 def substring(geom, start_dist, end_dist, normalized=False):
-    """Return a line segment between specified distances along a linear geometry.
+    """Return a line segment between specified distances along a LineString
 
     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 distances 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.
+
+    If the start distance equals the end distance, a Point is returned.
+
+    If the start distance is actually beyond the end distance, then the
+    reversed substring is returned such that the start distance is
+    at the first coordinate.
+
+    Parameters
+    ----------
+    geom : LineString
+        The geometry to get a substring of.
+    start_dist : float
+        The distance along `geom` of the start of the substring.
+    end_dist : float
+        The distance along `geom` of the end of the substring.
+    normalized : bool, False
+        Whether the distance parameters are interpreted as a
+        fraction of the geometry's length.
+
+    Returns
+    -------
+    Union[Point, LineString]
+        The substring between `start_dist` and `end_dist` or a Point
+        if they are at the same location.
+
+    Raises
+    ------
+    TypeError
+        If `geom` is not a LineString.
+
+    Examples
+    --------
+    >>> from shapely.geometry import LineString
+    >>> from shapely.ops import substring
+    >>> ls = LineString((i, 0) for i in range(6))
+    >>> ls.wkt
+    'LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)'
+    >>> substring(ls, start_dist=1, end_dist=3).wkt
+    'LINESTRING (1 0, 2 0, 3 0)'
+    >>> substring(ls, start_dist=3, end_dist=1).wkt
+    'LINESTRING (3 0, 2 0, 1 0)'
+    >>> substring(ls, start_dist=1, end_dist=-3).wkt
+    'LINESTRING (1 0, 2 0)'
+    >>> substring(ls, start_dist=0.2, end_dist=-0.6, normalized=True).wkt
+    'LINESTRING (1 0, 2 0)'
+
+    Returning a `Point` when `start_dist` and `end_dist` are at the
+    same location.
+
+    >>> substring(ls, 2.5, -2.5).wkt
+    'POINT (2.5 0)'
     """
 
-    assert(isinstance(geom, LineString))
+    if not isinstance(geom, LineString):
+        raise TypeError("Can only calculate a substring of LineString geometries. A %s was provided." % geom.type)
 
     # Filter out cases in which to return a point
     if start_dist == end_dist:
@@ -584,36 +639,52 @@ def substring(geom, start_dist, end_dist, normalized=False):
     elif normalized and -start_dist >= 1 and -end_dist >= 1:
         return geom.interpolate(0, normalized)
 
-    start_point = geom.interpolate(start_dist, normalized)
-    end_point = geom.interpolate(end_dist, normalized)
-
-    min_dist = min(start_dist, end_dist)
-    max_dist = max(start_dist, end_dist)
     if normalized:
-        min_dist *= geom.length
-        max_dist *= geom.length
+        start_dist *= geom.length
+        end_dist *= geom.length
 
-    if start_dist < end_dist:
-        vertex_list = [(start_point.x, start_point.y)]
-    else:
+    # Filter out cases where distances meet at a middle point from opposite ends.
+    if start_dist < 0 < end_dist and abs(start_dist) + end_dist == geom.length:
+        return geom.interpolate(end_dist)
+    elif end_dist < 0 < start_dist and abs(end_dist) + start_dist == geom.length:
+        return geom.interpolate(start_dist)
+
+    start_point = geom.interpolate(start_dist)
+    end_point = geom.interpolate(end_dist)
+
+    if start_dist < 0:
+        start_dist = geom.length + start_dist  # Values may still be negative,
+    if end_dist < 0:                           # but only in the out-of-range
+        end_dist = geom.length + end_dist      # sense, not the wrap-around sense.
+
+    reverse = start_dist > end_dist
+    if reverse:
+        start_dist, end_dist = end_dist, start_dist
+
+    if start_dist < 0:
+        start_dist = 0  # to avoid duplicating the first vertex
+
+    if reverse:
         vertex_list = [(end_point.x, end_point.y)]
+    else:
+        vertex_list = [(start_point.x, start_point.y)]
 
     coords = list(geom.coords)
     current_distance = 0
     for p1, p2 in zip(coords, coords[1:]):
-        if min_dist < current_distance < max_dist:
+        if start_dist < current_distance < end_dist:
             vertex_list.append(p1)
-        elif current_distance >= max_dist:
+        elif current_distance >= end_dist:
             break
 
         current_distance += ((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2) ** 0.5
 
-    if start_dist < end_dist:
-        vertex_list.append((end_point.x, end_point.y))
-    else:
+    if reverse:
         vertex_list.append((start_point.x, start_point.y))
         # reverse direction result
         vertex_list = reversed(vertex_list)
+    else:
+        vertex_list.append((end_point.x, end_point.y))
 
     return LineString(vertex_list)
 
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))
 
diff --git a/shapely/topology.py b/shapely/topology.py
index 1618046..d4e56b1 100644
--- a/shapely/topology.py
+++ b/shapely/topology.py
@@ -23,6 +23,7 @@ class Validating(object):
 class Delegating(Validating):
 
     def __init__(self, name):
+        self._name = name
         self.fn = lgeos.methods[name]
 
     def _check_topology(self, err, *geoms):
diff --git a/shapely/wkt.py b/shapely/wkt.py
index cf98671..b465ddd 100644
--- a/shapely/wkt.py
+++ b/shapely/wkt.py
@@ -49,10 +49,10 @@ def dumps(ob, trim=False, **kw):
         A geometry object of any type to be dumped to WKT.
     trim : bool, default False
         Remove excess decimals from the WKT.
-    rounding_precision (GEOS 3.3+) : int
+    rounding_precision : int
         Round output to the specified number of digits.
         Default behavior returns full precision.
-    output_dimension (GEOS 3.3+): int, default 3
+    output_dimension : int, default 3
         Force removal of dimensions above the one specified.
 
     Returns
@@ -74,10 +74,10 @@ def dump(ob, fp, **settings):
         A file-like object which implements a `write` method.
     trim : bool, default False
         Remove excess decimals from the WKT.
-    rounding_precision (GEOS 3.3+) : int
+    rounding_precision : int
         Round output to the specified number of digits.
         Default behavior returns full precision.
-    output_dimension (GEOS 3.3+): int, default 3
+    output_dimension : int, default 3
         Force removal of dimensions above the one specified.
 
     Returns


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..628ef2d
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,174 @@
+name: Tests (Linux)
+
+on: [push, pull_request]
+
+jobs:
+  test:
+    name: Python ${{ matrix.python }} (GEOS ${{ matrix.geos }}, speedups ${{ matrix.speedups }}, numpy ${{ matrix.numpy || 'not installed' }})
+    # runs-on: ubuntu-latest
+    runs-on: ubuntu-16.04
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          # 2015
+          - python: 3.5
+            geos: 3.5.2
+            numpy: 1.10.4
+            speedups: 1
+          - python: 3.5
+            geos: 3.5.2
+            numpy: 1.10.4
+            speedups: 0
+          - python: 3.5
+            geos: 3.5.2
+            speedups: 0
+          # 2017
+          - python: 3.6
+            geos: 3.6.4
+            numpy: 1.13.3
+            speeedups: 1
+          - python: 3.6
+            geos: 3.6.4
+            numpy: 1.13.3
+            speeedups: 0
+          - python: 3.6
+            geos: 3.6.4
+            speeedups: 0
+          # 2018
+          - python: 3.7
+            geos: 3.7.3
+            numpy: 1.15.4
+            speedups: 1
+          - python: 3.7
+            geos: 3.7.3
+            numpy: 1.15.4
+            speedups: 0
+          - python: 3.7
+            geos: 3.7.3
+            speedups: 0
+          # 2019
+          - python: 3.8
+            geos: 3.8.1
+            numpy: 1.17.5
+            speedups: 1
+          - python: 3.8
+            geos: 3.8.1
+            numpy: 1.17.5
+            speedups: 0
+          - python: 3.8
+            geos: 3.8.1
+            speedups: 0
+          # 2020
+          - python: 3.9
+            geos: 3.9.1
+            numpy: 1.19.5
+            speedups: 1
+          - python: 3.9
+            geos: 3.9.1
+            numpy: 1.19.5
+            speedups: 0
+          - python: 3.9
+            geos: 3.9.1
+            speedups: 0
+          # dev
+          - python: 3.9
+            geos: master
+            speedups: 1
+
+    env:
+      GEOS_VERSION: ${{ matrix.geos }}
+      GEOS_INSTALL: ${{ github.workspace }}/geosinstall/geos-${{ matrix.geos }}
+      SPEEDUPS: ${{ matrix.speedups }}
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Set up Python ${{ matrix.python }}
+        uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python }}
+
+      - name: Cache GEOS and pip packages
+        uses: actions/cache@v2
+        with:
+          key: ${{ runner.os }}-geos-${{ matrix.geos }}
+          path: |
+            ~/.cache/pip
+            ${{ github.workspace }}/geosinstall
+
+      - name: Install GEOS
+        shell: bash
+        run: |
+          ./ci/install_geos.sh
+          cd ${{ github.workspace }}
+
+      - name: Install python dependencies
+        shell: bash
+        run: |
+          pip install --disable-pip-version-check --upgrade pip
+          pip install --upgrade wheel
+          if [ "$GEOS_VERSION" = "master" ]; then
+            pip install --upgrade --pre Cython numpy;
+          elif [ "$SPEEDUPS" == "1" ]; then
+            pip install --install-option="--no-cython-compile" cython;
+            pip install --upgrade numpy==${{ matrix.numpy }};
+          elif [ -n "${{ matrix.numpy }}" ]; then
+            pip install --upgrade numpy==${{ matrix.numpy }};
+          else
+            pip uninstall --yes numpy;
+          fi
+          pip install --upgrade coveralls pytest-cov pytest
+          pip list
+
+      - name: Set environment variables
+        shell: bash
+        run: |
+          echo "LD_LIBRARY_PATH=$GEOS_INSTALL/lib" >> $GITHUB_ENV
+          echo $GEOS_INSTALL/bin >> $GITHUB_PATH
+
+      - name: Build and install Shapely
+        shell: bash
+        run: |
+          pip install -v -e .
+
+      - name: Overview of the Python environment (pip list)
+        shell: bash
+        run: |
+          pip list
+
+      - name: Run tests
+        shell: bash
+        continue-on-error: ${{ matrix.geos == 'master' }}
+        run: |
+          if [ "$SPEEDUPS" == "1" ]; then SPEEDUPS_FLAG=--with-speedups; else SPEEDUPS_FLAG=--without-speedups; fi
+          python -c "from shapely import geos; print(geos.geos_version_string)"
+          python -c "from shapely import speedups; print(speedups.enabled)"
+          python -m pytest --cov shapely --cov-report term-missing "${SPEEDUPS_FLAG}"
+
+      - name: Run doctests
+        if: ${{ matrix.python == '3.8' && matrix.numpy && matrix.speedups == '1' }}
+        shell: bash
+        run: |
+          python -m pytest shapely --doctest-modules
+
+      - name: Upload coverage
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COVERALLS_PARALLEL: true
+        shell: bash
+        run: |
+          coveralls --service=github || echo "!! intermittent coveralls failure"
+
+  coveralls:
+    name: Indicate completion to coveralls.io
+    needs: test
+    runs-on: ubuntu-latest
+    container: python:3-slim
+    steps:
+    - name: Finished
+      run: |
+        pip3 install --upgrade coveralls
+        coveralls --finish || echo "!! intermittent coveralls failure"
+      env:
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/tests/__init__.py b/tests/__init__.py
index 9ea07db..d7bcf52 100755
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -17,16 +17,12 @@ except ImportError:
     numpy = False
     numpy_version = 'not available'
 
-# Show some diagnostic information; handy for Travis CI
+# Show some diagnostic information; handy for CI
 print('Python version: ' + sys.version.replace('\n', ' '))
 print('GEOS version: ' + geos_version_string)
 print('Numpy version: ' + numpy_version)
 print('Cython speedups: ' + str(speedups.available))
 
-if lgeos.geos_version >= (3, 3, 0):
-    # Remove any WKT writer defaults to pass tests for all versions of GEOS
-    WKTWriter.defaults = {}
-
 
 shapely20_deprecated = pytest.mark.filterwarnings(
     "ignore::shapely.errors.ShapelyDeprecationWarning"
diff --git a/tests/test_collection.py b/tests/test_collection.py
index e707507..09b5c11 100644
--- a/tests/test_collection.py
+++ b/tests/test_collection.py
@@ -31,7 +31,7 @@ def test_empty(geom):
     assert geom.type == "GeometryCollection"
     assert geom.type == geom.geom_type
     assert geom.is_empty
-    assert len(geom) == 0
+    assert len(geom.geoms) == 0
     assert geom.geoms == []
 
 
@@ -53,7 +53,7 @@ def test_child_with_deleted_parent():
 def test_from_geojson(geometrycollection_geojson):
     geom = shape(geometrycollection_geojson)
     assert geom.geom_type == "GeometryCollection"
-    assert len(geom) == 2
+    assert len(geom.geoms) == 2
 
     geom_types = [g.geom_type for g in geom.geoms]
     assert "Point" in geom_types
@@ -95,3 +95,22 @@ def test_geometrycollection_adapter_deprecated(geometrycollection_geojson):
     with pytest.warns(ShapelyDeprecationWarning):
         asShape(d)
 
+
+@pytest.mark.parametrize('geom', [
+    GeometryCollection(),
+    shape({"type": "GeometryCollection", "geometries": []}),
+    shape({"type": "GeometryCollection", "geometries": [
+        {"type": "Point", "coordinates": ()},
+        {"type": "LineString", "coordinates": (())}
+    ]}),
+    wkt.loads('GEOMETRYCOLLECTION EMPTY'),
+])
+def test_len_empty_deprecated(geom):
+    with pytest.warns(ShapelyDeprecationWarning, match="__len__"):
+        assert len(geom) == 0
+
+
+def test_len_deprecated(geometrycollection_geojson):
+    geom = shape(geometrycollection_geojson)
+    with pytest.warns(ShapelyDeprecationWarning, match="__len__"):
+        assert len(geom) == 2
diff --git a/tests/test_create_inconsistent_dimensionality.py b/tests/test_create_inconsistent_dimensionality.py
index 90ad40c..a449c11 100644
--- a/tests/test_create_inconsistent_dimensionality.py
+++ b/tests/test_create_inconsistent_dimensionality.py
@@ -29,8 +29,16 @@ direct_cases = [
 
 
 wkt_cases = [
-    ('LINESTRING (1 1 1, 2 2)', 'LINESTRING (1.0000000000000000 1.0000000000000000, 2.0000000000000000 2.0000000000000000)'),
-    ('POLYGON ((0 0 0, 1 0 0, 1 1, 0 1 0, 0 0 0))', 'POLYGON ((0.0000000000000000 0.0000000000000000, 1.0000000000000000 0.0000000000000000, 1.0000000000000000 1.0000000000000000, 0.0000000000000000 1.0000000000000000, 0.0000000000000000 0.0000000000000000))')
+    # preserve 3rd dimension
+    ('MULTIPOINT (1 1 1, 2 2)', 'MULTIPOINT Z (1 1 1, 2 2 0)'),
+    ('MULTIPOINT (1 1, 2 2 2)', 'MULTIPOINT Z (1 1 0, 2 2 2)'),
+    ('LINESTRING (1 1 1, 2 2)', 'LINESTRING Z (1 1 1, 2 2 0)'),
+    ('POLYGON ((0 0 0, 1 0 0, 1 1, 0 1 0, 0 0 0))',
+     'POLYGON Z ((0 0 0, 1 0 0, 1 1 0, 0 1 0, 0 0 0))'),
+    # drop 3rd dimension
+    ('LINESTRING (1 1, 2 2 2)', 'LINESTRING (1 1, 2 2)'),
+    ('POLYGON ((0 0, 1 0 1, 1 1, 0 1, 0 0))',
+     'POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))'),
 ]
 
 
diff --git a/tests/test_empty_polygons.py b/tests/test_empty_polygons.py
index 89f42e7..d740a21 100644
--- a/tests/test_empty_polygons.py
+++ b/tests/test_empty_polygons.py
@@ -19,4 +19,4 @@ def test_multipolygon_empty_polygon():
 
 def test_multipolygon_empty_among_polygon():
     """An empty polygon passed to MultiPolygon() is ignored."""
-    assert len(MultiPolygon([Point(0,0).buffer(1.0), Polygon()])) == 1
+    assert len(MultiPolygon([Point(0,0).buffer(1.0), Polygon()]).geoms) == 1
diff --git a/tests/test_iterops.py b/tests/test_iterops.py
index 5341896..3962933 100644
--- a/tests/test_iterops.py
+++ b/tests/test_iterops.py
@@ -1,13 +1,17 @@
 """Test operator iterations
 """
-from . import unittest
+from . import unittest, shapely20_deprecated
 from shapely import iterops
+from shapely.errors import ShapelyDeprecationWarning
 from shapely.geometry import Point, Polygon
 from shapely.geos import TopologicalError
 
+import pytest
+
 
 class IterOpsTestCase(unittest.TestCase):
 
+    @shapely20_deprecated
     def test_iterops(self):
 
         coords = ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0))
@@ -43,7 +47,7 @@ class IterOpsTestCase(unittest.TestCase):
                  for p in points))),
             [[(0.5, 0.5)]])
 
-
+    @shapely20_deprecated
     def test_err(self):
         # bowtie polygon.
         coords = ((0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0), (0.0, 0.0))
@@ -55,6 +59,7 @@ class IterOpsTestCase(unittest.TestCase):
             all([isinstance(x, Polygon)
                  for x in iterops.intersects(polygon, points, True)]))
 
+    @shapely20_deprecated
     def test_topological_error(self):
         p1 = [(339, 346), (459, 346), (399, 311), (340, 277), (399, 173),
               (280, 242), (339, 415), (280, 381), (460, 207), (339, 346)]
@@ -67,3 +72,12 @@ class IterOpsTestCase(unittest.TestCase):
 
         with self.assertRaises(TopologicalError):
             list(iterops.within(polygon1, [polygon2]))
+
+    def test_iterops_deprecated(self):
+
+        coords = ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0))
+        polygon = Polygon(coords)
+        points = [Point(0.5, 0.5), Point(2.0, 2.0)]
+
+        with pytest.warns(ShapelyDeprecationWarning, match="'contains' function"):
+            list(iterops.contains(polygon, points, True))
diff --git a/tests/test_linear_referencing.py b/tests/test_linear_referencing.py
index 7d423a3..3f670e1 100644
--- a/tests/test_linear_referencing.py
+++ b/tests/test_linear_referencing.py
@@ -1,5 +1,4 @@
 from . import unittest
-from shapely.geos import geos_version
 from shapely.geometry import Point, LineString, MultiLineString
 
 
@@ -13,34 +12,28 @@ class LinearReferencingTestCase(unittest.TestCase):
             list(self.line1.coords), list(self.line2.coords)
         ])
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_line1_project(self):
         self.assertEqual(self.line1.project(self.point), 1.0)
         self.assertEqual(self.line1.project(self.point, normalized=True), 0.5)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_line2_project(self):
         self.assertEqual(self.line2.project(self.point), 1.0)
         self.assertAlmostEqual(
             self.line2.project(self.point, normalized=True), 0.16666666666, 8)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_multiline_project(self):
         self.assertEqual(self.multiline.project(self.point), 1.0)
         self.assertEqual(
             self.multiline.project(self.point, normalized=True), 0.125)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_not_supported_project(self):
         with self.assertRaises(TypeError):
             self.point.buffer(1.0).project(self.point)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_not_on_line_project(self):
         # Points that aren't on the line project to 0.
         self.assertEqual(self.line1.project(Point(-10, -10)), 0.0)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_line1_interpolate(self):
         self.assertTrue(self.line1.interpolate(0.5).equals(Point(0.5, 0.0)))
         self.assertTrue(self.line1.interpolate(-0.5).equals(Point(1.5, 0.0)))
@@ -49,20 +42,17 @@ class LinearReferencingTestCase(unittest.TestCase):
         self.assertTrue(
             self.line1.interpolate(-0.5, normalized=True).equals(Point(1, 0)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_line2_interpolate(self):
         self.assertTrue(self.line2.interpolate(0.5).equals(Point(3.0, 0.5)))
         self.assertTrue(
             self.line2.interpolate(0.5, normalized=True).equals(Point(3, 3)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_multiline_interpolate(self):
         self.assertTrue(self.multiline.interpolate(0.5).equals(Point(0.5, 0)))
         self.assertTrue(
             self.multiline.interpolate(0.5, normalized=True).equals(
                 Point(3.0, 2.0)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_line_ends_interpolate(self):
         # Distances greater than length of the line or less than
         # zero yield the line's ends.
diff --git a/tests/test_make_valid.py b/tests/test_make_valid.py
index 78c7fd6..2166467 100644
--- a/tests/test_make_valid.py
+++ b/tests/test_make_valid.py
@@ -8,7 +8,7 @@ from tests.conftest import requires_geos_38
 def test_make_valid_invalid_input():
     geom = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
     valid = make_valid(geom)
-    assert len(valid) == 2
+    assert len(valid.geoms) == 2
     assert all(geom.type == 'Polygon' for geom in valid.geoms)
 
 
diff --git a/tests/test_multipoint.py b/tests/test_multipoint.py
index b9449e9..afbffaf 100644
--- a/tests/test_multipoint.py
+++ b/tests/test_multipoint.py
@@ -140,3 +140,9 @@ def test_getitem_deprecated():
     geom = MultiPoint([[5.0, 6.0], [7.0, 8.0]])
     with pytest.warns(ShapelyDeprecationWarning, match="__getitem__"):
         part = geom[0]
+
+
+def test_len_deprecated():
+    geom = MultiPoint([[5.0, 6.0], [7.0, 8.0]])
+    with pytest.warns(ShapelyDeprecationWarning, match="__len__"):
+        assert len(geom) == 2
diff --git a/tests/test_operations.py b/tests/test_operations.py
index a02dfb4..23dfd08 100644
--- a/tests/test_operations.py
+++ b/tests/test_operations.py
@@ -3,7 +3,7 @@ import pytest
 from shapely.geometry import Point, LineString, Polygon, MultiPoint, \
                              GeometryCollection
 from shapely.wkt import loads
-from shapely.geos import TopologicalError, geos_version
+from shapely.geos import TopologicalError
 
 class OperationsTestCase(unittest.TestCase):
 
@@ -73,15 +73,13 @@ class OperationsTestCase(unittest.TestCase):
         with pytest.raises(TopologicalError):
             invalid_polygon.relate(invalid_polygon)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_hausdorff_distance(self):
         point = Point(1, 1)
         line = LineString([(2, 0), (2, 4), (3, 4)])
 
         distance = point.hausdorff_distance(line)
         self.assertEqual(distance, point.distance(Point(3, 4)))
-    
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
+
     def test_interpolate(self):
         # successful interpolation
         test_line = LineString(((1,1),(1,2)))
@@ -94,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
diff --git a/tests/test_parallel_offset.py b/tests/test_parallel_offset.py
index ae7a2f1..7dbbe55 100644
--- a/tests/test_parallel_offset.py
+++ b/tests/test_parallel_offset.py
@@ -1,10 +1,9 @@
 from . import unittest
-from shapely.geos import geos_version
 from shapely.geometry import LineString, LinearRing
 from shapely.wkt import loads
 
 class OperationsTestCase(unittest.TestCase):
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
+
     def test_parallel_offset_linestring(self):
         line1 = LineString([(0, 0), (10, 0)])
         left = line1.parallel_offset(5, 'left')
@@ -26,7 +25,6 @@ class OperationsTestCase(unittest.TestCase):
                          resolution=1),
                          LineString([(0, 2), (7, 2), (7, -5)]))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_parallel_offset_linear_ring(self):
         lr1 = LinearRing([(0, 0), (5, 0), (5, 5), (0, 5), (0, 0)])
         self.assertEqual(lr1.parallel_offset(2, 'left', resolution=1),
diff --git a/tests/test_polygonize.py b/tests/test_polygonize.py
index 219504b..a67ef8d 100755
--- a/tests/test_polygonize.py
+++ b/tests/test_polygonize.py
@@ -1,5 +1,4 @@
 from . import unittest
-from shapely.geos import geos_version
 from shapely.geometry import Point, LineString, Polygon
 from shapely.geometry.base import dump_coords
 from shapely.ops import polygonize, polygonize_full
@@ -20,7 +19,6 @@ class PolygonizeTestCase(unittest.TestCase):
         result = list(polygonize(lines))
         self.assertTrue(all([isinstance(x, Polygon) for x in result]))
 
-    @unittest.skipIf(geos_version < (3, 3, 0), 'GEOS 3.3.0 required')
     def test_polygonize_full(self):
 
         lines2 = [
@@ -34,7 +32,7 @@ class PolygonizeTestCase(unittest.TestCase):
             ]
 
         result2, dangles, cuts, invalids = polygonize_full(lines2)
-        self.assertEqual(len(result2), 2)
+        self.assertEqual(len(result2.geoms), 2)
         self.assertTrue(all([isinstance(x, Polygon) for x in result2.geoms]))
         self.assertEqual(list(dangles.geoms), [])
         self.assertTrue(all([isinstance(x, LineString) for x in cuts.geoms]))
diff --git a/tests/test_prepared.py b/tests/test_prepared.py
index f0af2d0..a5d8c0b 100644
--- a/tests/test_prepared.py
+++ b/tests/test_prepared.py
@@ -1,12 +1,8 @@
 import pytest
-from shapely.geos import geos_version
 from shapely.geometry import Point, Polygon
 from shapely.prepared import PreparedGeometry, prep
 
 
-pytestmark = pytest.mark.skipif(geos_version < (3, 1, 0), reason="GEOS 3.1.0 required")
-
-
 def test_prepared_geometry():
     polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
     p = PreparedGeometry(polygon)
diff --git a/tests/test_shared_paths.py b/tests/test_shared_paths.py
index 29b1de5..b033060 100644
--- a/tests/test_shared_paths.py
+++ b/tests/test_shared_paths.py
@@ -2,10 +2,8 @@ from . import unittest
 
 from shapely.geometry import Point, LineString, Polygon, MultiLineString, \
                              GeometryCollection
-from shapely.geos import geos_version
 from shapely.ops import shared_paths
 
-@unittest.skipIf(geos_version < (3, 3, 0), 'GEOS 3.3.0 required')
 class SharedPaths(unittest.TestCase):
     def test_shared_paths_forward(self):
         g1 = LineString([(0, 0), (10, 0), (10, 5), (20, 5)])
@@ -13,10 +11,10 @@ class SharedPaths(unittest.TestCase):
         result = shared_paths(g1, g2)
         
         self.assertTrue(isinstance(result, GeometryCollection))
-        self.assertTrue(len(result) == 2)
+        self.assertTrue(len(result.geoms) == 2)
         a, b = result.geoms
         self.assertTrue(isinstance(a, MultiLineString))
-        self.assertTrue(len(a) == 1)
+        self.assertTrue(len(a.geoms) == 1)
         self.assertEqual(a.geoms[0].coords[:], [(5, 0), (10, 0)])
         self.assertTrue(b.is_empty)
 
@@ -26,10 +24,10 @@ class SharedPaths(unittest.TestCase):
         result = shared_paths(g1, g2)
         
         self.assertTrue(isinstance(result, GeometryCollection))
-        self.assertTrue(len(result) == 2)
+        self.assertTrue(len(result.geoms) == 2)
         a, b = result.geoms
         self.assertTrue(isinstance(b, MultiLineString))
-        self.assertTrue(len(b) == 1)
+        self.assertTrue(len(b.geoms) == 1)
         self.assertEqual(b.geoms[0].coords[:], [(5, 0), (10, 0)])
         self.assertTrue(a.is_empty)
     
diff --git a/tests/test_snap.py b/tests/test_snap.py
index ad1abab..188ac1d 100644
--- a/tests/test_snap.py
+++ b/tests/test_snap.py
@@ -1,10 +1,9 @@
 from . import unittest
 
 from shapely.geometry import LineString, Polygon
-from shapely.geos import geos_version
 from shapely.ops import snap
 
-@unittest.skipIf(geos_version < (3, 3, 0), 'GEOS 3.3.0 required')
+
 class Snap(unittest.TestCase):
     def test_snap(self):
         
diff --git a/tests/test_split.py b/tests/test_split.py
index 3eaa1c8..86c77e4 100644
--- a/tests/test_split.py
+++ b/tests/test_split.py
@@ -10,7 +10,7 @@ class TestSplitGeometry(unittest.TestCase):
 	def helper(self, geom, splitter, expected_chunks):
 		s = split(geom, splitter)
 		self.assertEqual(s.type, "GeometryCollection")
-		self.assertEqual(len(s), expected_chunks)
+		self.assertEqual(len(s.geoms), expected_chunks)
 		if expected_chunks > 1:
 			# split --> expected collection that when merged is again equal to original geometry
 			if s.geoms[0].type == 'LineString':
@@ -122,6 +122,11 @@ class TestSplitLine(TestSplitGeometry):
 		self.assertTrue(splitter.touches(self.ls))
 		self.helper(self.ls, splitter, 1)
 
+		# splitter boundary touches interior of line --> return 2 segments
+		splitter = LineString([(0, 1), (1, 1)])  # touches at (1, 1)
+		self.assertTrue(splitter.touches(self.ls))
+		self.helper(self.ls, splitter, 2)
+
 	def test_split_line_with_multiline(self):
 		# crosses at one point --> return 2 segments
 		splitter = MultiLineString([[(0, 1), (1, 0)], [(0, 0), (2, -2)]])
diff --git a/tests/test_substring.py b/tests/test_substring.py
index acc5c7b..5c1658d 100644
--- a/tests/test_substring.py
+++ b/tests/test_substring.py
@@ -2,7 +2,6 @@ import json
 
 from . import unittest
 from shapely.ops import substring
-from shapely.geos import geos_version
 from shapely.geometry import Point, LineString, shape
 
 
@@ -12,76 +11,101 @@ class SubstringTestCase(unittest.TestCase):
         self.point = Point(1, 1)
         self.line1 = LineString(([0, 0], [2, 0]))
         self.line2 = LineString(([3, 0], [3, 6], [4.5, 6]))
+        self.line3 = LineString(((0, i) for i in range(5)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_startpoint(self):
         self.assertTrue(substring(self.line1, -500, -600).equals(Point(0, 0)))
         self.assertTrue(substring(self.line1, -500, -500).equals(Point(0, 0)))
         self.assertTrue(substring(self.line1, -1, -1.1, True).equals(Point(0, 0)))
         self.assertTrue(substring(self.line1, -1.1, -1.1, True).equals(Point(0, 0)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_endpoint(self):
         self.assertTrue(substring(self.line1, 500, 600).equals(Point(2, 0)))
         self.assertTrue(substring(self.line1, 500, 500).equals(Point(2, 0)))
         self.assertTrue(substring(self.line1, 1, 1.1, True).equals(Point(2, 0)))
         self.assertTrue(substring(self.line1, 1.1, 1.1, True).equals(Point(2, 0)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_midpoint(self):
         self.assertTrue(substring(self.line1, 0.5, 0.5).equals(Point(0.5, 0)))
         self.assertTrue(substring(self.line1, -0.5, -0.5).equals(Point(1.5, 0)))
         self.assertTrue(substring(self.line1, 0.5, 0.5, True).equals(Point(1, 0)))
         self.assertTrue(substring(self.line1, -0.5, -0.5, True).equals(Point(1, 0)))
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
+        # Coming from opposite ends
+        self.assertTrue(substring(self.line1, 1.5, -0.5).equals(Point(1.5, 0)))
+        self.assertTrue(substring(self.line1, -0.5, 1.5).equals(Point(1.5, 0)))
+        self.assertTrue(substring(self.line1, -0.7, 0.3, True).equals(Point(0.6, 0)))
+        self.assertTrue(substring(self.line1, 0.3, -0.7, True).equals(Point(0.6, 0)))
+
     def test_return_startsubstring(self):
-        self.assertTrue(substring(self.line1, -500, 0.6).equals(LineString(([0, 0], [0.6, 0]))))
-        self.assertTrue(substring(self.line1, -1.1, 0.6, True).equals(LineString(([0, 0], [1.2, 0]))))
+        self.assertEqual(substring(self.line1, -500, 0.6).wkt, LineString(([0, 0], [0.6, 0])).wkt)
+        self.assertEqual(substring(self.line1, -1.1, 0.6, True).wkt, LineString(([0, 0], [1.2, 0])).wkt)
+
+    def test_return_startsubstring_reversed(self):
+        # not normalized
+        self.assertEqual(substring(self.line1, -1, -500).wkt, LineString(([1, 0], [0, 0])).wkt)
+        self.assertEqual(substring(self.line3, 3.5, 0).wkt, LineString(([0, 3.5], [0, 3], [0, 2], [0, 1], [0, 0])).wkt)
+        self.assertEqual(substring(self.line3, -1.5, -500).wkt, LineString(([0, 2.5], [0, 2], [0, 1], [0, 0])).wkt)
+        # normalized
+        self.assertEqual(substring(self.line1, -0.5, -1.1, True).wkt, LineString(([1.0, 0], [0, 0])).wkt)
+        self.assertEqual(substring(self.line3, 0.5, 0, True).wkt, LineString(([0, 2.0], [0, 1], [0, 0])).wkt)
+        self.assertEqual(substring(self.line3, -0.5, -1.1, True).wkt, LineString(([0, 2.0], [0, 1], [0, 0])).wkt)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_endsubstring(self):
-        self.assertTrue(substring(self.line1, 0.6, 500).equals(LineString(([0.6, 0], [2, 0]))))
-        self.assertTrue(substring(self.line1, 0.6, 1.1, True).equals(LineString(([1.2, 0], [2, 0]))))
+        self.assertEqual(substring(self.line1, 0.6, 500).wkt, LineString(([0.6, 0], [2, 0])).wkt)
+        self.assertEqual(substring(self.line1, 0.6, 1.1, True).wkt, LineString(([1.2, 0], [2, 0])).wkt)
+
+    def test_return_endsubstring_reversed(self):
+        # not normalized
+        self.assertEqual(substring(self.line1, 500, -1).wkt, LineString(([2, 0], [1, 0])).wkt)
+        self.assertEqual(substring(self.line3, 4, 2.5).wkt, LineString(([0, 4], [0, 3], [0, 2.5])).wkt)
+        self.assertEqual(substring(self.line3, 500, -1.5).wkt, LineString(([0, 4], [0, 3], [0, 2.5])).wkt)
+        # normalized
+        self.assertEqual(substring(self.line1, 1.1, -0.5, True).wkt, LineString(([2, 0], [1.0, 0])).wkt)
+        self.assertEqual(substring(self.line3, 1, 0.5, True).wkt, LineString(([0, 4], [0, 3], [0, 2.0])).wkt)
+        self.assertEqual(substring(self.line3, 1.1, -0.5, True).wkt, LineString(([0, 4], [0, 3], [0, 2.0])).wkt)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
-    def test_return_midsubstring(self):
-        self.assertTrue(substring(self.line1, 0.5, 0.6).equals(LineString(([0.5, 0], [0.6, 0]))))
-        self.assertTrue(substring(self.line1, 0.6, 0.5).equals(LineString(([0.6, 0], [0.5, 0]))))
-        self.assertTrue(substring(self.line1, -0.5, -0.6).equals(LineString(([1.5, 0], [1.4, 0]))))
-        self.assertTrue(substring(self.line1, -0.6, -0.5).equals(LineString(([1.4, 0], [1.5, 0]))))
-        self.assertTrue(substring(self.line1, 0.5, 0.6, True).equals(LineString(([1, 0], [1.2, 0]))))
-        self.assertTrue(substring(self.line1, 0.6, 0.5, True).equals(LineString(([1.2, 0], [1, 0]))))
-        self.assertTrue(substring(self.line1, -0.5, -0.6, True).equals(LineString(([1, 0], [0.8, 0]))))
-        self.assertTrue(substring(self.line1, -0.6, -0.5, True).equals(LineString(([0.8, 0], [1, 0]))))
-
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_midsubstring(self):
-        self.assertTrue(substring(self.line1, 0.5, 0.6).equals(LineString(([0.5, 0], [0.6, 0]))))
-        self.assertTrue(substring(self.line1, 0.6, 0.5).equals(LineString(([0.6, 0], [0.5, 0]))))
-        self.assertTrue(substring(self.line1, -0.5, -0.6).equals(LineString(([1.5, 0], [1.4, 0]))))
-        self.assertTrue(substring(self.line1, -0.6, -0.5).equals(LineString(([1.4, 0], [1.5, 0]))))
-        self.assertTrue(substring(self.line1, 0.5, 0.6, True).equals(LineString(([1, 0], [1.2, 0]))))
-        self.assertTrue(substring(self.line1, 0.6, 0.5, True).equals(LineString(([1.2, 0], [1, 0]))))
-        self.assertTrue(substring(self.line1, -0.5, -0.6, True).equals(LineString(([1, 0], [0.8, 0]))))
-        self.assertTrue(substring(self.line1, -0.6, -0.5, True).equals(LineString(([0.8, 0], [1, 0]))))
-
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
+        self.assertEqual(substring(self.line1, 0.5, 0.6).wkt, LineString(([0.5, 0], [0.6, 0])).wkt)
+        self.assertEqual(substring(self.line1, -0.6, -0.5).wkt, LineString(([1.4, 0], [1.5, 0])).wkt)
+        self.assertEqual(substring(self.line1, 0.5, 0.6, True).wkt, LineString(([1, 0], [1.2, 0])).wkt)
+        self.assertEqual(substring(self.line1, -0.6, -0.5, True).wkt, LineString(([0.8, 0], [1, 0])).wkt)
+
+    def test_return_midsubstring_reversed(self):
+        self.assertEqual(substring(self.line1, 0.6, 0.5).wkt, LineString(([0.6, 0], [0.5, 0])).wkt)
+        self.assertEqual(substring(self.line1, -0.5, -0.6).wkt, LineString(([1.5, 0], [1.4, 0])).wkt)
+        self.assertEqual(substring(self.line1, 0.6, 0.5, True).wkt, LineString(([1.2, 0], [1, 0])).wkt)
+        self.assertEqual(substring(self.line1, -0.5, -0.6, True).wkt, LineString(([1, 0], [0.8, 0])).wkt)
+
+        # with vertices
+        # not normalized
+        self.assertEqual(substring(self.line3, 3.5, 2.5).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)    # (+, +)
+        self.assertEqual(substring(self.line3, -0.5, -1.5).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)  # (-, -)
+        self.assertEqual(substring(self.line3, 3.5, -1.5).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)   # (+, -)
+        self.assertEqual(substring(self.line3, -0.5, 2.5).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)   # (-, +)
+
+        # normalized
+        self.assertEqual(substring(self.line3, 0.875, 0.625, True).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)
+        self.assertEqual(substring(self.line3, -0.125, -0.375, True).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)
+        self.assertEqual(substring(self.line3, 0.875, -0.375, True).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)
+        self.assertEqual(substring(self.line3, -0.125, 0.625, True).wkt, LineString(([0, 3.5], [0, 3], [0, 2.5])).wkt)
+
     def test_return_substring_with_vertices(self):
-        self.assertTrue(substring(self.line2, 1, 7).equals(LineString(([3, 1], [3, 6], [4, 6]))))
-        self.assertTrue(substring(self.line2, 0.2, 0.9, True).equals(LineString(([3, 1.5], [3, 6], [3.75, 6]))))
-        self.assertTrue(substring(self.line2, 0, 0.9, True).equals(LineString(([3, 0], [3, 6], [3.75, 6]))))
-        self.assertTrue(substring(self.line2, 0.2, 1, True).equals(LineString(([3, 1.5], [3, 6], [4.5, 6]))))
+        self.assertEqual(substring(self.line2, 1, 7).wkt, LineString(([3, 1], [3, 6], [4, 6])).wkt)
+        self.assertEqual(substring(self.line2, 0.2, 0.9, True).wkt, LineString(([3, 1.5], [3, 6], [3.75, 6])).wkt)
+        self.assertEqual(substring(self.line2, 0, 0.9, True).wkt, LineString(([3, 0], [3, 6], [3.75, 6])).wkt)
+        self.assertEqual(substring(self.line2, 0.2, 1, True).wkt, LineString(([3, 1.5], [3, 6], [4.5, 6])).wkt)
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_substring_issue682(self):
         assert list(substring(self.line2, 0.1, 0).coords) == [(3.0, 0.1), (3.0, 0.0)]
 
-    @unittest.skipIf(geos_version < (3, 2, 0), 'GEOS 3.2.0 required')
     def test_return_substring_issue848(self):
         line = shape(json.loads(data_issue_848))
         cut_line = substring(line, 0.7, 0.8, normalized=True)
         assert len(cut_line.coords) == 53
 
+    def test_raise_type_error(self):
+        self.assertRaises(TypeError, substring, Point(0, 0), 0, 0)
+
 
 data_issue_848 = '{"type": "LineString", "coordinates": [[-87.71314, 41.96793], [-87.71312, 41.96667], [-87.71311, 41.96643], [-87.7131, 41.96635], [-87.71309, 41.9663], [-87.71303, 41.96621], [-87.71298, 41.96615], [-87.71292, 41.96611], [-87.7128, 41.96607], [-87.71268, 41.96605], [-87.71255, 41.96605], [-87.7124, 41.96605], [-87.71219, 41.96605], [-87.71173, 41.96606], [-87.71108, 41.96607], [-87.71027, 41.96607], [-87.70884, 41.96609], [-87.70763, 41.96611], [-87.70645, 41.96612], [-87.70399, 41.96613], [-87.70267, 41.96614], [-87.70166, 41.96615], [-87.70075, 41.96615], [-87.69954, 41.96615], [-87.69873, 41.96616], [-87.69789, 41.96618], [-87.69675, 41.9662], [-87.69502, 41.96621], [-87.69411, 41.96621], [-87.69145, 41.96623], [-87.69026, 41.96624], [-87.68946, 41.96625], [-87.6885, 41.96625], [-87.68718, 41.96628], [-87.68545, 41.96631], [-87.68399, 41.96632], [-87.68271, 41.96635], [-87.68159, 41.96636], [-87.68034, 41.96638], [-87.67863, 41.96641], [-87.67766, 41.96642], [-87.67741, 41.96641], [-87.67722, 41.9664], [-87.67695, 41.96638], [-87.67665, 41.96632], [-87.67638, 41.96623], [-87.67613, 41.96612], [-87.67589, 41.96596], [-87.6757, 41.96579], [-87.67557, 41.96565], [-87.67544, 41.96547], [-87.67539, 41.96536], [-87.6753, 41.96519], [-87.67524, 41.96503], [-87.67523, 41.96491], [-87.67522, 41.96477], [-87.67521, 41.96457], [-87.6752, 41.96434], [-87.67519, 41.96371], [-87.67517, 41.96175], [-87.67513, 41.96077], [-87.67505, 41.95798], [-87.67501, 41.95666], [-87.67497, 41.95513], [-87.67496, 41.95452], [-87.67491, 41.95392], [-87.67487, 41.95302], [-87.67485, 41.95202], [-87.67484, 41.95101], [-87.67479, 41.94959], [-87.67476, 41.94859], [-87.67474, 41.94703], [-87.67468, 41.94596], [-87.67466, 41.94513], [-87.67463, 41.94494], [-87.67457, 41.94474], [-87.6745, 41.94455], [-87.67442, 41.94438], [-87.6743, 41.94424], [-87.67419, 41.94414], [-87.67405, 41.94404], [-87.67386, 41.94393], [-87.67367, 41.94386], [-87.67348, 41.9438], [-87.67334, 41.94376], [-87.67311, 41.94373], [-87.67289, 41.9437], [-87.67263, 41.94369], [-87.67234, 41.94369], [-87.6715, 41.9437], [-87.67088, 41.94371], [-87.66938, 41.94373], [-87.66749, 41.94377], [-87.66585, 41.94378], [-87.66508, 41.94379], [-87.66361, 41.94381], [-87.6591, 41.94391], [-87.65767, 41.94391], [-87.65608, 41.94393], [-87.6555, 41.94394], [-87.65521, 41.94394], [-87.65503, 41.94393], [-87.65488, 41.9439], [-87.6547, 41.94386], [-87.65454, 41.9438], [-87.65441, 41.94375], [-87.65425, 41.94364], [-87.6541, 41.94351], [-87.654, 41.94342], [-87.65392, 41.94331], [-87.65382, 41.94319], [-87.65375, 41.94306], [-87.65367, 41.94292], [-87.65361, 41.9428], [-87.65355, 41.94269], [-87.65351, 41.94257], [-87.65347, 41.94238], [-87.65345, 41.94218], [-87.65338, 41.93975], [-87.65337, 41.93939], [-87.65337, 41.93893], [-87.65336, 41.93865], [-87.65333, 41.93763], [-87.65331, 41.93717], [-87.65328, 41.93627], [-87.65327, 41.93603], [-87.65323, 41.93532], [-87.65322, 41.93491], [-87.6532, 41.93445], [-87.65314, 41.93312], [-87.65313, 41.93273], [-87.6531, 41.93218], [-87.65307, 41.93151], [-87.65305, 41.9309], [-87.65302, 41.9303], [-87.65299, 41.92951], [-87.65296, 41.9287], [-87.65295, 41.92842], [-87.65294, 41.92768], [-87.65292, 41.92715], [-87.65289, 41.92599], [-87.65288, 41.92537], [-87.65287, 41.92505], [-87.65282, 41.92352], [-87.65276, 41.92172], [-87.65274, 41.92113], [-87.65264, 41.91822], [-87.65264, 41.91808], [-87.65262, 41.91763], [-87.65261, 41.91718], [-87.65255, 41.91563], [-87.6525, 41.91406], [-87.65242, 41.91377], [-87.65234, 41.91362], [-87.65223, 41.91351], [-87.65208, 41.91339], [-87.65183, 41.91322], [-87.65093, 41.9126], [-87.65017, 41.91203], [-87.64985, 41.9118], [-87.64971, 41.91171], [-87.64957, 41.91164], [-87.64948, 41.9116], [-87.64939, 41.91158], [-87.6492, 41.91153], [-87.649, 41.9115], [-87.64883, 41.9115], [-87.64863, 41.9115], [-87.64792, 41.91151], [-87.64781, 41.9115], [-87.64768, 41.91146], [-87.64756, 41.91139], [-87.64745, 41.91122], [-87.6474, 41.91112], [-87.64739, 41.91101], [-87.64738, 41.91086], [-87.64736, 41.91071], [-87.64734, 41.91061], [-87.64728, 41.91051], [-87.64718, 41.91044], [-87.64709, 41.9104], [-87.64697, 41.91036], [-87.64682, 41.91034], [-87.64664, 41.91033], [-87.64646, 41.91033], [-87.6458, 41.91034], [-87.64523, 41.91034], [-87.64348, 41.91036], [-87.64255, 41.91039], [-87.641, 41.9104], [-87.64038, 41.9104], [-87.63975, 41.9104], [-87.6393, 41.91041], [-87.63814, 41.91042], [-87.63798, 41.91041], [-87.63787, 41.91039], [-87.63771, 41.91034], [-87.63757, 41.91027], [-87.63746, 41.91021], [-87.63736, 41.91011], [-87.6373, 41.90999], [-87.63727, 41.90986], [-87.63726, 41.90973], [-87.63725, 41.90951], [-87.63723, 41.90874], [-87.63718, 41.90758], [-87.63713, 41.90607], [-87.63711, 41.90543], [-87.63702, 41.90381], [-87.63702, 41.90368], [-87.63701, 41.90334], [-87.63699, 41.90322], [-87.63694, 41.90312], [-87.63688, 41.90299], [-87.63682, 41.90292], [-87.63671, 41.90279], [-87.63659, 41.90265], [-87.63653, 41.90255], [-87.63649, 41.90245], [-87.63646, 41.90235], [-87.63647, 41.90221], [-87.63647, 41.90211], [-87.6365, 41.90202], [-87.63653, 41.9019], [-87.63659, 41.90177], [-87.63666, 41.90156], [-87.63669, 41.90143], [-87.6367, 41.90131], [-87.6367, 41.90119], [-87.63664, 41.90029], [-87.63664, 41.90008], [-87.63662, 41.89975], [-87.63658, 41.89892], [-87.63657, 41.89867], [-87.63654, 41.89761], [-87.63654, 41.89738], [-87.63649, 41.89726], [-87.63641, 41.89715], [-87.63634, 41.89708], [-87.63623, 41.89699], [-87.63595, 41.89677], [-87.63583, 41.89667], [-87.63574, 41.89654], [-87.63569, 41.89645], [-87.63568, 41.89633], [-87.63565, 41.89542], [-87.63563, 41.89434], [-87.6356, 41.89327], [-87.63558, 41.89261], [-87.63554, 41.89147], [-87.63553, 41.89051], [-87.63548, 41.8903], [-87.6354, 41.89021], [-87.63533, 41.89012], [-87.63524, 41.89007], [-87.63508, 41.89001], [-87.63493, 41.88997], [-87.63475, 41.88994], [-87.63462, 41.88991], [-87.63447, 41.88989], [-87.63436, 41.88984], [-87.63425, 41.88979], [-87.63414, 41.8897], [-87.63407, 41.88962], [-87.63402, 41.88952], [-87.63399, 41.88943], [-87.63397, 41.88897], [-87.63396, 41.88707], [-87.63391, 41.88572], [-87.63389, 41.88441], [-87.63385, 41.8827], [-87.63384, 41.88144], [-87.63378, 41.88014], [-87.63374, 41.87872], [-87.63369, 41.87726], [-87.63369, 41.87706], [-87.63365, 41.87695], [-87.63359, 41.87691], [-87.63353, 41.87688], [-87.63345, 41.87686], [-87.63338, 41.87685], [-87.63263, 41.87685], [-87.63173, 41.87686], [-87.62925, 41.87689], [-87.62821, 41.87691], [-87.62757, 41.87693], [-87.6265, 41.87696], [-87.62635, 41.87696], [-87.62603, 41.87697], [-87.62605, 41.87831], [-87.6261, 41.87951], [-87.62616, 41.88203], [-87.62619, 41.88322], [-87.62622, 41.88443], [-87.62626, 41.88534], [-87.62625, 41.88552], [-87.62625, 41.88557], [-87.62627, 41.88562], [-87.6263, 41.88566], [-87.62635, 41.88569], [-87.62642, 41.88572], [-87.6265, 41.88573], [-87.62655, 41.88574], [-87.62661, 41.88574], [-87.62683, 41.88574], [-87.62784, 41.88574], [-87.62887, 41.88574], [-87.62948, 41.88574], [-87.62982, 41.88574], [-87.62992, 41.88574], [-87.63011, 41.88574], [-87.6302, 41.88574], [-87.63089, 41.88574], [-87.63204, 41.88574], [-87.63285, 41.88573], [-87.63391, 41.88572], [-87.63396, 41.88707], [-87.63397, 41.88897], [-87.63399, 41.88943], [-87.63402, 41.88952], [-87.63407, 41.88962], [-87.63414, 41.8897], [-87.63425, 41.88979], [-87.63436, 41.88984], [-87.63447, 41.88989], [-87.63462, 41.88991], [-87.63475, 41.88994], [-87.63493, 41.88997], [-87.63508, 41.89001], [-87.63524, 41.89007], [-87.63533, 41.89012], [-87.6354, 41.89021], [-87.63548, 41.8903], [-87.63553, 41.89051], [-87.63554, 41.89147], [-87.63558, 41.89261], [-87.6356, 41.89327], [-87.63563, 41.89434], [-87.63565, 41.89542], [-87.63568, 41.89633], [-87.63569, 41.89645], [-87.63574, 41.89654], [-87.63583, 41.89667], [-87.63595, 41.89677], [-87.63623, 41.89699], [-87.63634, 41.89708], [-87.63641, 41.89715], [-87.63649, 41.89726], [-87.63654, 41.89738], [-87.63654, 41.89761], [-87.63657, 41.89867], [-87.63658, 41.89892], [-87.63662, 41.89975], [-87.63664, 41.90008], [-87.63664, 41.90029], [-87.6367, 41.90119], [-87.6367, 41.90131], [-87.63669, 41.90143], [-87.63666, 41.90156], [-87.63659, 41.90177], [-87.63653, 41.9019], [-87.6365, 41.90202], [-87.63647, 41.90211], [-87.63647, 41.90221], [-87.63646, 41.90235], [-87.63649, 41.90245], [-87.63653, 41.90255], [-87.63659, 41.90265], [-87.63671, 41.90279], [-87.63682, 41.90292], [-87.63688, 41.90299], [-87.63694, 41.90312], [-87.63699, 41.90322], [-87.63701, 41.90334], [-87.63702, 41.90368], [-87.63702, 41.90381], [-87.63711, 41.90543], [-87.63713, 41.90607], [-87.63718, 41.90758], [-87.63723, 41.90874], [-87.63725, 41.90951], [-87.63726, 41.90973], [-87.63727, 41.90986], [-87.6373, 41.90999], [-87.63736, 41.91011], [-87.63746, 41.91021], [-87.63757, 41.91027], [-87.63771, 41.91034], [-87.63787, 41.91039], [-87.63798, 41.91041], [-87.63814, 41.91042], [-87.6393, 41.91041], [-87.63975, 41.9104], [-87.64038, 41.9104], [-87.641, 41.9104], [-87.64255, 41.91039], [-87.64348, 41.91036], [-87.64523, 41.91034], [-87.6458, 41.91034], [-87.64646, 41.91033], [-87.64664, 41.91033], [-87.64682, 41.91034], [-87.64697, 41.91036], [-87.64709, 41.9104], [-87.64718, 41.91044], [-87.64728, 41.91051], [-87.64734, 41.91061], [-87.64736, 41.91071], [-87.64738, 41.91086], [-87.64739, 41.91101], [-87.6474, 41.91112], [-87.64745, 41.91122], [-87.64756, 41.91139], [-87.64768, 41.91146], [-87.64781, 41.9115], [-87.64792, 41.91151], [-87.64863, 41.9115], [-87.64883, 41.9115], [-87.649, 41.9115], [-87.6492, 41.91153], [-87.64939, 41.91158], [-87.64948, 41.9116], [-87.64957, 41.91164], [-87.64971, 41.91171], [-87.64985, 41.9118], [-87.65017, 41.91203], [-87.65093, 41.9126], [-87.65183, 41.91322], [-87.65208, 41.91339], [-87.65223, 41.91351], [-87.65234, 41.91362], [-87.65242, 41.91377], [-87.6525, 41.91406], [-87.65255, 41.91563], [-87.65261, 41.91718], [-87.65262, 41.91763], [-87.65264, 41.91808], [-87.65264, 41.91822], [-87.65274, 41.92113], [-87.65276, 41.92172], [-87.65282, 41.92352], [-87.65287, 41.92505], [-87.65288, 41.92537], [-87.65289, 41.92599], [-87.65292, 41.92715], [-87.65294, 41.92768], [-87.65295, 41.92842], [-87.65296, 41.9287], [-87.65299, 41.92951], [-87.65302, 41.9303], [-87.65305, 41.9309], [-87.65307, 41.93151], [-87.6531, 41.93218], [-87.65313, 41.93273], [-87.65314, 41.93312], [-87.6532, 41.93445], [-87.65322, 41.93491], [-87.65323, 41.93532], [-87.65327, 41.93603], [-87.65328, 41.93627], [-87.65331, 41.93717], [-87.65333, 41.93763], [-87.65336, 41.93865], [-87.65337, 41.93893], [-87.65337, 41.93939], [-87.65338, 41.93975], [-87.65345, 41.94218], [-87.65347, 41.94238], [-87.65351, 41.94257], [-87.65355, 41.94269], [-87.65361, 41.9428], [-87.65367, 41.94292], [-87.65375, 41.94306], [-87.65382, 41.94319], [-87.65392, 41.94331], [-87.654, 41.94342], [-87.6541, 41.94351], [-87.65425, 41.94364], [-87.65441, 41.94375], [-87.65454, 41.9438], [-87.6547, 41.94386], [-87.65488, 41.9439], [-87.65503, 41.94393], [-87.65521, 41.94394], [-87.6555, 41.94394], [-87.65608, 41.94393], [-87.65767, 41.94391], [-87.6591, 41.94391], [-87.66361, 41.94381], [-87.66508, 41.94379], [-87.66585, 41.94378], [-87.66749, 41.94377], [-87.66938, 41.94373], [-87.67088, 41.94371], [-87.6715, 41.9437], [-87.67234, 41.94369], [-87.67263, 41.94369], [-87.67289, 41.9437], [-87.67311, 41.94373], [-87.67334, 41.94376], [-87.67348, 41.9438], [-87.67367, 41.94386], [-87.67386, 41.94393], [-87.67405, 41.94404], [-87.67419, 41.94414], [-87.6743, 41.94424], [-87.67442, 41.94438], [-87.6745, 41.94455], [-87.67457, 41.94474], [-87.67463, 41.94494], [-87.67466, 41.94513], [-87.67468, 41.94596], [-87.67474, 41.94703], [-87.67476, 41.94859], [-87.67479, 41.94959], [-87.67484, 41.95101], [-87.67485, 41.95202], [-87.67487, 41.95302], [-87.67491, 41.95392], [-87.67496, 41.95452], [-87.67497, 41.95513], [-87.67501, 41.95666], [-87.67505, 41.95798], [-87.67513, 41.96077], [-87.67517, 41.96175], [-87.67519, 41.96371], [-87.6752, 41.96434], [-87.67521, 41.96457], [-87.67522, 41.96477], [-87.67523, 41.96491], [-87.67524, 41.96503], [-87.6753, 41.96519], [-87.67539, 41.96536], [-87.67544, 41.96547], [-87.67557, 41.96565], [-87.6757, 41.96579], [-87.67589, 41.96596], [-87.67613, 41.96612], [-87.67638, 41.96623], [-87.67665, 41.96632], [-87.67695, 41.96638], [-87.67722, 41.9664], [-87.67741, 41.96641], [-87.67766, 41.96642], [-87.67863, 41.96641], [-87.68034, 41.96638], [-87.68159, 41.96636], [-87.68271, 41.96635], [-87.68399, 41.96632], [-87.68545, 41.96631], [-87.68718, 41.96628], [-87.6885, 41.96625], [-87.68946, 41.96625], [-87.69026, 41.96624], [-87.69145, 41.96623], [-87.69411, 41.96621], [-87.69502, 41.96621], [-87.69675, 41.9662], [-87.69789, 41.96618], [-87.69873, 41.96616], [-87.69954, 41.96615], [-87.70075, 41.96615], [-87.70166, 41.96615], [-87.70267, 41.96614], [-87.70399, 41.96613], [-87.70645, 41.96612], [-87.70763, 41.96611], [-87.70884, 41.96609], [-87.71027, 41.96607], [-87.71108, 41.96607], [-87.71173, 41.96606], [-87.71219, 41.96605], [-87.7124, 41.96605], [-87.71255, 41.96605], [-87.71268, 41.96605], [-87.7128, 41.96607], [-87.71292, 41.96611], [-87.71298, 41.96615], [-87.71303, 41.96621], [-87.71309, 41.9663], [-87.7131, 41.96635], [-87.71311, 41.96643], [-87.71312, 41.96667], [-87.71314, 41.96793]]}'
diff --git a/tests/test_svg.py b/tests/test_svg.py
index fcc2679..dac3e71 100644
--- a/tests/test_svg.py
+++ b/tests/test_svg.py
@@ -174,7 +174,8 @@ class SvgTestCase(unittest.TestCase):
         self.assertSVG(GeometryCollection(), '<g />')
         # Valid
         self.assertSVG(
-            Point(7, 3).union(LineString([(4, 2), (8, 4)])),
+            GeometryCollection(
+                [Point(7, 3), LineString([(4, 2), (8, 4)])]),
             '<g><circle cx="7.0" cy="3.0" r="3.0" stroke="#555555" '
             'stroke-width="1.0" fill="#66cc99" opacity="0.6" />'
             '<polyline fill="none" stroke="#66cc99" stroke-width="2.0" '
diff --git a/tests/test_union.py b/tests/test_union.py
index 4482fc2..8ab13e4 100644
--- a/tests/test_union.py
+++ b/tests/test_union.py
@@ -4,7 +4,6 @@ import random
 from itertools import islice
 from functools import partial
 from shapely.errors import ShapelyDeprecationWarning
-from shapely.geos import geos_version
 from shapely.geometry import Point, MultiPolygon
 from shapely.ops import cascaded_union, unary_union
 
@@ -54,14 +53,12 @@ class UnionTestCase(unittest.TestCase):
             list(islice(halton(7), 20, 120)),
         )
 
-    @unittest.skipIf(geos_version < (3, 3, 0), 'GEOS 3.3.0 required')
     def test_unary_union(self):
         patches = [Point(xy).buffer(0.05) for xy in self.coords]
         u = unary_union(patches)
         self.assertEqual(u.geom_type, 'MultiPolygon')
         self.assertAlmostEqual(u.area, 0.71857254056)
 
-    @unittest.skipIf(geos_version < (3, 3, 0), 'GEOS 3.3.0 required')
     def test_unary_union_multi(self):
         # Test of multipart input based on comment by @schwehr at
         # https://github.com/Toblerity/Shapely/issues/47#issuecomment-21809308
diff --git a/tests/test_voronoi_diagram.py b/tests/test_voronoi_diagram.py
index 0d3ede6..79ec2c0 100644
--- a/tests/test_voronoi_diagram.py
+++ b/tests/test_voronoi_diagram.py
@@ -22,7 +22,7 @@ def test_no_regions():
     mp = MultiPoint(points=[(0.5, 0.5)])
     regions = voronoi_diagram(mp)
 
-    assert len(regions) == 0
+    assert len(regions.geoms) == 0
 
 
 @requires_geos_35
@@ -30,7 +30,7 @@ def test_two_regions():
     mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
     regions = voronoi_diagram(mp)
 
-    assert len(regions) == 2
+    assert len(regions.geoms) == 2
 
 
 @requires_geos_35
@@ -38,7 +38,7 @@ def test_edges():
     mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
     regions = voronoi_diagram(mp, edges=True)
 
-    assert len(regions) == 1
+    assert len(regions.geoms) == 1
     assert all(r.type == 'LineString' for r in regions.geoms)
 
 
@@ -49,7 +49,7 @@ def test_smaller_envelope():
 
     regions = voronoi_diagram(mp, envelope=poly)
 
-    assert len(regions) == 2
+    assert len(regions.geoms) == 2
     assert sum(r.area for r in regions.geoms) > poly.area
 
 
@@ -63,7 +63,7 @@ def test_larger_envelope():
 
     regions = voronoi_diagram(mp, envelope=poly)
 
-    assert len(regions) == 2
+    assert len(regions.geoms) == 2
     assert sum(r.area for r in regions.geoms) == poly.area
 
 
@@ -72,7 +72,7 @@ def test_from_polygon():
     poly = load_wkt('POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))')
     regions = voronoi_diagram(poly)
 
-    assert len(regions) == 4
+    assert len(regions.geoms) == 4
 
 
 @requires_geos_35
@@ -80,7 +80,7 @@ def test_from_polygon_with_enough_tolerance():
     poly = load_wkt('POLYGON ((0 0, 0.5 0, 0.5 0.5, 0 0.5, 0 0))')
     regions = voronoi_diagram(poly, tolerance=1.0)
 
-    assert len(regions) == 2
+    assert len(regions.geoms) == 2
 
 
 @requires_geos_35
@@ -131,5 +131,5 @@ def test_from_multipoint_without_tolerace_without_floating_point_coordinates():
     """But it's fine without it."""
     mp = load_wkt('MULTIPOINT (0 0, 1 0, 1 2, 0 1)')
     regions = voronoi_diagram(mp)
-    assert len(regions) == 4
+    assert len(regions.geoms) == 4
 
diff --git a/tests/test_wkb.py b/tests/test_wkb.py
index 22de544..ddc2f46 100644
--- a/tests/test_wkb.py
+++ b/tests/test_wkb.py
@@ -1,6 +1,12 @@
+import binascii
+
+import pytest
+
+from shapely import wkt
 from shapely.wkb import dumps, loads
 from shapely.geometry import Point
-import binascii
+from shapely.geos import geos_version
+
 
 
 def bin2hex(value):
@@ -49,3 +55,20 @@ def test_loads_srid():
     # replace geometry srid with another
     result = dumps(geom, srid=27700)
     assert bin2hex(result) == "0101000020346C0000333333333333F33F3333333333330B40"
+
+
+requires_geos_39 = pytest.mark.xfail(
+    geos_version < (3, 9, 0), reason="GEOS >= 3.9.0 is required", strict=True)
+
+
+@requires_geos_39
+def test_point_empty():
+    g = wkt.loads("POINT EMPTY")
+    assert g.wkb_hex == "0101000000000000000000F87F000000000000F87F"
+
+
+@requires_geos_39
+def test_point_z_empty():
+    g = wkt.loads("POINT Z EMPTY")
+    assert g.wkb_hex == \
+        "0101000080000000000000F87F000000000000F87F000000000000F87F"
diff --git a/tests/test_wkt.py b/tests/test_wkt.py
index 945cb14..3816227 100644
--- a/tests/test_wkt.py
+++ b/tests/test_wkt.py
@@ -18,12 +18,12 @@ def pipi4():
 
 def test_wkt(pipi):
     """.wkt and wkt.dumps() both do not trim by default."""
-    assert pipi.wkt == "POINT ({0:.16f} {1:.16f})".format(pi, -pi)
+    assert pipi.wkt == "POINT ({0:.15f} {1:.15f})".format(pi, -pi)
 
 
 def test_wkt(pipi4):
     """.wkt and wkt.dumps() both do not trim by default."""
-    assert pipi4.wkt == "POINT ({0:.16f} {1:.16f})".format(pi*4, -pi*4)
+    assert pipi4.wkt == "POINT ({0:.14f} {1:.14f})".format(pi*4, -pi*4)
 
 
 def test_dumps(pipi4):