diff --git a/shapely/constructive.py b/shapely/constructive.py index 35e808fc3..390e0dcc2 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 including the outer boundaries of the + coverage or simplification of only the inner (shared) edges. + + 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), simplifies 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..d09b50be9 100644 --- a/shapely/tests/test_constructive.py +++ b/shapely/tests/test_constructive.py @@ -982,3 +982,66 @@ 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) + # Anything other than MultiPolygon or a GeometryCollection is returned as-is + if shapely.get_type_id(geometry) not in (6, 7): + assert actual[0].equals(geometry) + + +@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12") +def test_coverage_simplify_multipolygon(): + mp = MultiPolygon( + [ + Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]), + Polygon([(2, 2), (2, 3), (3, 3), (3, 2), (2, 2)]), + ] + ) + actual = shapely.coverage_simplify(mp, 1) + assert actual.equals( + shapely.from_wkt( + "GEOMETRYCOLLECTION (POLYGON ((0 1, 1 1, 1 0, 0 1)), " + "POLYGON ((2 3, 3 3, 3 2, 2 3)))" + ) + ) + + +@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12") +def test_coverage_simplify_collection(): + mp = shapely.GeometryCollection( + [ + shapely.Polygon([(0, 0), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)]), + shapely.Polygon([(0, 10), (10, 5), (20, 10), (20, 20), (0, 20), (0, 10)]), + ] + ) + low_tolerance = shapely.coverage_simplify(mp, 1) + mid_tolerance = shapely.coverage_simplify(mp, 8) + high_tolerance = shapely.coverage_simplify(mp, 10) + + assert low_tolerance.equals(mp.normalize()) + assert mid_tolerance.equals( + shapely.from_wkt( + "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), " + "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10)))" + ) + ) + assert high_tolerance.equals( + shapely.from_wkt( + "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 20 0, 20 10)), " + "POLYGON ((20 10, 0 10, 0 20, 20 10)))" + ) + ) + + no_boundary = shapely.coverage_simplify(mp, 10, simplify_boundary=False) + assert no_boundary.equals( + shapely.from_wkt( + "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), " + "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10)))" + ) + ) diff --git a/src/ufuncs.c b/src/ufuncs.c index 9d08fa20c..c409249e2 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] = 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; + } + } else { + geom_arr[i] = GEOSGeom_clone_r(ctx, in1); + } + } + } + + + 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; }