Refactor Geometry classes to subclass the C extension type (#983)

.github/workflows/tests.yml

@@ -15,61 +15,26 @@ jobs:
             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
           # 2021
           - python: "3.10"
             geos: 3.10.0
             numpy: 1.21.3
             speedups: 1
-          - python: "3.10"
-            geos: 3.10.0
-            numpy: 1.21.3
-            speedups: 0
-          - python: "3.10"
-            geos: 3.10.0
-            speedups: 0
           # dev
           - python: "3.10"
             geos: main
@@ -102,11 +67,17 @@ jobs:
           ./ci/install_geos.sh
           cd ${{ github.workspace }}
 
+      - name: Set environment variables for PyGEOS installation
+        shell: bash
+        run: |
+          echo "GEOS_CONFIG=$GEOS_INSTALL/bin/geos-config" >> $GITHUB_ENV
+
       - name: Install python dependencies
         shell: bash
         run: |
           pip install --disable-pip-version-check --upgrade pip
           pip install --upgrade wheel
+          pip install git+https://github.com/pygeos/pygeos.git -v
           if [ "$GEOS_VERSION" = "main" ]; then
             pip install --upgrade --pre Cython numpy;
           elif [ "$SPEEDUPS" == "1" ]; then

appveyor.yml

@@ -55,17 +55,27 @@ build_script:
   - python -c "import sys; print('Python ' + sys.version)"
   - python -c "import ctypes, os; print(ctypes.CDLL(os.environ['GEOS_LIBRARY_PATH']))"
 
-  - ps: 'Write-Host "Building extension" -ForegroundColor Magenta'
+  - ps: 'Write-Host "Installing dev requirements" -ForegroundColor Magenta'
   - cd %APPVEYOR_BUILD_FOLDER%
   - pip install -r requirements-dev.txt
+
+  - ps: 'Write-Host "Installing PyGEOS" -ForegroundColor Magenta'
+  - set GEOS_LIBRARY_PATH=%GEOS_INSTALL%\lib
+  - set GEOS_INCLUDE_PATH=%GEOS_INSTALL%\include
+  - pip install git+https://github.com/pygeos/pygeos.git -v
+
+  - ps: 'Write-Host "Building extension" -ForegroundColor Magenta'
+  - set GEOS_LIBRARY_PATH=%GEOS_INSTALL%\bin\geos_c.dll
   - python setup.py build_ext --inplace -I%CYTHONPF%\\include -lgeos_c -L%CYTHONPF%\\lib
 
   - ps: 'Write-Host "Running pytest" -ForegroundColor Magenta'
   - pytest
   - python -c "import shapely; print(shapely.__version__)"
   - python -c "from shapely.geos import geos_version_string; print(geos_version_string)"
-  - python -c "from shapely import speedups; assert speedups.enabled"
-  - python -c "from shapely import vectorized"
+  # this fails on python>=3.8 with pygeos due to wrong DLL detection,
+  # temporarily skipping since the speedups module will be removed later anyway
+  #- python -c "from shapely import speedups; assert speedups.enabled"
+  #- python -c "from shapely import vectorized"
 
 
 cache:

shapely/algorithms/cga.py

@@ -1,3 +1,5 @@
+import pygeos
+
 
 def signed_area(ring):
     """Return the signed area enclosed by a ring in linear time using the
@@ -8,9 +10,13 @@ def signed_area(ring):
     ys.append(ys[1])
     return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0
 
-def is_ccw_impl(name):
+def is_ccw_impl(name=None):
     """Predicate implementation"""
     def is_ccw_op(ring):
         return signed_area(ring) >= 0.0
-    return is_ccw_op
+
+    if pygeos.geos_version >= (3, 7, 0):
+        return pygeos.is_ccw
+    else:
+        return is_ccw_op
 

shapely/coords.py

@@ -25,147 +25,48 @@ class CoordinateSequence:
 
     """
 
-    # Attributes
-    # ----------
-    # _cseq : c_void_p
-    #     Ctypes pointer to GEOS coordinate sequence
-    # _ndim : int
-    #     Number of dimensions (2 or 3, generally)
-    # __p__ : object
-    #     Parent (Shapely) geometry
-    _cseq = None
-    _ndim = None
-    __p__ = None
-
-    def __init__(self, parent):
-        self.__p__ = parent
-
-    def _update(self):
-        if self.__p__.is_empty:
-            self._ndim = None
-            self._cseq = None
-        else:
-            self._ndim = self.__p__._ndim
-            self._cseq = lgeos.GEOSGeom_getCoordSeq(self.__p__._geom)
+    def __init__(self, coords):
+        self._coords = coords
 
     def __len__(self):
-        self._update()
-        if self._cseq is None:
-            return 0
-        else:
-            cs_len = c_uint(0)
-            if self._cseq:
-                lgeos.GEOSCoordSeq_getSize(self._cseq, byref(cs_len))
-            return cs_len.value
+        return self._coords.shape[0]
 
     def __iter__(self):
-        self._update()
-        dx = c_double()
-        dy = c_double()
-        dz = c_double()
-        has_z = self._ndim == 3
         for i in range(self.__len__()):
-            lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(dx))
-            lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(dy))
-            if has_z:
-                lgeos.GEOSCoordSeq_getZ(self._cseq, i, byref(dz))
-                yield (dx.value, dy.value, dz.value)
-            else:
-                yield (dx.value, dy.value)
+            yield tuple(self._coords[i].tolist())
 
     def __getitem__(self, key):
-        self._update()
-        dx = c_double()
-        dy = c_double()
-        dz = c_double()
         m = self.__len__()
-        has_z = self._ndim == 3
         if isinstance(key, int):
-            if self._cseq is None or key + m < 0 or key >= m:
+            if key + m < 0 or key >= m:
                 raise IndexError("index out of range")
             if key < 0:
                 i = m + key
             else:
                 i = key
-            lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(dx))
-            lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(dy))
-            if has_z:
-                lgeos.GEOSCoordSeq_getZ(self._cseq, i, byref(dz))
-                return (dx.value, dy.value, dz.value)
-            else:
-                return (dx.value, dy.value)
+            return tuple(self._coords[i].tolist())
         elif isinstance(key, slice):
-            if self._cseq is None:
-                return []
-            else:
-                res = []
-                start, stop, stride = key.indices(m)
-                for i in range(start, stop, stride):
-                    lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(dx))
-                    lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(dy))
-                    if has_z:
-                        lgeos.GEOSCoordSeq_getZ(self._cseq, i, byref(dz))
-                        res.append((dx.value, dy.value, dz.value))
-                    else:
-                        res.append((dx.value, dy.value))
-                return res
+            res = []
+            start, stop, stride = key.indices(m)
+            for i in range(start, stop, stride):
+                res.append(tuple(self._coords[i].tolist()))
+            return res
         else:
             raise TypeError("key must be an index or slice")
 
-    @property
-    def _ctypes(self):
-        self._update()
-        has_z = self._ndim == 3
-        n = self._ndim or 0
-        if n == 0:
-            # ignore with NumPy 1.21 __array_interface__
-            raise AttributeError("empty geometry sequence")
-        m = self.__len__()
-        array_type = c_double * (m * n)
-        data = array_type()
-        temp = c_double()
-        for i in range(m):
-            lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(temp))
-            data[n*i] = temp.value
-            lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(temp))
-            data[n*i+1] = temp.value
-            if has_z:
-                lgeos.GEOSCoordSeq_getZ(self._cseq, i, byref(temp))
-                data[n*i+2] = temp.value
-        return data
-
-    def array_interface(self):
-        """Provide the Numpy array protocol."""
-        if sys.byteorder == 'little':
-            typestr = '<f8'
-        elif sys.byteorder == 'big':
-            typestr = '>f8'
-        else:
-            raise ValueError(
-                "Unsupported byteorder: neither little nor big-endian")
-        ai = {
-            'version': 3,
-            'typestr': typestr,
-            'data': self._ctypes,
-            }
-        ai.update({'shape': (len(self), self._ndim)})
-        return ai
-
-    __array_interface__ = property(array_interface)
+    def __array__(self, dtype=None):
+        return self._coords
 
     @property
     def xy(self):
         """X and Y arrays"""
-        self._update()
         m = self.__len__()
         x = array('d')
         y = array('d')
-        temp = c_double()
         for i in range(m):
-            lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(temp))
-            x.append(temp.value)
-            lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(temp))
-            y.append(temp.value)
+            xy = self._coords[i].tolist()
+            x.append(xy[0])
+            y.append(xy[1])
         return x, y