Refactor GeoInterface implementation. (#161)

* Refactor GeoInterface implementation.

* Also fix ambiguity for current GeoInterface.

* Remove latest T<:.

src/geo_interface.jl

@@ -19,55 +19,67 @@ GeoInterface.geomtrait(::MultiPolygon) = MultiPolygonTrait()
 GeoInterface.geomtrait(::GeometryCollection) = GeometryCollectionTrait()
 GeoInterface.geomtrait(geom::PreparedGeometry) = GeoInterface.geomtrait(geom.ownedby)
 
-GeoInterface.ngeom(::AbstractGeometryTrait, geom::AbstractGeometry) =
-    isEmpty(geom) ? 0 : numGeometries(geom)
-GeoInterface.ngeom(::AbstractPointTrait, geom::AbstractGeometry) = 0
+GeoInterface.ngeom(::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) =
+    isEmpty(geom) ? 0 : Int(numGeometries(geom))
+GeoInterface.ngeom(::LineStringTrait, geom::LineString) = Int(numPoints(geom))
+GeoInterface.ngeom(::LinearRingTrait, geom::LinearRing) = Int(numPoints(geom))
+GeoInterface.ngeom(::PolygonTrait, geom::Polygon) = Int(numInteriorRings(geom)) + 1
+GeoInterface.ngeom(t::AbstractGeometryTrait, geom::PreparedGeometry) =
+    GeoInterface.ngeom(t, geom.ownedby)
+GeoInterface.ngeom(::AbstractPointTrait, geom::Point) = 0
+GeoInterface.ngeom(::AbstractPointTrait, geom::PreparedGeometry) = 0
 
-function GeoInterface.getgeom(::AbstractGeometryTrait, geom::AbstractGeometry, i)
+function GeoInterface.getgeom(
+    ::AbstractGeometryCollectionTrait,
+    geom::AbstractMultiGeometry,
+    i,
+)
     getGeometry(geom, i)
 end
-
-GeoInterface.getgeom(::AbstractPointTrait, geom::AbstractGeometry, i) = nothing
-GeoInterface.ngeom(::AbstractGeometryTrait, geom::Union{LineString,LinearRing}) =
-    numPoints(geom)
-GeoInterface.ngeom(t::AbstractPointTrait, geom::Union{LineString,LinearRing}) = 0
-GeoInterface.getgeom(::AbstractGeometryTrait, geom::Union{LineString,LinearRing}, i) =
-    Point(getPoint(geom, i))
-GeoInterface.getgeom(::AbstractPointTrait, geom::Union{LineString,LinearRing}, i) = nothing
-
-GeoInterface.ngeom(::AbstractGeometryTrait, geom::Polygon) = numInteriorRings(geom) + 1
-GeoInterface.ngeom(::AbstractPointTrait, geom::Polygon) = 0
-function GeoInterface.getgeom(::AbstractGeometryTrait, geom::Polygon, i)
+GeoInterface.getgeom(::MultiPointTrait, geom::MultiPoint, i) = getGeometry(geom, i)::Point
+GeoInterface.getgeom(::MultiLineStringTrait, geom::MultiLineString, i) =
+    getGeometry(geom, i)::LineString
+GeoInterface.getgeom(::MultiPolygonTrait, geom::MultiPolygon, i) =
+    getGeometry(geom, i)::Polygon
+GeoInterface.getgeom(
+    ::Union{LineStringTrait,LinearRingTrait},
+    geom::Union{LineString,LinearRing},
+    i,
+) = Point(getPoint(geom, i))
+GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry) = nothing
+function GeoInterface.getgeom(::PolygonTrait, geom::Polygon, i::Int)
     if i == 1
         LinearRing(exteriorRing(geom))
     else
         LinearRing(interiorRing(geom, i - 1))
     end
 end
-GeoInterface.getgeom(::AbstractPointTrait, geom::Polygon, i) = nothing
-
-GeoInterface.ngeom(t::AbstractGeometryTrait, geom::PreparedGeometry) =
-    GeoInterface.ngeom(t, geom.ownedby)
-GeoInterface.ngeom(t::AbstractPointTrait, geom::PreparedGeometry) = 0
 GeoInterface.getgeom(t::AbstractGeometryTrait, geom::PreparedGeometry, i) =
     GeoInterface.getgeom(t, geom.ownedby, i)
 GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry, i) = 0
 
-GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) =
-    isEmpty(geom) ? 0 : getCoordinateDimension(geom)
-GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) =
-    getCoordinates(getCoordSeq(geom), 1)[i]
+GeoInterface.coordinates(t::AbstractPointTrait, geom::Point) = collect(getcoord(t, geom))
+GeoInterface.coordinates(t::AbstractPointTrait, geom::AbstractGeometry) = nothing
+GeoInterface.coordinates(t::AbstractGeometryTrait, geom::AbstractGeometry) =
+    [GeoInterface.coordinates(x) for x in getgeom(t, geom)]
+GeoInterface.coordinates(t::AbstractGeometryCollectionTrait, geom::AbstractMultiGeometry) =
+    [GeoInterface.coordinates(x) for x in getgeom(t, geom)]
 
+GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) =
+    isEmpty(geom) ? 0 : Int(getCoordinateDimension(geom))
 GeoInterface.ncoord(t::AbstractGeometryTrait, geom::PreparedGeometry) =
     GeoInterface.ncoord(t, geom.ownedby)
+
+GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) =
+    getCoordinates(getCoordSeq(geom), 1)[i]
 GeoInterface.getcoord(t::AbstractGeometryTrait, geom::PreparedGeometry, i) =
     GeoInterface.getcoord(t, geom.ownedby, i)
 
 # FIXME this doesn't work for 3d geoms, Z is missing
 function GeoInterface.extent(::AbstractGeometryTrait, geom::AbstractGeometry)
     # minx, miny, maxx, maxy = getExtent(geom)
     env = envelope(geom)
-    return Extent(X = (getXMin(env), getXMax(env)), Y = (getYMin(env), getYMax(env)))
+    return Extent(; X = (getXMin(env), getXMax(env)), Y = (getYMin(env), getYMax(env)))
 end
 
 GI.convert(::Type{Point}, ::PointTrait, geom::Point; context = nothing) = geom
@@ -180,7 +192,6 @@ function _geom_to_coord_seq(geom, context)
     return seq
 end
 
-
 GeoInterface.distance(
     ::AbstractGeometryTrait,
     ::AbstractGeometryTrait,
@@ -268,7 +279,6 @@ GeoInterface.union(
 
 GeoInterfaceRecipes.@enable_geo_plots AbstractGeometry
 
-
 # -----
 # LibGeos operations for any GeoInterface.jl compatible geometries
 # -----

src/geos_types.jl

@@ -1,5 +1,13 @@
+"""All Geometries in LibGEOS are an AbstractGeometry."""
 abstract type AbstractGeometry end
 
+"""
+All MultiGeometries in LibGEOS are an AbstractMultiGeometry.
+
+Used to specialize on methods that only work on MultiGeometries, like `ngeom`.
+"""
+abstract type AbstractMultiGeometry <: AbstractGeometry end
+
 function Base.show(io::IO, geo::AbstractGeometry)
     compact = get(io, :compact, false)
     if compact
@@ -41,7 +49,7 @@ mutable struct Point <: AbstractGeometry
         Point(createPoint(x, y, z, context), context)
 end
 
-mutable struct MultiPoint <: AbstractGeometry
+mutable struct MultiPoint <: AbstractMultiGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a multipoint from a pointer - only makes sense if it is a pointer to a multipoint
@@ -111,7 +119,7 @@ mutable struct LineString <: AbstractGeometry
     end
 end
 
-mutable struct MultiLineString <: AbstractGeometry
+mutable struct MultiLineString <: AbstractMultiGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a multiline string from a multilinestring or a linestring pointer, else error
@@ -185,7 +193,6 @@ mutable struct LinearRing <: AbstractGeometry
         finalizer(destroyGeom, ring)
         ring
     end
-
 end
 
 mutable struct Polygon <: AbstractGeometry
@@ -227,7 +234,7 @@ mutable struct Polygon <: AbstractGeometry
     ) = Polygon(createPolygon(exterior, holes, context), context)
 end
 
-mutable struct MultiPolygon <: AbstractGeometry
+mutable struct MultiPolygon <: AbstractMultiGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create multipolygon using a multipolygon or polygon pointer, else error
@@ -274,7 +281,7 @@ mutable struct MultiPolygon <: AbstractGeometry
     )
 end
 
-mutable struct GeometryCollection <: AbstractGeometry
+mutable struct GeometryCollection <: AbstractMultiGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a geometric collection from a pointer to a geometric collection, else error
@@ -423,7 +430,7 @@ function compare(
         ng1 = ngeom(geo1)
         ng2 = ngeom(geo2)
         ng1 == ng2 || return false
-        for i = 1:ng1
+        for i in 1:ng1
             compare(cmp, getgeom(geo1, i), getgeom(geo2, i), ctx) || return false
         end
     end
@@ -443,7 +450,7 @@ function compare_coord_seqs(cmp, geo1, geo2, ctx)
     np1 == np2 || return false
     coords1 = Vector{Float64}(undef, ncoords1)
     coords2 = Vector{Float64}(undef, ncoords1)
-    for i = 1:np1
+    for i in 1:np1
         coordinates!(coords1, geo1, i, ctx)
         coordinates!(coords2, geo2, i, ctx)
         cmp(coords1, coords2) || return false
@@ -468,7 +475,7 @@ function compare_coord_seqs(cmp::IsApprox, geo1, geo2, ctx)
     s1 = 0.0
     s2 = 0.0
     s12 = 0.0
-    for i = 1:np1
+    for i in 1:np1
         coordinates!(coords1, geo1, i, ctx)
         coordinates!(coords2, geo2, i, ctx)
         if ncoords1 == 2
@@ -501,7 +508,7 @@ function Base.hash(geo::AbstractGeometry, h::UInt)::UInt
     if has_coord_seq(geo)
         return hash_coord_seq(geo, h)
     else
-        for i = 1:ngeom(geo)
+        for i in 1:ngeom(geo)
             h = hash(getgeom(geo, i), h)
         end
     end
@@ -514,14 +521,13 @@ function hash_coord_seq(geo::HasCoordSeq, h::UInt)::UInt
     end
     buf = Vector{Float64}(undef, nc)
     ctx = get_context(geo)
-    for i = 1:npoints(geo)
+    for i in 1:npoints(geo)
         coordinates!(buf, geo, i, ctx)
         h = hash(buf, h)
     end
     return h
 end
 
-
 # teach ccall how to get the pointer to pass to libgeos
 # this way the Julia compiler will track the lifetimes for us
 Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractGeometry) = x.ptr

test/test_geo_interface.jl

@@ -11,7 +11,7 @@ const LG = LibGEOS
     @test GeoInterface.ncoord(pt) == 2
     @test GeoInterface.getcoord(pt, 1) ≈ 1.0
     @test GeoInterface.testgeometry(pt)
-    @test GeoInterface.extent(pt) == Extent(X = (1.0, 1.0), Y = (2.0, 2.0))
+    @test GeoInterface.extent(pt) == Extent(X=(1.0, 1.0), Y=(2.0, 2.0))
     plot(pt)
 
     pt = LibGEOS.Point(1.0, 2.0, 3.0)
@@ -32,6 +32,11 @@ const LG = LibGEOS
     @test GeoInterface.geomtrait(pt) == PointTrait()
     @test GeoInterface.testgeometry(pt)
 
+    @inferred GeoInterface.ncoord(pt)
+    @inferred GeoInterface.ngeom(pt)
+    @inferred GeoInterface.getgeom(pt)
+    @inferred GeoInterface.coordinates(pt)
+
     pt = LibGEOS.readgeom("POINT EMPTY")
     @test GeoInterface.coordinates(pt) ≈ Float64[] atol = 1e-5
     @test GeoInterface.geomtrait(pt) == PointTrait()
@@ -49,6 +54,11 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(mpt)
     plot(mpt)
 
+    @inferred GeoInterface.ncoord(mpt)
+    @inferred GeoInterface.ngeom(mpt)
+    @inferred GeoInterface.getgeom(mpt)
+    @inferred GeoInterface.coordinates(mpt)
+
     coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2]]
     ls = LibGEOS.LineString(coords)
     @test GeoInterface.coordinates(ls) == coords
@@ -60,6 +70,11 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(ls)
     plot(ls)
 
+    @inferred GeoInterface.ncoord(ls)
+    @inferred GeoInterface.ngeom(ls)
+    @inferred GeoInterface.getgeom(ls)
+    @inferred GeoInterface.coordinates(ls)
+
     ls = LibGEOS.readgeom("LINESTRING EMPTY")
     @test GeoInterface.coordinates(ls) == []
     @test GeoInterface.geomtrait(ls) == LineStringTrait()
@@ -72,6 +87,11 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(mls)
     plot(mls)
 
+    @inferred GeoInterface.ncoord(mls)
+    @inferred GeoInterface.ngeom(mls)
+    @inferred GeoInterface.getgeom(mls)
+    @inferred GeoInterface.coordinates(mls)
+
     coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2], [8, 1]]
     lr = LibGEOS.LinearRing(coords)
     @test GeoInterface.coordinates(lr) == coords
@@ -84,6 +104,11 @@ const LG = LibGEOS
     # Cannot convert LinearRingTrait to series data for plotting
     # plot(lr)
 
+    @inferred GeoInterface.ncoord(lr)
+    @inferred GeoInterface.ngeom(lr)
+    @inferred GeoInterface.getgeom(lr)
+    @inferred GeoInterface.coordinates(lr)
+
     coords = Vector{Vector{Float64}}[
         Vector{Float64}[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]],
         Vector{Float64}[[1, 8], [2, 8], [2, 9], [1, 9], [1, 8]],
@@ -99,6 +124,11 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(polygon)
     plot(polygon)
 
+    @inferred GeoInterface.ncoord(polygon)
+    @inferred GeoInterface.ngeom(polygon)
+    @inferred GeoInterface.getgeom(polygon)
+    @inferred GeoInterface.coordinates(polygon)
+
     polygon = LibGEOS.readgeom("POLYGON EMPTY")
     @test GeoInterface.coordinates(polygon) == [[]]
     @test GeoInterface.geomtrait(polygon) == PolygonTrait()
@@ -115,13 +145,18 @@ const LG = LibGEOS
     ]]]
     @test GeoInterface.geomtrait(multipolygon) == MultiPolygonTrait()
     @test GeoInterface.testgeometry(multipolygon)
-    @test GeoInterface.extent(multipolygon) == Extent(X = (0.0, 10.0), Y = (0.0, 10.0))
+    @test GeoInterface.extent(multipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0))
     plot(multipolygon)
 
+    @inferred GeoInterface.ncoord(multipolygon)
+    @inferred GeoInterface.ngeom(multipolygon)
+    @inferred GeoInterface.getgeom(multipolygon)
+    @inferred GeoInterface.coordinates(multipolygon)
+
     pmultipolygon = LibGEOS.prepareGeom(multipolygon)
     @test GeoInterface.geomtrait(pmultipolygon) == MultiPolygonTrait()
     @test GeoInterface.testgeometry(pmultipolygon)
-    @test GeoInterface.extent(pmultipolygon) == Extent(X = (0.0, 10.0), Y = (0.0, 10.0))
+    @test GeoInterface.extent(pmultipolygon) == Extent(X=(0.0, 10.0), Y=(0.0, 10.0))
     LibGEOS.destroyGeom(pmultipolygon)
 
     geomcollection = LibGEOS.readgeom(
@@ -196,6 +231,11 @@ const LG = LibGEOS
     @test GeoInterface.testgeometry(geomcollection)
     plot(geomcollection)
 
+    @inferred GeoInterface.ncoord(geomcollection)
+    @inferred GeoInterface.ngeom(geomcollection)
+    @inferred GeoInterface.getgeom(geomcollection)
+    # @inferred GeoInterface.coordinates(geomcollection)  # can't be inferred
+
     geomcollection = LibGEOS.readgeom(
         "GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))",
     )