Geointerface conversions

src/LibGEOS.jl

@@ -6,6 +6,8 @@ using GeoInterfaceRecipes
 using Extents
 using CEnum
 
+const GI = GeoInterface
+
 export Point,
     MultiPoint,
     LineString,

src/geo_interface.jl

@@ -1,5 +1,14 @@
 GeoInterface.isgeometry(::Type{<:AbstractGeometry}) = true
 
+geointerface_geomtype(::PointTrait) = Point
+geointerface_geomtype(::MultiPointTrait) = MultiPoint
+geointerface_geomtype(::LineStringTrait) = LineString
+geointerface_geomtype(::MultiLineStringTrait) = MultiLineString
+geointerface_geomtype(::LinearRingTrait) = LinearRing
+geointerface_geomtype(::PolygonTrait) = Polygon
+geointerface_geomtype(::MultiPolygonTrait) = MultiPolygon
+geointerface_geomtype(::GeometryCollectionTrait) = GeometryCollection
+
 GeoInterface.geomtrait(::Point) = PointTrait()
 GeoInterface.geomtrait(::MultiPoint) = MultiPointTrait()
 GeoInterface.geomtrait(::LineString) = LineStringTrait()
@@ -44,6 +53,10 @@ GeoInterface.getgeom(t::AbstractGeometryTrait, geom::PreparedGeometry, i) =
     GeoInterface.getgeom(t, geom.ownedby, i)
 GeoInterface.getgeom(t::AbstractPointTrait, geom::PreparedGeometry, i) = 0
 
+GeoInterface.x(::AbstractPointTrait, geom::AbstractGeometry) = getX(geom.ptr, 1, get_context(geom))
+GeoInterface.y(::AbstractPointTrait, geom::AbstractGeometry) = getY(geom.ptr, 1, get_context(geom))
+GeoInterface.z(::AbstractPointTrait, geom::AbstractGeometry) = getZ(geom.ptr, 1, get_context(geom))
+
 GeoInterface.ncoord(::AbstractGeometryTrait, geom::AbstractGeometry) =
     isEmpty(geom) ? 0 : getCoordinateDimension(geom)
 GeoInterface.getcoord(::AbstractGeometryTrait, geom::AbstractGeometry, i) =
@@ -60,40 +73,69 @@ function GeoInterface.extent(::AbstractGeometryTrait, geom::AbstractGeometry)
     return Extent(X = (getXMin(env), getXMax(env)), Y = (getYMin(env), getYMax(env)))
 end
 
-function Base.convert(::Type{T}, geom::T) where {T<:AbstractGeometry}
-    return geom
-end
-
-function Base.convert(::Type{T}, geom::X) where {T<:AbstractGeometry,X}
-    return Base.convert(T, GeoInterface.geomtrait(geom), geom)
+GI.convert(::Type{Point}, ::PointTrait, geom::Point; context=nothing) = geom
+function GI.convert(::Type{Point}, ::PointTrait, geom; context=get_global_context())
+    if GI.is3d(geom)
+        return Point(GI.x(geom), GI.y(geom), GI.z(geom), context)
+    else
+        return Point(GI.x(geom), GI.y(geom), context)
+    end
 end
-
-function Base.convert(::Type{Point}, type::PointTrait, geom)
-    return Point(GeoInterface.coordinates(geom))
+GI.convert(::Type{MultiPoint}, ::MultiPointTrait, geom::MultiPoint; context=nothing) = geom
+function GI.convert(::Type{MultiPoint}, t::MultiPointTrait, geom; context=get_global_context())
+    points = Point[GI.convert(Point, PointTrait(), p) for p in GI.getpoint(t, geom)]
+    return MultiPoint(points, context)
 end
-function Base.convert(::Type{MultiPoint}, type::MultiPointTrait, geom)
-    return MultiPoint(GeoInterface.coordinates(geom))
+GI.convert(::Type{LineString}, ::LineStringTrait, geom::LineString; context=nothing) = geom
+function GI.convert(::Type{LineString}, ::LineStringTrait, geom; context=get_global_context())
+    # Faster to make a CoordSeq directly here
+    seq = _geom_to_coord_seq(geom, context)
+    return LineString(createLineString(seq, context), context)
 end
-function Base.convert(::Type{LineString}, type::LineStringTrait, geom)
-    return LineString(GeoInterface.coordinates(geom))
+GI.convert(::Type{LinearRing}, ::LinearRingTrait, geom::LinearRing; context=nothing) = geom
+function GI.convert(::Type{LinearRing}, ::LinearRingTrait, geom; context=get_global_context())
+    # Faster to make a CoordSeq directly here
+    seq = _geom_to_coord_seq(geom, context)
+    return LinearRing(createLinearRing(seq, context), context)
 end
-function Base.convert(::Type{MultiLineString}, type::MultiLineStringTrait, geom)
-    return MultiLineString(GeoInterface.coordinates(geom))
+GI.convert(::Type{MultiLineString}, ::MultiLineStringTrait, geom::MultiLineString; context=nothing) = geom
+function GI.convert(::Type{MultiLineString}, ::MultiLineStringTrait, geom; context=get_global_context())
+    linestrings = LineString[GI.convert(LineString, LineStringTrait(), g; context) for g in getgeom(geom)]
+    return MultiLineString(linestrings)
 end
-function Base.convert(::Type{Polygon}, type::PolygonTrait, geom)
-    return Polygon(GeoInterface.coordinates(geom))
+GI.convert(::Type{Polygon}, ::PolygonTrait, geom::Polygon; context=nothing) = geom
+function GI.convert(::Type{Polygon}, ::PolygonTrait, geom; context=get_global_context())
+    exterior = GI.convert(LinearRing, GI.LinearRingTrait(), GI.getexterior(geom); context)
+    holes = LinearRing[GI.convert(LinearRing, GI.LinearRingTrait(), g; context) for g in GI.gethole(geom)]
+    return Polygon(exterior, holes)
 end
-function Base.convert(::Type{MultiPolygon}, type::MultiPolygonTrait, geom)
-    return MultiPolygon(GeoInterface.coordinates(geom))
+GI.convert(::Type{MultiPolygon}, ::MultiPolygonTrait, geom::MultiPolygon; context=nothing) = geom
+function GI.convert(::Type{MultiPolygon}, ::MultiPolygonTrait, geom; context=get_global_context())
+    polygons = Polygon[GI.convert(Polygon, PolygonTrait(), g; context) for g in GI.getgeom(geom)]
+    return MultiPolygon(polygons)
 end
 
-function Base.convert(t::Type{<:AbstractGeometry}, type::AbstractGeometryTrait, geom)
+function GI.convert(t::Type{<:AbstractGeometry}, ::AbstractGeometryTrait, geom; context=nothing)
     error(
-        "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait to a $t (yet). Please report an issue.",
+        "Cannot convert an object of $(of(geom)) with the $(of()) trait to a $t (yet). Please report an issue.",
     )
-    return f(GeoInterface.coordinates(geom))
 end
 
+function _geom_to_coord_seq(geom, context)
+    npoint = GI.npoint(geom)
+    ndim = GI.is3d(geom) ? 3 : 2
+    seq = createCoordSeq(npoint, context; ndim)
+    for (i, p) in enumerate(GI.getpoint(geom))
+        if ndim == 2
+            setCoordSeq!(seq, i, (GI.x(p), GI.y(p)), context)
+        else
+            setCoordSeq!(seq, i, (GI.x(p), GI.y(p), GI.z(p)), context)
+        end
+    end
+    return seq
+end
+
+
 GeoInterface.distance(
     ::AbstractGeometryTrait,
     ::AbstractGeometryTrait,
@@ -180,3 +222,42 @@ GeoInterface.union(
 ) = union(a, b)
 
 GeoInterfaceRecipes.@enable_geo_plots AbstractGeometry
+
+
+# -----
+# LibGeos operations for any GeoInterface.jl compatible geometries
+# -----
+
+# Internal convert method that avoids the overhead of `convert(LibGEOS, geom)`
+to_geos(geom) = to_geos(GI.geomtrait(geom), geom)
+to_geos(trait, geom) = GI.convert(geointerface_geomtype(trait), trait, geom)
+
+# These methods are all the same with 1 or two geometries, some arguments, and maybe keywords.
+# We define them with `@eval` to avoid all the boilerplate code.
+
+buffer(obj, dist::Real, args...; kw...) = buffer(to_geos(obj), dist::Real, args...; kw...)
+bufferWithStyle(obj, dist::Real; kw...) = bufferWithStyle(to_geos(obj), dist; kw...)
+
+# 1 geom methods
+for f in (
+    :area, :geomLength, :envelope, :minimumRotatedRectangle, :convexhull, :boundary,
+    :unaryUnion, :pointOnSurface, :centroid, :node, :simplify, :topologyPreserveSimplify, :uniquePoints,
+    :delaunayTriangulationEdges, :delaunayTriangulation, :constrainedDelaunayTriangulation,
+)
+    # We convert the geometry to a GEOS geometry and forward it to the geos method
+    @eval $f(geom, args...; kw...) = $f(to_geos(geom), args...; kw...)
+    @eval $f(geom::AbstractGeometry, args...; kw...) =
+        throw(MethodError($f, (geom, args...)))
+end
+
+# 2 geom methods
+for f in (
+    :project, :projectNormalized, :intersection, :difference, :symmetricDifference, :union, :sharedPaths,
+    :snap, :distance, :hausdorffdistance, :nearestPoints, :disjoint, :touches, :intersects, :crosses,
+    :within, :contains, :overlaps, :equalsexact, :covers, :coveredby, :equals,
+)
+    # We convert the geometries to GEOS geometries and forward them to the geos method
+    @eval $f(geom1, geom2, args...; kw...) = $f(to_geos(geom1), to_geos(geom2), args...; kw...)
+    @eval $f(geom1::AbstractGeometry, geom2::AbstractGeometry, args...; kw...) =
+        throw(MethodError($f, (geom1, geom2, args...)))
+end

src/geos_functions.jl

@@ -145,10 +145,11 @@ function getDimensions(ptr::GEOSCoordSeq, context::GEOSContext = get_global_cont
 end
 
 # convenience functions
+# Use Tuple where possible
 function setCoordSeq!(
     ptr::GEOSCoordSeq,
     i::Integer,
-    coords::Vector{Float64},
+    coords::Union{Vector{<:Real},Tuple},
     context::GEOSContext = get_global_context(),
 )
     ndim = length(coords)
@@ -158,9 +159,9 @@ function setCoordSeq!(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
         )
     end
-    setX!(ptr, i, coords[1], context)
-    setY!(ptr, i, coords[2], context)
-    ndim >= 3 && setZ!(ptr, i, coords[3], context)
+    setX!(ptr, i, Float64(coords[1]), context)
+    setY!(ptr, i, Float64(coords[2]), context)
+    ndim >= 3 && setZ!(ptr, i, Float64(coords[3]), context)
     ptr
 end
 
@@ -233,7 +234,7 @@ end
 
 function getX(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ncoords = getSize(ptr, context)
-    xcoords = Array{Float64}(undef, ncoords)
+    xcoords = Vector{Float64}(undef, ncoords)
     start = pointer(xcoords)
     floatsize = sizeof(Float64)
     for i = 0:ncoords-1
@@ -258,7 +259,7 @@ end
 
 function getY(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ncoords = getSize(ptr, context)
-    ycoords = Array{Float64}(undef, ncoords)
+    ycoords = Vector{Float64}(undef, ncoords)
     start = pointer(ycoords)
     floatsize = sizeof(Float64)
     for i = 0:ncoords-1
@@ -777,7 +778,10 @@ function within(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_conte
     result != 0x00
 end
 
-function Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1))
+Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1)) =
+    contains(obj1, obj2, context)
+
+function contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1))
     result = GEOSContains_r(context, obj1, obj2)
     if result == 0x02
         error("LibGEOS: Error in GEOSContains")
@@ -842,7 +846,10 @@ function destroyPreparedGeom(obj::PreparedGeometry, context::GEOSContext = get_g
     GEOSPreparedGeom_destroy_r(context, obj)
 end
 
-function Base.contains(
+Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1)) =
+    contains(obj1, obj2, context)
+
+function contains(
     obj1::PreparedGeometry,
     obj2::Geometry,
     context::GEOSContext = get_context(obj1)

src/geos_types.jl

@@ -89,12 +89,17 @@ mutable struct LineString <: AbstractGeometry
         finalizer(destroyGeom, line)
         line
     end
-    #create a linestring from a list of coordinates
+    # create a linestring from a vector of points
     function LineString(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context())
         line = new(createLineString(coords, context), context)
         finalizer(destroyGeom, line)
         line
     end
+    function LineString(coords::Vector{Point}, context::GEOSContext = get_global_context())
+        line = new(createLineString(coords, context), context)
+        finalizer(destroyGeom, line)
+        line
+    end
 end
 
 mutable struct MultiLineString <: AbstractGeometry
@@ -124,6 +129,13 @@ mutable struct MultiLineString <: AbstractGeometry
                 GEOSGeom[createLineString(coords, context) for coords in multiline],
                 context),
             context)
+    MultiLineString(multiline::Vector{LineString}, context::GEOSContext = get_global_context()) =
+        MultiLineString(
+            createCollection(
+                GEOS_MULTILINESTRING,
+                GEOSGeom[ls.ptr for ls in multiline],
+                context),
+            context)
 end
 
 mutable struct LinearRing <: AbstractGeometry
@@ -256,6 +268,13 @@ mutable struct GeometryCollection <: AbstractGeometry
                 collection,
                 context),
             context)
+    GeometryCollection(collection::Vector{<:AbstractGeometry}, context::GEOSContext = get_global_context()) =
+        GeometryCollection(
+            createCollection(
+                GEOS_GEOMETRYCOLLECTION,
+                GEOSGeom[geom.ptr for geom in collection],
+                context),
+            context)
 end
 
 const Geometry = Union{
@@ -434,7 +453,7 @@ typesalt(::Type{Polygon}            ) = 0xa5c895d62ef56723
 
 function Base.hash(geo::AbstractGeometry, h::UInt)::UInt
     h = hash(typesalt(typeof(geo)), h)
-    if has_coord_seq(geo) 
+    if has_coord_seq(geo)
         return hash_coord_seq(geo, h)
     else
         for i in 1:ngeom(geo)