diff --git a/shapely/constructive.py b/shapely/constructive.py index 35e808fc3..faeb2180a 100644 --- a/shapely/constructive.py +++ b/shapely/constructive.py @@ -34,6 +34,7 @@ "oriented_envelope", "minimum_rotated_rectangle", "minimum_bounding_circle", + "coverage_simplify", ] @@ -91,7 +92,7 @@ def buffer( join_style="round", mitre_limit=5.0, single_sided=False, - **kwargs + **kwargs, ): """ Computes the buffer of a geometry for positive and negative buffer distance. @@ -185,7 +186,7 @@ def buffer( np.intc(join_style), mitre_limit, np.bool_(single_sided), - **kwargs + **kwargs, ) @@ -251,7 +252,7 @@ def offset_curve( np.intc(quad_segs), np.intc(join_style), np.double(mitre_limit), - **kwargs + **kwargs, ) @@ -330,7 +331,7 @@ def clip_by_rect(geometry, xmin, ymin, xmax, ymax, **kwargs): np.double(ymin), np.double(xmax), np.double(ymax), - **kwargs + **kwargs, ) @@ -861,6 +862,43 @@ def simplify(geometry, tolerance, preserve_topology=True, **kwargs): return lib.simplify(geometry, tolerance, **kwargs) +@requires_geos("3.12.0") +@multithreading_enabled +def coverage_simplify(geometry, tolerance, simplify_boundary=True, **kwargs): + """Returns a simplified version of an input geometry using the coverage simplification + + Assumes that the geometry forms a polygonal coverage, i.e. that the boundaries of + the polygons are shared with neighboring polygons. Under this assumption, the + function simplifies the edges using the Visvalingam-Whyatt algorithm, while + preserving a valid coverage. In the most simplified case, polygons are reduced to + triangles. + + The function allows simplification of all edges inclduing the boundaries of the + coverage or simplification of the inner (shared) edges only. + + Invalid polygons within the collections and geometry types other than MultiPolygons + or GeometryCollections of Polygons and MultiPolygons are returned unchanged. + + If the geometry is polygonal but does not form a valid coverage due to overlaps, + it will be simplified but it may result in invalid topology. + + + Parameters + ---------- + geometry : Geometry or array_like + tolerance : float or array_like + The degree of simplification roughly equal to the square root of the area + of triangles that will be removed. + simplify_boundary : bool, optional + By default (True), simplifiies both internal edges of the coverage as well + as its boundary. If set to False, only simplifies internal edges. + **kwargs + See :ref:`NumPy ufunc docs ` for other keyword arguments. + + """ + return lib.coverage_simplify(geometry, tolerance, simplify_boundary, **kwargs) + + @multithreading_enabled def snap(geometry, reference, tolerance, **kwargs): """Snaps an input geometry to reference geometry's vertices. diff --git a/shapely/tests/test_constructive.py b/shapely/tests/test_constructive.py index 793ebf801..e50f6057d 100644 --- a/shapely/tests/test_constructive.py +++ b/shapely/tests/test_constructive.py @@ -982,3 +982,11 @@ def test_concave_hull_kwargs(): result3 = shapely.concave_hull(mp, ratio=0) result4 = shapely.concave_hull(mp, ratio=1) assert shapely.get_num_coordinates(result4) < shapely.get_num_coordinates(result3) + + +@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12") +@pytest.mark.parametrize("geometry", all_types) +def test_coverage_simplify_geom_types(geometry): + actual = shapely.coverage_simplify([geometry, geometry], 0.0) + assert actual.shape == (2,) + assert isinstance(actual[0], Geometry) diff --git a/src/ufuncs.c b/src/ufuncs.c index 9d08fa20c..63ba9652b 100644 --- a/src/ufuncs.c +++ b/src/ufuncs.c @@ -2457,7 +2457,7 @@ static void points_func(char** args, const npy_intp* dimensions, const npy_intp* GEOS_INIT_THREADS; - char *ip1 = args[0]; + char *ip1 = args[0]; npy_intp is1 = steps[0], cs1 = steps[3]; npy_intp n = dimensions[0], n_c1 = dimensions[1]; npy_intp i; @@ -3474,6 +3474,80 @@ static PyUFuncGenericFunction to_geojson_funcs[1] = {&to_geojson_func}; #endif // GEOS_SINCE_3_10_0 +#if GEOS_SINCE_3_12_0 +static char coverage_simplify_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT}; +static void coverage_simplify_func(char** args, const npy_intp* dimensions, const npy_intp* steps, + void* data) { + GEOSGeometry* in1 = NULL; + GEOSGeometry** geom_arr; + + CHECK_NO_INPLACE_OUTPUT(3); + + // allocate a temporary array to store output GEOSGeometry objects + geom_arr = malloc(sizeof(void*) * dimensions[0]); + CHECK_ALLOC(geom_arr); + + GEOS_INIT_THREADS; + + TERNARY_LOOP { + CHECK_SIGNALS_THREADS(i); + if (errstate == PGERR_PYSIGNAL) { + destroy_geom_arr(ctx, geom_arr, i - 1); + break; + } + // get the geometry: return on error + if (!get_geom(*(GeometryObject**)ip1, &in1)) { + errstate = PGERR_NOT_A_GEOMETRY; + destroy_geom_arr(ctx, geom_arr, i - 1); + break; + } + double in2 = *(double*)ip2; + npy_bool in3 = !(*(npy_bool*)ip3); + if ((in1 == NULL) || npy_isnan(in2)) { + // in case of a missing value: return NULL (None) + geom_arr[i] = NULL; + } + else if + (GEOSGeomTypeId_r(ctx, in1) != GEOS_MULTIPOLYGON && + GEOSGeomTypeId_r(ctx, in1) != GEOS_GEOMETRYCOLLECTION) { + geom_arr[i] = GEOSGeom_clone_r(ctx, in1); + } else { + int isValid = 1; + if (GEOSGeomTypeId_r(ctx, in1) == GEOS_GEOMETRYCOLLECTION) { + int numGeoms = GEOSGetNumGeometries_r(ctx, in1); + for (int j = 0; j < numGeoms; j++) { + GEOSGeometry* geom = GEOSGetGeometryN_r(ctx, in1, j); + if (GEOSGeomTypeId_r(ctx, geom) != GEOS_POLYGON && GEOSGeomTypeId_r(ctx, geom) != GEOS_MULTIPOLYGON) { + isValid = 0; + break; + } + } + } + if (!isValid) { + geom_arr[i] = GEOSGeom_clone_r(ctx, in1); + } else { + geom_arr[i] = GEOSCoverageSimplifyVW_r(ctx, in1, in2, (int)in3); + if (geom_arr[i] == NULL) { + errstate = PGERR_GEOS_EXCEPTION; + destroy_geom_arr(ctx, geom_arr, i - 1); + break; + } + } + } + } + + + GEOS_FINISH_THREADS; + + // fill the numpy array with PyObjects while holding the GIL + if (errstate == PGERR_SUCCESS) { + geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]); + } + free(geom_arr); +} +static PyUFuncGenericFunction coverage_simplify_funcs[1] = {&coverage_simplify_func}; + +#endif // GEOS_SINCE_3_12_0 /* TODO polygonizer functions TODO prepared geometry predicate functions @@ -3733,6 +3807,10 @@ int init_ufuncs(PyObject* m, PyObject* d) { DEFINE_CUSTOM(concave_hull, 3); #endif +#if GEOS_SINCE_3_12_0 + DEFINE_CUSTOM(coverage_simplify, 3); +#endif + Py_DECREF(ufunc); return 0; }