diff --git a/GeometryOpsCore/src/other_primitives.jl b/GeometryOpsCore/src/other_primitives.jl index b0ac8c6b7..c165f525e 100644 --- a/GeometryOpsCore/src/other_primitives.jl +++ b/GeometryOpsCore/src/other_primitives.jl @@ -159,10 +159,10 @@ on geometries from other packages and specify how to rebuild them. rebuild(geom, child_geoms; kw...) = rebuild(GI.trait(geom), geom, child_geoms; kw...) function rebuild(trait::GI.AbstractTrait, geom, child_geoms; crs=GI.crs(geom), extent=nothing) T = GI.geointerface_geomtype(trait) - if GI.is3d(geom) - # The Boolean type parameters here indicate "3d-ness" and "measure" coordinate, respectively. - return T{true,false}(child_geoms; crs, extent) - else - return T{false,false}(child_geoms; crs, extent) - end + # Check the dimensionality of the first child geometry, since it may have changed + # NOTE that without this, 2D to 3D conversions will fail + hasZ = GI.is3d(first(child_geoms)) + hasM = GI.ismeasured(first(child_geoms)) + + return T{hasZ,hasM}(child_geoms; crs, extent) end diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index e0cd484bd..44897ff3b 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -70,6 +70,7 @@ include("transformations/segmentize.jl") include("transformations/simplify.jl") include("transformations/tuples.jl") include("transformations/transform.jl") +include("transformations/forcedims.jl") include("transformations/correction/geometry_correction.jl") include("transformations/correction/closed_ring.jl") include("transformations/correction/intersecting_polygons.jl") diff --git a/src/transformations/forcedims.jl b/src/transformations/forcedims.jl new file mode 100644 index 000000000..d614273b2 --- /dev/null +++ b/src/transformations/forcedims.jl @@ -0,0 +1,36 @@ +#= +# Force dimensions (xy, xyz) + +These functions force the geometry to be 2D or 3D. They work on any geometry, vector of geometries, feature collection, or table! + +They're implemented by `apply` pretty simply. +=# + +export forcexy, forcexyz + +""" + forcexy(geom) + +Force the geometry to be 2D. Works on any geometry, vector of geometries, feature collection, or table! +""" +function forcexy(geom) + return apply(GI.PointTrait(), geom) do point + (GI.x(point), GI.y(point)) + end +end + +""" + forcexyz(geom, z = 0) + +Force the geometry to be 3D. Works on any geometry, vector of geometries, feature collection, or table! + +The `z` parameter is the default z value - if a point has no z value, it will be set to this value. +If it does, then the z value will be kept. +""" +function forcexyz(geom, z = 0) + return apply(GI.PointTrait(), geom) do point + x, y = GI.x(point), GI.y(point) + z = GI.is3d(geom) ? GI.z(point) : z + (x, y, z) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d4d25637a..3fc144484 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,6 +27,7 @@ include("helpers.jl") @safetestset "Simplify" begin include("transformations/simplify.jl") end @safetestset "Segmentize" begin include("transformations/segmentize.jl") end @safetestset "Transform" begin include("transformations/transform.jl") end +@safetestset "Force Dimensions" begin include("transformations/forcedims.jl") end @safetestset "Geometry Correction" begin include("transformations/correction/geometry_correction.jl") end @safetestset "Closed Rings" begin include("transformations/correction/closed_ring.jl") end @safetestset "Intersecting Polygons" begin include("transformations/correction/intersecting_polygons.jl") end diff --git a/test/transformations/forcedims.jl b/test/transformations/forcedims.jl new file mode 100644 index 000000000..61e4f9e3e --- /dev/null +++ b/test/transformations/forcedims.jl @@ -0,0 +1,35 @@ +import GeometryOps as GO +import GeoInterface as GI +using ..TestHelpers + +using Test + +geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), +GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) + +@testset_implementations "force dimensions" begin + + # Test forcexy on 3D geometry + geom3d = GO.transform($geom) do p + (GI.x(p), GI.y(p), 1.0) + end + @test GI.is3d(geom3d) + geom2d = GO.forcexy(geom3d) + @test !GI.is3d(geom2d) + @test GO.equals(geom2d, geom) + + # Test forcexyz with default z + geom3d_default = GO.forcexyz($geom) + @test GI.is3d(geom3d_default) + points3d = collect(GO.flatten(GI.PointTrait, geom3d_default)) + @test all(p -> GI.z(p) == 0, points3d) + + # Test forcexyz with custom z + geom3d_custom = GO.forcexyz($geom, 5.0) + @test GI.is3d(geom3d_custom) + points3d_custom = collect(GO.flatten(GI.PointTrait, geom3d_custom)) + @test all(p -> GI.z(p) == 5.0, points3d_custom) + + # Test forcexyz preserves existing z values + @test GO.equals(GO.forcexyz(geom3d), geom3d) +end diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl index 40a700e65..770e39786 100644 --- a/test/transformations/transform.jl +++ b/test/transformations/transform.jl @@ -13,3 +13,17 @@ geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), f = CoordinateTransformations.Translation(3.5, 1.5) @test GO.transform(f, $geom) == translated end + +@testset_implementations "transform 2D to 3D" begin + flat_points_raw = collect(GO.flatten(GI.PointTrait, $geom)) + flat_points_transformed = map(flat_points_raw) do p + (GI.x(p), GI.y(p), hypot(GI.x(p), GI.y(p))) + end + + geom_transformed = GO.transform($geom) do p + (GI.x(p), GI.y(p), hypot(GI.x(p), GI.y(p))) + end + @test collect(GO.flatten(GI.PointTrait, geom_transformed)) == flat_points_transformed + @test GI.is3d(geom_transformed) + @test !GI.ismeasured(geom_transformed) +end