diff --git a/src/geo_interface.jl b/src/geo_interface.jl index 32214e1..c9d84ac 100644 --- a/src/geo_interface.jl +++ b/src/geo_interface.jl @@ -332,6 +332,8 @@ for f in ( :envelope, :minimumRotatedRectangle, :convexhull, + :concavehull, + :concavehull_of_polygons, :boundary, :unaryUnion, :pointOnSurface, diff --git a/src/geos_functions.jl b/src/geos_functions.jl index ce8233f..d1e2db1 100644 --- a/src/geos_functions.jl +++ b/src/geos_functions.jl @@ -677,6 +677,53 @@ function convexhull(obj::Geometry, context::GEOSContext = get_context(obj)) geomFromGEOS(result, context) end +""" + concavehull(geom, ratio; bylength = false, allow_holes = true) + +Compute the concave hull of a geometry, according to the chi-shape generated from +`geom` with the given `ratio`. + +If `bylength` is true, the chi-shape is generated from the length of the geometry, +calling `GEOSConcaveHullByLength`; otherwise it is generated from the area, calling +`GEOSConcaveHull`. + +If `allow_holes` is true, the concave hull may have holes. If false, then holes are disallowed. + +See also [`convexhull`](@ref), [`concavehull_of_polygons`](@ref). +""" +function concavehull(obj::Geometry; ratio::Union{Real, Nothing} = nothing, length_ratio::Union{Real, Nothing} = nothing, allow_holes::Bool = true, context::GEOSContext = get_context(obj)) + result = if isnothing(ratio) && isnothing(length_ratio) + throw(ArgumentError("Either `ratio` or `length_ratio` must be provided to `LibGEOS.concavehull`, none were provided.")) + elseif !isnothing(ratio) && !isnothing(length_ratio) + throw(ArgumentError("Only one of `ratio` or `length_ratio` may be provided to `LibGEOS.concavehull`; got both.")) + elseif !isnothing(ratio) && isnothing(length_ratio) + GEOSConcaveHull_r(context, obj, ratio, allow_holes) + elseif isnothing(ratio) && !isnothing(length_ratio) + GEOSConcaveHullByLength_r(context, obj, ratio, allow_holes) + end + if result == C_NULL + error("LibGEOS: Error in GEOSConcaveHull") + end + geomFromGEOS(result, context) +end + +""" + concavehull_of_polygons(obj::Geometry, ratio; tight = false, allow_holes = true) + +Compute the concave hull of a geometry, according to the chi-shape generated from +`obj` with the given `ratio`. + +If `tight` is true, the chi-shape is generated from the tightest fit of the geometry, +otherwise it is generated from the area. +""" +function concavehull_of_polygons(obj::Geometry, ratio::Real; tight::Bool = false, allow_holes::Bool = true, context::GEOSContext = get_context(obj)) + result = GEOSConcaveHullOfPolygons_r(context, obj, ratio, tight, allow_holes) + if result == C_NULL + error("LibGEOS: Error in GEOSConcaveHullOfPolygons") + end + geomFromGEOS(result, context) +end + function difference( obj1::Geometry, obj2::Geometry, diff --git a/test/test_geo_interface.jl b/test/test_geo_interface.jl index d13ec02..d7aa110 100644 --- a/test/test_geo_interface.jl +++ b/test/test_geo_interface.jl @@ -318,6 +318,7 @@ const LG = LibGEOS LG.delaunayTriangulationEdges, LG.delaunayTriangulation, LG.constrainedDelaunayTriangulation, + Base.Fix2(LG.concavehull, 0.3), # these have different signatures # LG.simplify, LG.topologyPreserveSimplify, ) @@ -358,12 +359,16 @@ const LG = LibGEOS @test geom isa MultiPoint @test GeoInterface.coordinates(geom) == coords for f in one_arg_functions - @test f(LibGEOS.MultiPoint(coords)) == f(GeoInterface.MultiPoint(coords)) + @testset let current_function = f + @test f(LibGEOS.MultiPoint(coords)) == f(GeoInterface.MultiPoint(coords)) + end end coords2 = [[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]] for f in two_arg_functions - @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) == - f(GeoInterface.LineString(coords), GeoInterface.LineString(coords)) + @testset let current_function = f + @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) == + f(GeoInterface.LineString(coords), GeoInterface.LineString(coords)) + end end coords = [[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]] @@ -371,12 +376,16 @@ const LG = LibGEOS @test geom isa LineString @test GeoInterface.coordinates(geom) == coords for f in one_arg_functions - @test f(LibGEOS.LineString(coords)) == f(GeoInterface.LineString(coords)) + @testset let current_function = f + @test f(LibGEOS.LineString(coords)) == f(GeoInterface.LineString(coords)) + end end coords2 = [[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]] for f in two_arg_functions - @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) == - f(GeoInterface.LineString(coords), GeoInterface.LineString(coords)) + @testset let current_function = f + @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) == + f(GeoInterface.LineString(coords), GeoInterface.LineString(coords)) + end end coords = [[[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]]] @@ -384,13 +393,17 @@ const LG = LibGEOS @test geom isa MultiLineString @test GeoInterface.coordinates(geom) == coords for f in one_arg_functions - @test f(LibGEOS.MultiLineString(coords)) == - f(GeoInterface.MultiLineString(coords)) + @testset let current_function = f + @test f(LibGEOS.MultiLineString(coords)) == + f(GeoInterface.MultiLineString(coords)) + end end coords2 = [[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]] for f in two_arg_functions - @test f(LibGEOS.MultiLineString(coords), LibGEOS.MultiLineString(coords2)) == - f(GeoInterface.MultiLineString(coords), LibGEOS.MultiLineString(coords2)) + @testset let current_function = f + @test f(LibGEOS.MultiLineString(coords), LibGEOS.MultiLineString(coords2)) == + f(GeoInterface.MultiLineString(coords), LibGEOS.MultiLineString(coords2)) + end end coords = [[[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]]] @@ -401,12 +414,16 @@ const LG = LibGEOS @test GeoInterface.nhole(geom) == 0 @test GeoInterface.coordinates(geom) == coords for f in one_arg_functions - @test f(LibGEOS.Polygon(coords)) == f(GeoInterface.Polygon(coords)) + @testset let current_function = f + @test f(LibGEOS.Polygon(coords)) == f(GeoInterface.Polygon(coords)) + end end coords2 = [[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]] for f in two_arg_functions - @test f(LibGEOS.Polygon(coords), LibGEOS.Polygon(coords2)) == - f(GeoInterface.Polygon(coords), LibGEOS.Polygon(coords2)) + @testset let current_function = f + @test f(LibGEOS.Polygon(coords), LibGEOS.Polygon(coords2)) == + f(GeoInterface.Polygon(coords), LibGEOS.Polygon(coords2)) + end end pgeom = LibGEOS.prepareGeom(geom) @@ -418,12 +435,16 @@ const LG = LibGEOS @test geom isa MultiPolygon @test GeoInterface.coordinates(geom) == coords for f in one_arg_functions - @test f(LibGEOS.MultiPolygon(coords)) == f(GeoInterface.MultiPolygon(coords)) + @testset let current_function = f + @test f(LibGEOS.MultiPolygon(coords)) == f(GeoInterface.MultiPolygon(coords)) + end end coords2 = [[[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]]] for f in two_arg_functions - @test f(LibGEOS.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) == - f(GeoInterface.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) + @testset let current_function = f + @test f(LibGEOS.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) == + f(GeoInterface.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) + end end end