迁移检测线串与预准备几何体相交的Cython代码至Shapely2.0遇段错误

我编写了自定义Cython代码用于检测线串与预准备几何体的相交情况，从Shapely 1.8迁移至2.0后，以下代码行触发段错误：

result[i] = <np.uint8_t> GEOSPreparedIntersects_r(geos_handle, geom1, geom2)

完整Cython代码 #

#!python
#cython: language_level=3

#include <stdio.h>
#include <stddef.h>
#include <stdint.h>

cimport cython
from libc.stdint cimport uintptr_t
import shapely.prepared
import numpy as np
cimport numpy as np

__all__ = ['two_points_intersect_geom', "DTYPE"]

np.import_array()
DTYPE = np.float32

ctypedef np.float32_t DTYPE_t

cdef extern from "geos_c.h":
    ctypedef void *GEOSContextHandle_t
    ctypedef struct GEOSGeometry
    ctypedef struct GEOSCoordSequence
    ctypedef struct GEOSPreparedGeometry

    GEOSCoordSequence *GEOSCoordSeq_create_r(GEOSContextHandle_t, unsigned int, unsigned int) nogil
    int GEOSCoordSeq_getSize_r(GEOSContextHandle_t, GEOSCoordSequence *, unsigned int *) nogil
    int GEOSCoordSeq_setX_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
    int GEOSCoordSeq_setY_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
    int GEOSCoordSeq_setZ_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
    GEOSGeometry *GEOSGeom_createLineString_r(GEOSContextHandle_t, GEOSCoordSequence *) nogil

    char GEOSPreparedIntersects_r(GEOSContextHandle_t, const GEOSPreparedGeometry *, const GEOSGeometry *) nogil
    char GEOSIntersects_r(GEOSContextHandle_t, const GEOSGeometry *, const GEOSGeometry *) nogil

cdef GEOSContextHandle_t get_geos_context_handle():
    # Note: This requires that lgeos is defined, so needs to be imported as:
    from shapely.geos import lgeos
    cdef uintptr_t handle = lgeos.geos_handle
    return <GEOSContextHandle_t> handle

cdef GEOSPreparedGeometry *geos_from_prepared(shapely_geom) except *:
    """Get the Prepared GEOS geometry pointer from the given shapely geometry."""
    cdef uintptr_t geos_geom = shapely_geom._geom
    return <GEOSPreparedGeometry *> geos_geom

@cython.boundscheck(False)
@cython.wraparound(False)
def two_points_intersect_geom(np.ndarray[DTYPE_t, ndim=3] latlon, geometry):
    """
    Example:
        import numpy as np
        import cartopy.feature as cfeature
        from shapely.ops import unary_union
        from shapely.prepared import prep
        land = prep(unary_union(list(cfeature.NaturalEarthFeature('physical', 'land', '50m').geometries())))
        latlon = np.array([
            [[0, 0], [0, 10]],
            [[0, 0], [0, -10]],
        ], dtype=float)
        two_points_intersect_geom(latlon, land)
    """
    cdef GEOSCoordSequence *coord_sequence
    cdef GEOSPreparedGeometry *geom1
    cdef GEOSGeometry *geom2
    cdef double lat, lon
    cdef int n_point_pairs = len(latlon)
    cdef int seqSize = 2
    cdef int seqDim = 2
    cdef int i, j
    cdef np.ndarray[np.uint8_t, ndim=1, cast=True] result = np.empty(n_point_pairs, dtype=np.uint8)

    if not isinstance(geometry, shapely.prepared.PreparedGeometry):
        geometry = shapely.prepared.prep(geometry)

    geos_handle = get_geos_context_handle()
    geom1 = geos_from_prepared(geometry)
    for i in range(n_point_pairs):
        coord_sequence = GEOSCoordSeq_create_r(geos_handle, seqSize, seqDim)
        for j in range(2):
            lat = latlon[i][j][0]
            lon = latlon[i][j][1]
            d = GEOSCoordSeq_setX_r(geos_handle, coord_sequence, j, lat)
            d = GEOSCoordSeq_setY_r(geos_handle, coord_sequence, j, lon)
        geom2 = GEOSGeom_createLineString_r(geos_handle, coord_sequence)
        result[i] = <np.uint8_t> GEOSPreparedIntersects_r(geos_handle, geom1, geom2)
    return result.view(dtype=np.bool_)

测试代码 #

geom = Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
land = prep(geom)

latlon = np.array(
    [
        [[-0.5, 0.5], [0.5, 0.5]],
        [[10, 10], [20, 20]],
    ],
    dtype=DTYPE_CYTHON,
)
ret = two_points_intersect_geom(latlon, land)
self.assertListEqual([True, False], list(ret))

排查情况 #

已确认段错误并非来自geos_handle、geom2或结果赋值操作，问题根源在geom1，怀疑获取预准备几何体GEOS指针的函数存在问题。查阅Shapely与GEOS的变更记录后未找到失效原因，且当前使用的是未弃用的_geom属性，而非已弃用的__geom__。

Shapely 2.0对PreparedGeometry的内部结构做了核心变更：_geom属性不再指向GEOSPreparedGeometry指针，而是指向原始的GEOSGeometry对象。要获取预准备几何体的指针，需要访问新增的_prepared_geom内部属性。

修改geos_from_prepared函数 #

将原函数中的_geom替换为_prepared_geom：

cdef GEOSPreparedGeometry *geos_from_prepared(shapely_geom) except *:
    """Get the Prepared GEOS geometry pointer from the given shapely geometry."""
    cdef uintptr_t geos_geom = shapely_geom._prepared_geom
    return <GEOSPreparedGeometry *> geos_geom

额外注意事项 #

• 必须确保输入的shapely_geom确实是PreparedGeometry对象，普通Geometry没有_prepared_geom属性，误访问会引发错误
• _prepared_geom是Shapely 2.0的内部属性，虽目前未被弃用，但后续版本可能发生变化，建议优先使用Shapely提供的公共API，尽量避免直接操作内部指针

内容的提问来源于stack exchange，提问作者Tom McLean