From d8e9986f90a8feac845d35331c7cbd029cf230b9 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Tue, 2 Jan 2024 09:45:37 +0100
Subject: [PATCH] Add is3d method (#195)

Co-authored-by: Maarten Pronk <8655030+evetion@users.noreply.github.com>
---
 src/geo_interface.jl       |  1 +
 src/geos_functions.jl      |  2 ++
 test/test_geo_interface.jl | 13 +++++++++++++
 3 files changed, 16 insertions(+)

diff --git a/src/geo_interface.jl b/src/geo_interface.jl
index 231ff58..244d834 100644
--- a/src/geo_interface.jl
+++ b/src/geo_interface.jl
@@ -31,6 +31,7 @@ GeoInterface.ngeom(t::AbstractGeometryTrait, geom::PreparedGeometry) =
 GeoInterface.ngeom(::AbstractPointTrait, geom::Point) = 0
 GeoInterface.ngeom(::AbstractPointTrait, geom::PreparedGeometry) = 0
 
+GI.is3d(::AbstractGeometryTrait, geom::AbstractGeometry) = hasZ(geom)
 GI.getexterior(::AbstractPolygonTrait, geom::Polygon) = exteriorRing(geom)
 GI.gethole(::AbstractPolygonTrait, geom::Polygon, n) = interiorRing(geom, n)
 
diff --git a/src/geos_functions.jl b/src/geos_functions.jl
index 14e7334..ce8233f 100644
--- a/src/geos_functions.jl
+++ b/src/geos_functions.jl
@@ -1188,6 +1188,8 @@ function hasZ(obj::Geometry, context::GEOSContext = get_context(obj))
     end
     result != 0x00
 end
+hasZ(obj::PreparedGeometry, context::GEOSContext = get_context(obj)) =
+    hasZ(obj.ownedby, context)
 
 # Call only on LINESTRING (return 2 on exception, 1 on true, 0 on false)
 function isClosed(obj::LineString, context::GEOSContext = get_context(obj))
diff --git a/test/test_geo_interface.jl b/test/test_geo_interface.jl
index dae6d2f..6d667d9 100644
--- a/test/test_geo_interface.jl
+++ b/test/test_geo_interface.jl
@@ -11,6 +11,7 @@ const LG = LibGEOS
     @test GeoInterface.ncoord(pt) == 2
     @test GeoInterface.getcoord(pt, 1) ≈ 1.0
     @test GeoInterface.testgeometry(pt)
+    @test GeoInterface.is3d(pt) == false
     @test GeoInterface.extent(pt) == Extent(X=(1.0, 1.0), Y=(2.0, 2.0))
     plot(pt)
 
@@ -22,6 +23,7 @@ const LG = LibGEOS
     @test GeoInterface.ncoord(pt) == 3
     @test GeoInterface.getcoord(pt, 3) ≈ 3.0
     @test GeoInterface.testgeometry(pt)
+    @test GeoInterface.is3d(pt) == true
     # This doesn't return the Z extent
     @test_broken GeoInterface.extent(pt) ==
                  Extent(X = (1.0, 1.0), Y = (2.0, 2.0), Z = (3.0, 3.0))
@@ -33,6 +35,7 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(pt)
     @test GeoInterface.x(pt) == 1
     @test GeoInterface.y(pt) == 2
+    @test GeoInterface.is3d(pt) == false
     @test isnan(GeoInterface.z(pt))
 
     @inferred GeoInterface.ncoord(pt)
@@ -47,6 +50,7 @@ const LG = LibGEOS
     @test GeoInterface.coordinates(pt) ≈ Float64[] atol = 1e-5
     @test GeoInterface.geomtrait(pt) == PointTrait()
     @test GeoInterface.ncoord(pt) == 0
+    @test GeoInterface.is3d(pt) == false
     @test GeoInterface.testgeometry(pt)
 
     mpt = LibGEOS.readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
@@ -58,6 +62,7 @@ const LG = LibGEOS
     @test p isa LibGEOS.Point
     @test GeoInterface.coordinates(p) == [10, 0]
     @test GeoInterface.testgeometry(mpt)
+    @test GeoInterface.is3d(mpt) == false
     plot(mpt)
 
     @inferred GeoInterface.ncoord(mpt)
@@ -74,6 +79,7 @@ const LG = LibGEOS
     @test p isa LibGEOS.Point
     @test GeoInterface.coordinates(p) == [9, 2]
     @test GeoInterface.testgeometry(ls)
+    @test GeoInterface.is3d(ls) == false
     plot(ls)
 
     @inferred GeoInterface.ncoord(ls)
@@ -91,6 +97,7 @@ const LG = LibGEOS
     @test GeoInterface.geomtrait(mls) == MultiLineStringTrait()
     @test GeoInterface.ngeom(mls) == 2
     @test GeoInterface.testgeometry(mls)
+    @test GeoInterface.is3d(mls) == false
     plot(mls)
 
     @inferred GeoInterface.ncoord(mls)
@@ -106,6 +113,7 @@ const LG = LibGEOS
     p = GeoInterface.getgeom(lr, 3)
     @test p isa LibGEOS.Point
     @test GeoInterface.coordinates(p) == [9, 2]
+    @test GeoInterface.is3d(lr) == false
     @test GeoInterface.testgeometry(lr)
     # Cannot convert LinearRingTrait to series data for plotting
     # plot(lr)
@@ -128,6 +136,7 @@ const LG = LibGEOS
     @test ls isa LibGEOS.LinearRing
     @test GeoInterface.coordinates(ls) == coords[2]
     @test GeoInterface.testgeometry(polygon)
+    @test GeoInterface.is3d(polygon) == false
     plot(polygon)
 
     @inferred GeoInterface.ncoord(polygon)
@@ -151,6 +160,7 @@ const LG = LibGEOS
     ]]]
     @test GeoInterface.geomtrait(multipolygon) == MultiPolygonTrait()
     @test GeoInterface.testgeometry(multipolygon)
+    @test GeoInterface.is3d(multipolygon) == false
     @test GeoInterface.extent(multipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0))
     plot(multipolygon)
 
@@ -162,6 +172,7 @@ const LG = LibGEOS
     pmultipolygon = LibGEOS.prepareGeom(multipolygon)
     @test GeoInterface.geomtrait(pmultipolygon) == MultiPolygonTrait()
     @test GeoInterface.testgeometry(pmultipolygon)
+    @test GeoInterface.is3d(pmultipolygon) == false
     @test GeoInterface.extent(pmultipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0))
     LibGEOS.destroyGeom(pmultipolygon)
 
@@ -235,6 +246,7 @@ const LG = LibGEOS
     end
     @test GeoInterface.geomtrait(geomcollection) == GeometryCollectionTrait()
     @test GeoInterface.testgeometry(geomcollection)
+    @test GeoInterface.is3d(geomcollection) == false
     plot(geomcollection)
 
     @inferred GeoInterface.ncoord(geomcollection)
@@ -267,6 +279,7 @@ const LG = LibGEOS
     geomcollection = LibGEOS.readgeom("GEOMETRYCOLLECTION EMPTY")
     @test GeoInterface.ngeom(geomcollection) == 0
     @test GeoInterface.geomtrait(geomcollection) == GeometryCollectionTrait()
+    @test GeoInterface.is3d(geomcollection) == false
     @test GeoInterface.testgeometry(geomcollection)
 
     @testset "Conversion" begin