Implement GeoInterface isempty, x, and y. (#191)

* Implement GeoInterface isempty.

* Implement GI X and Y.

* Ignore GeoInterfaceRecipes ambiguity.

* Add tests. Speed up some Point methods.

* Fix test.

---------

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>

src/geo_interface.jl

@@ -19,6 +19,8 @@ GeoInterface.geomtrait(::MultiPolygon) = MultiPolygonTrait()
 GeoInterface.geomtrait(::GeometryCollection) = GeometryCollectionTrait()
 GeoInterface.geomtrait(geom::PreparedGeometry) = GeoInterface.geomtrait(geom.ownedby)
 
+GeoInterface.isempty(::AbstractGeometryTrait, geom::AbstractGeometry) = isEmpty(geom)
+
 GeoInterface.ngeom(::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) =
     isEmpty(geom) ? 0 : Int(numGeometries(geom))
 GeoInterface.ngeom(::LineStringTrait, geom::LineString) = Int(numPoints(geom))
@@ -45,7 +47,7 @@ GeoInterface.getgeom(
     ::Union{LineStringTrait,LinearRingTrait},
     geom::Union{LineString,LinearRing},
     i,
-) = Point(getPoint(geom, i))
+) = getPoint(geom, i)
 GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry) = nothing
 function GeoInterface.getgeom(::PolygonTrait, geom::Polygon, i::Int)
     if i == 1
@@ -70,11 +72,15 @@ GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) =
 GeoInterface.ncoord(t::AbstractGeometryTrait, geom::PreparedGeometry) =
     GeoInterface.ncoord(t, geom.ownedby)
 
-GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) =
+GeoInterface.getcoord(::AbstractPointTrait, geom::Point, i) =
     getCoordinates(getCoordSeq(geom), 1)[i]
-GeoInterface.getcoord(t::AbstractGeometryTrait, geom::PreparedGeometry, i) =
+GeoInterface.getcoord(t::AbstractPointTrait, geom::PreparedGeometry, i) =
     GeoInterface.getcoord(t, geom.ownedby, i)
 
+GeoInterface.x(::AbstractPointTrait, point::AbstractGeometry) = getGeomX(point)
+GeoInterface.y(::AbstractPointTrait, point::AbstractGeometry) = getGeomY(point)
+GeoInterface.z(::AbstractPointTrait, point::AbstractGeometry) = getGeomZ(point)
+
 # FIXME this doesn't work for 3d geoms, Z is missing
 function GeoInterface.extent(::AbstractGeometryTrait, geom::AbstractGeometry)
     # minx, miny, maxx, maxy = getExtent(geom)

src/geos_functions.jl

@@ -347,14 +347,19 @@ function getCoordinates!(out::Vector{Float64}, ptr::GEOSCoordSeq, i::Integer, co
     end
     ndim = getDimensions(ptr, context)
     @assert length(out) == ndim
-    GC.@preserve out begin
-        start = pointer(out)
-        floatsize = sizeof(Float64)
-        GEOSCoordSeq_getX_r(context, ptr, i - 1, start)
-        GEOSCoordSeq_getY_r(context, ptr, i - 1, start + floatsize)
-        if ndim == 3
-            GEOSCoordSeq_getZ_r(context, ptr, i - 1, start + 2 * floatsize)
-        end
+    start = pointer(out)
+    floatsize = sizeof(Float64)
+    if ndim == 3
+        GEOSCoordSeq_getXYZ_r(
+            context,
+            ptr,
+            i - 1,
+            start,
+            start + floatsize,
+            start + 2 * floatsize,
+        )
+    else
+        GEOSCoordSeq_getXY_r(context, ptr, i - 1, start, start + floatsize)
     end
     out
 end
@@ -365,13 +370,12 @@ function getCoordinates(
     context::GEOSContext = get_global_context(),
 )
     ndim = getDimensions(ptr, context)
-    out = Array{Float64}(undef, ndim)
+    out = Vector{Float64}(undef, ndim)
     getCoordinates!(out, ptr, i, context)
     out
 end
 
 function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
-    ndim = getDimensions(ptr, context)
     ncoords = getSize(ptr, context)
     coordseq = Vector{Float64}[]
     sizehint!(coordseq, ncoords)
@@ -1369,6 +1373,15 @@ function getGeomY(obj::Point, context::GEOSContext = get_context(obj))
     out[]
 end
 
+function getGeomZ(obj::Point, context::GEOSContext = get_context(obj))
+    out = Ref{Float64}()
+    result = GEOSGeomGetZ_r(context, obj, out)
+    if result == -1
+        error("LibGEOS: Error in GEOSGeomGetY")
+    end
+    out[]
+end
+
 # Return NULL on exception, Geometry must be a Polygon.
 # Returned object is a pointer to internal storage: it must NOT be destroyed directly.
 function interiorRing(obj::Polygon, n::Integer, context::GEOSContext = get_context(obj))

test/runtests.jl

@@ -1,3 +1,5 @@
+
+using GeoInterface, GeoInterfaceRecipes, Extents
 using GeoInterface, Extents
 using Test, LibGEOS, RecipesBase
 import Aqua
@@ -25,4 +27,5 @@ end
     include("test_invalid_geometry.jl")
     include("test_strtree.jl")
     include("test_misc.jl")
+
 end

test/test_geo_interface.jl

@@ -31,11 +31,17 @@ const LG = LibGEOS
     @test GeoInterface.coordinates(pt) ≈ [1, 2] atol = 1e-5
     @test GeoInterface.geomtrait(pt) == PointTrait()
     @test GeoInterface.testgeometry(pt)
+    @test GeoInterface.x(pt) == 1
+    @test GeoInterface.y(pt) == 2
+    @test isnan(GeoInterface.z(pt))
 
     @inferred GeoInterface.ncoord(pt)
     @inferred GeoInterface.ngeom(pt)
     @inferred GeoInterface.getgeom(pt)
     @inferred GeoInterface.coordinates(pt)
+    @inferred GeoInterface.x(pt)
+    @inferred GeoInterface.y(pt)
+    @inferred GeoInterface.z(pt)
 
     pt = LibGEOS.readgeom("POINT EMPTY")
     @test GeoInterface.coordinates(pt) ≈ Float64[] atol = 1e-5

test/test_geos_operations.jl

@@ -45,12 +45,14 @@ end
     expected = readgeom("LINESTRING (130 240, 650 240)")
     output = convexhull(input)
     @test !isEmpty(output)
+    @test !GeoInterface.isempty(output)
     @test writegeom(output) == writegeom(expected)
 
     # LibGEOS.delaunayTriangulationTest
     g1 = readgeom("POLYGON EMPTY")
     g2 = delaunayTriangulationEdges(g1)
     @test isEmpty(g1)
+    @test GeoInterface.isempty(g1)
     @test isEmpty(g2)
     @test GeoInterface.geomtrait(g2) == MultiLineStringTrait()