Check the dimensions of child geoms when rebuilding

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

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")

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

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

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

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