diff --git a/src/geointerface.jl b/src/geointerface.jl
index c8ebfd3a..61a6a171 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -10,6 +10,7 @@ const lookup_method = Dict{DataType,Function}(
     GeoInterface.MultiLineStringTrait => createmultilinestring,
     GeoInterface.PolygonTrait => createpolygon,
     GeoInterface.MultiPolygonTrait => createmultipolygon,
+    GeoInterface.GeometryCollectionTrait => creategeomcollection,
 )
 
 let pointtypes = (wkbPoint, wkbPoint25D, wkbPointM, wkbPointZM),
@@ -280,18 +281,25 @@ let pointtypes = (wkbPoint, wkbPoint25D, wkbPointM, wkbPointZM),
         return toWKT(geom)
     end
 
-    function GeoInterface.convert(
-        ::Type{T},
-        type::GeometryTraits,
-        geom,
-    ) where {T<:IGeometry}
-        f = get(lookup_method, typeof(type), nothing)
-        isnothing(f) && error(
-            "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
-        )
-        return f(GeoInterface.coordinates(geom))
+    function GeoInterface.convert(::Type{T}, trait::GeometryTraits, geom) where {T<:IGeometry}
+        f = get(lookup_method, typeof(trait), nothing)
+        isnothing(f) && _convert_error(geom, trait)
+        coords = GeoInterface.coordinates(geom)
+        return f(coords)
     end
 
+    function GeoInterface.convert(::Type{T}, trait::GeoInterface.GeometryCollectionTrait, collection) where {T<:IGeometry}
+        gdal_collection = creategeomcollection()
+        for geom in GeoInterface.getgeom(collection)
+            addgeom!(gdal_collection, to_gdal(geom))
+        end
+        return gdal_collection
+    end
+
+    @noinline _convert_error(geom, trait) = error(
+        "Cannot convert an object of $(typeof(geom)) with the $(typeof(trait)) trait (yet). Please report an issue at https://github.com/yeesian/ArchGDAL.jl if you need this",
+    )
+
     function GeoInterface.geomtrait(
         geom::Union{map(T -> AbstractGeometry{T}, pointtypes)...},
     )
diff --git a/src/ogr/geometry.jl b/src/ogr/geometry.jl
index a8471748..01e27410 100644
--- a/src/ogr/geometry.jl
+++ b/src/ogr/geometry.jl
@@ -1,3 +1,15 @@
+# Local convert method converts any GeoInterface geometry
+# to the equivalent ArchGDAL geometry.
+# `geom` can be points, so this is a hot path that needs to be fast.
+function to_gdal(geom)
+    GeoInterface.isgeometry(geom) || _not_a_geom_error()
+    trait = GeoInterface.geomtrait(geom)
+    typ = geointerface_geomtype(trait)
+    return GeoInterface.convert(typ, trait, geom)
+end
+
+@noinline _not_a_geom_error() =
+    throw(ArgumentError("Ojbect of type $(typeof(geom)) is not a GeoInterface compatible geometry"))
 
 """
     fromWKB(data)
@@ -110,6 +122,7 @@ function clone(geom::AbstractGeometry{T}) where {T}
         return IGeometry{T}(GDAL.ogr_g_clone(geom))
     end
 end
+clone(geom) = clone(to_gdal(geom))
 
 function unsafe_clone(geom::AbstractGeometry{T}) where {T}
     if geom.ptr == C_NULL
@@ -171,6 +184,7 @@ with the geometry multiple times, by storing caches of calculated geometry infor
 function preparegeom(geom::AbstractGeometry{T}) where {T}
     return IPreparedGeometry{T}(GDAL.ogrcreatepreparedgeometry(geom))
 end
+preparegeom(geom) = preparegeom(to_gdal(geom))
 
 function unsafe_preparegeom(geom::AbstractGeometry{T}) where {T}
     return PreparedGeometry{T}(GDAL.ogrcreatepreparedgeometry(geom))
@@ -263,6 +277,7 @@ function envelope(geom::AbstractGeometry)::GDAL.OGREnvelope
     GDAL.ogr_g_getenvelope(geom, envelope)
     return envelope[]
 end
+envelope(geom) = envelope(to_gdal(geom))
 
 """
     envelope3d(geom::AbstractGeometry)
@@ -274,6 +289,7 @@ function envelope3d(geom::AbstractGeometry)::GDAL.OGREnvelope3D
     GDAL.ogr_g_getenvelope3d(geom, envelope)
     return envelope[]
 end
+envelope3d(geom) = envelope3d(to_gdal(geom))
 
 """
     boundingbox(geom::AbstractGeometry)
@@ -293,6 +309,7 @@ function boundingbox(geom::AbstractGeometry)::IGeometry
         [MinX, MaxY],
     ])
 end
+boundingbox(geom) = boundingbox(to_gdal(geom))
 
 """
     toWKB(geom::AbstractGeometry, order::OGRwkbByteOrder = wkbNDR)
@@ -312,6 +329,7 @@ function toWKB(
     @ogrerr result "Failed to export geometry to WKB"
     return buffer
 end
+toWKB(geom, order::OGRwkbByteOrder=wkbNDR) = toWKB(to_gdal(geom), order)
 
 """
     toISOWKB(geom::AbstractGeometry, order::OGRwkbByteOrder = wkbNDR)
@@ -331,6 +349,7 @@ function toISOWKB(
     @ogrerr result "Failed to export geometry to ISO WKB"
     return buffer
 end
+toISOWKB(geom, order::OGRwkbByteOrder=wkbNDR) = toISOWKB(to_gdal(geom), order)
 
 """
     wkbsize(geom::AbstractGeometry)
@@ -338,6 +357,7 @@ end
 Returns size (in bytes) of related binary representation.
 """
 wkbsize(geom::AbstractGeometry)::Integer = GDAL.ogr_g_wkbsize(geom)
+wkbsize(geom) = wkbsize(to_gdal(geom))
 
 """
     toWKT(geom::AbstractGeometry)
@@ -352,6 +372,7 @@ function toWKT(geom::AbstractGeometry)::String
     GDAL.vsifree(pointer(wkt_ptr[]))
     return wkt
 end
+toWKT(geom) = toWKT(to_gdal(geom))
 
 """
     toISOWKT(geom::AbstractGeometry)
@@ -366,6 +387,7 @@ function toISOWKT(geom::AbstractGeometry)::String
     GDAL.vsifree(pointer(isowkt_ptr[]))
     return wkt
 end
+toISOWKT(geom) = toISOWKT(to_gdal(geom))
 
 """
     getgeomtype(geom::AbstractGeometry)
@@ -386,6 +408,7 @@ function geomname(geom::AbstractGeometry)::Union{String,Missing}
         GDAL.ogr_g_getgeometryname(geom)
     end
 end
+geomname(geom) = geomname(to_gdal(geom))
 
 """
     flattento2d!(geom::AbstractGeometry)
@@ -400,6 +423,18 @@ function flattento2d!(geom::G)::G where {G<:AbstractGeometry}
     return geom
 end
 
+"""
+    flattento2d(geom)
+
+Convert any GeoInterface compatible geometry to strictly 2D.
+
+A new geometry object will always be created.
+"""
+function flattento2d(geom::G)::G where {G<:AbstractGeometry}
+    flattento2d(clone(geom))
+end
+flattento2d(geom) = flattento2d!(to_gdal(geom))
+
 """
     closerings!(geom::AbstractGeometry)
 
@@ -413,6 +448,18 @@ function closerings!(geom::G)::G where {G<:AbstractGeometry}
     return geom
 end
 
+"""
+    closerings(geom)
+
+Convert any GeoInterface compatible geometry to strictly 2D.
+
+A new geometry object will always be created.
+"""
+function closerings(geom::G)::G where {G<:AbstractGeometry}
+    closerings(clone(geom))
+end
+closerings(geom) = closerings!(to_gdal(geom))
+
 """
     fromGML(data)
 
@@ -430,19 +477,21 @@ fromGML(data)::IGeometry = IGeometry(GDAL.ogr_g_createfromgml(data))
 unsafe_fromGML(data)::Geometry = Geometry(GDAL.ogr_g_createfromgml(data))
 
 """
-    toGML(geom::AbstractGeometry)
+    toGML(geom)
 
-Convert a geometry into GML format.
+Convert any GeoInterface compatible geometry into GML format.
 """
 toGML(geom::AbstractGeometry)::String = GDAL.ogr_g_exporttogml(geom)
+toGML(geom) = toGML(to_gdal(geom))
 
 """
-    toKML(geom::AbstractGeometry, altitudemode = C_NULL)
+    toKML(geom, altitudemode = C_NULL)
 
-Convert a geometry into KML format.
+Convert any GeoInterface compatible geometry into KML format.
 """
 toKML(geom::AbstractGeometry, altitudemode = C_NULL)::String =
     GDAL.ogr_g_exporttokml(geom, altitudemode)
+toKML(geom) = toKML(to_gdal(geom))
 # ↑ * `altitudemode`: value to write in altitudeMode element, or NULL.
 
 """
@@ -467,9 +516,11 @@ A GeoJSON fragment or NULL in case of error.
 """
 toJSON(geom::AbstractGeometry; kwargs...)::String =
     GDAL.ogr_g_exporttojsonex(geom, String["$k=$v" for (k, v) in kwargs])
+toJSON(geom; kwargs...) = toJSON(to_gdal(geom); kwargs...)
 
 toJSON(geom::AbstractGeometry, options::Vector{String})::String =
     GDAL.ogr_g_exporttojsonex(geom, options)
+toJSON(geom, options::Vector{String}) = toJSON(to_gdal(geom), options)
 
 """
     fromJSON(data::String)
@@ -550,6 +601,22 @@ function transform!(
     return geom
 end
 
+"""
+    transform(geom, coordtransform::CoordTransform)
+
+Apply arbitrary coordinate transformation to geometry,
+always creating a new object.
+
+### Parameters
+* `geom`: any GeoInterface compatible geometry.
+* `coordtransform`: transformation to apply.
+"""
+function transform(geom::G, coordtransform::CoordTransform)::G where {G<:AbstractGeometry}
+    transform!(clone(geom), coordtransform)
+end
+transform(geom, coordtransform::CoordTransform) =
+    transform!(to_gdal(geom), coordtransform)
+
 # """
 # Transform geometry to new spatial reference system.
 
@@ -590,14 +657,15 @@ Compute a simplified geometry.
 """
 simplify(geom::AbstractGeometry, tol::Real)::IGeometry =
     IGeometry(GDAL.ogr_g_simplify(geom, tol))
+simplify(geom, tol::Real) = simplify(to_gdal(geom), tol)
 
 unsafe_simplify(geom::AbstractGeometry, tol::Real)::Geometry =
     Geometry(GDAL.ogr_g_simplify(geom, tol))
 
 """
-    simplifypreservetopology(geom::AbstractGeometry, tol::Real)
+    simplifypreservetopology(geom, tol::Real)
 
-Simplify the geometry while preserving topology.
+Simplify a GeoInterface compatible geometry while preserving topology.
 
 ### Parameters
 * `geom`: the geometry.
@@ -605,12 +673,13 @@ Simplify the geometry while preserving topology.
 """
 simplifypreservetopology(geom::AbstractGeometry, tol::Real)::IGeometry =
     IGeometry(GDAL.ogr_g_simplifypreservetopology(geom, tol))
+simplifypreservetopology(geom, tol::Real) = simplifypreservetopology(to_gdal(geom), tol)
 
 unsafe_simplifypreservetopology(geom::AbstractGeometry, tol::Real)::Geometry =
     Geometry(GDAL.ogr_g_simplifypreservetopology(geom, tol))
 
 """
-    delaunaytriangulation(geom::AbstractGeometry, tol::Real, onlyedges::Bool)
+    delaunaytriangulation(geom, tol::Real, onlyedges::Bool)
 
 Return a Delaunay triangulation of the vertices of the geometry.
 
@@ -627,6 +696,7 @@ function delaunaytriangulation(
 )::IGeometry
     return IGeometry(GDAL.ogr_g_delaunaytriangulation(geom, tol, onlyedges))
 end
+delaunaytriangulation(geom, tol::Real, onlyedges) = delaunaytriangulation(to_gdal(geom), tol, onlyedges)
 
 function unsafe_delaunaytriangulation(
     geom::AbstractGeometry,
@@ -654,7 +724,24 @@ function segmentize!(geom::G, maxlength::Real)::G where {G<:AbstractGeometry}
 end
 
 """
-    intersects(g1::AbstractGeometry, g2::AbstractGeometry)
+    segmentize(geom, maxlength::Real)
+
+Modify a GeoInterface compatible geometry such it has no segment longer than the given distance.
+
+Interpolated points will have Z and M values (if needed) set to 0. Distance
+computation is performed in 2d only
+
+### Parameters
+* `geom`: the geometry to segmentize
+* `maxlength`: the maximum distance between 2 points after segmentization
+"""
+function segmentize(geom::G, maxlength::Real)::G where {G<:AbstractGeometry}
+    segmentize!(clone(geom), maxlength)
+end
+segmentize(geom, maxlength::Real) = segmentize!(to_gdal(geom), maxlength)
+
+"""
+    intersects(g1, g2)
 
 Returns whether the geometries intersect
 
@@ -664,75 +751,82 @@ boxes) of the two geometries overlap.
 """
 intersects(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_intersects(g1, g2))
+intersects(g1, g2)::Bool = intersects(to_gdal(g1), to_gdal(g2))
 
 intersects(g1::AbstractPreparedGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogrpreparedgeometryintersects(g1, g2))
 
 """
-    equals(g1::AbstractGeometry, g2::AbstractGeometry)
+    equals(g1, g2)
 
 Returns `true` if the geometries are equivalent.
 """
 equals(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_equals(g1, g2))
+equals(g1, g2)::Bool = equals(to_gdal(g1), to_gdal(g2))
 
 function Base.:(==)(g1::AbstractGeometry, g2::AbstractGeometry)
     return equals(g1, g2)
 end
 
 """
-    disjoint(g1::AbstractGeometry, g2::AbstractGeometry)
+    disjoint(g1, g2)
 
 Returns `true` if the geometries are disjoint.
 """
 disjoint(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_disjoint(g1, g2))
+disjoint(g1, g2)::Bool = disjoint(to_gdal(g1), to_gdal(g2))
 
 """
-    touches(g1::AbstractGeometry, g2::AbstractGeometry)
+    touches(g1, g2)
 
 Returns `true` if the geometries are touching.
 """
 touches(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_touches(g1, g2))
+touches(g1, g2)::Bool = touches(to_gdal(g1), to_gdal(g2))
 
 """
-    crosses(g1::AbstractGeometry, g2::AbstractGeometry)
+    crosses(g1, g2)
 
 Returns `true` if the geometries are crossing.
 """
 crosses(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_crosses(g1, g2))
+crosses(g1, g2)::Bool = crosses(to_gdal(g1), to_gdal(g2))
 
 """
-    within(g1::AbstractGeometry, g2::AbstractGeometry)
+    within(g1, g2)
 
 Returns `true` if g1 is contained within g2.
 """
 within(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_within(g1, g2))
+within(g1, g2)::Bool = within(to_gdal(g1), to_gdal(g2))
 
 """
-    contains(g1::AbstractGeometry, g2::AbstractGeometry)
+    contains(g1, g2)
 
 Returns `true` if g1 contains g2.
 """
 contains(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_contains(g1, g2))
-
 contains(g1::AbstractPreparedGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogrpreparedgeometrycontains(g1, g2))
+contains(g1, g2)::Bool = contains(to_gdal(g1), to_gdal(g2))
 
 """
-    overlaps(g1::AbstractGeometry, g2::AbstractGeometry)
+    overlaps(g1, g2)
 
 Returns `true` if the geometries overlap.
 """
 overlaps(g1::AbstractGeometry, g2::AbstractGeometry)::Bool =
     Bool(GDAL.ogr_g_overlaps(g1, g2))
+overlaps(g1, g2)::Bool = overlaps(to_gdal(g1), to_gdal(g2))
 
 """
-    boundary(geom::AbstractGeometry)
+    boundary(geom)
 
 Returns the boundary of the geometry.
 
@@ -741,12 +835,13 @@ geometry on which the method is invoked.
 """
 boundary(geom::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_boundary(geom))
+boundary(geom) = boundary(to_gdal(geom))
 
 unsafe_boundary(geom::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_boundary(geom))
 
 """
-    convexhull(geom::AbstractGeometry)
+    convexhull(geom)
 
 Returns the convex hull of the geometry.
 
@@ -755,12 +850,13 @@ geometry on which the method is invoked.
 """
 convexhull(geom::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_convexhull(geom))
+convexhull(geom) = convexhull(to_gdal(geom))
 
 unsafe_convexhull(geom::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_convexhull(geom))
 
 """
-    buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30)
+    buffer(geom, dist::Real, quadsegs::Integer = 30)
 
 Compute buffer of geometry.
 
@@ -784,6 +880,8 @@ accuracy of the result.
 """
 buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer = 30)::IGeometry =
     IGeometry(GDAL.ogr_g_buffer(geom, dist, quadsegs))
+buffer(geom, dist::Real, quadsegs::Integer = 30) =
+    buffer(to_gdal(geom), dist, quadsegs)
 
 function unsafe_buffer(
     geom::AbstractGeometry,
@@ -794,7 +892,7 @@ function unsafe_buffer(
 end
 
 """
-    intersection(g1::AbstractGeometry, g2::AbstractGeometry)
+    intersection(g1, g2)
 
 Returns a new geometry representing the intersection of the geometries, or NULL
 if there is no intersection or an error occurs.
@@ -805,6 +903,7 @@ two geometries intersect.
 """
 intersection(g1::AbstractGeometry, g2::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_intersection(g1, g2))
+intersection(g1, g2) = intersection(to_gdal(g1), to_gdal(g2))
 
 unsafe_intersection(g1::AbstractGeometry, g2::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_intersection(g1, g2))
@@ -816,6 +915,7 @@ Returns a new geometry representing the union of the geometries.
 """
 union(g1::AbstractGeometry, g2::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_union(g1, g2))
+union(g1, g2) = union(to_gdal(g1), to_gdal(g2))
 
 unsafe_union(g1::AbstractGeometry, g2::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_union(g1, g2))
@@ -832,12 +932,13 @@ multisurfaces (multipolygons).
 """
 pointonsurface(geom::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_pointonsurface(geom))
+pointonsurface(g) = pointonsurface(to_gdal(g))
 
 unsafe_pointonsurface(geom::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_pointonsurface(geom))
 
 """
-    difference(g1::AbstractGeometry, g2::AbstractGeometry)
+    difference(g1, g2)
 
 Generates a new geometry which is the region of this geometry with the region of
 the other geometry removed.
@@ -848,43 +949,53 @@ if the difference is empty.
 """
 difference(g1::AbstractGeometry, g2::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_difference(g1, g2))
+difference(g1, g2) = difference(to_gdal(g1), to_gdal(g2))
 
 unsafe_difference(g1::AbstractGeometry, g2::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_difference(g1, g2))
 
 """
-    symdifference(g1::AbstractGeometry, g2::AbstractGeometry)
+    symdifference(g1, g2)
 
 Returns a new geometry representing the symmetric difference of the geometries
 or NULL if the difference is empty or an error occurs.
 """
 symdifference(g1::AbstractGeometry, g2::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_symdifference(g1, g2))
+symdifference(g1, g2) = symdifference(to_gdal(g1), to_gdal(g2))
 
 unsafe_symdifference(g1::AbstractGeometry, g2::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_symdifference(g1, g2))
 
 """
-    distance(g1::AbstractGeometry, g2::AbstractGeometry)
+    distance(g1, g2)
 
 Returns the distance between the geometries or -1 if an error occurs.
 """
 distance(g1::AbstractGeometry, g2::AbstractGeometry)::Float64 =
     GDAL.ogr_g_distance(g1, g2)
+distance(g1, g2) = distance(to_gdal(g1), to_gdal(g2))
 
 """
-    geomlength(geom::AbstractGeometry)
+    geomlength(geom)
 
 Returns the length of the geometry, or 0.0 for unsupported geometry types.
 """
 geomlength(geom::AbstractGeometry)::Float64 = GDAL.ogr_g_length(geom)
+geomlength(geom) = geomlength(to_gdal(geom))
+# TODO Add `length` method to match GeoInterface naming conventions?
+# this will mean using `Base.length` everywhere
 
 """
-    geomarea(geom::AbstractGeometry)
+    geomarea(geom)
 
 Returns the area of the geometry or 0.0 for unsupported geometry types.
 """
 geomarea(geom::AbstractGeometry)::Float64 = GDAL.ogr_g_area(geom)
+geomarea(geom) = geomarea(to_gdal(geom))
+
+# Add `area` method to match GeoInterface naming conventions
+area(geom) = geomarea(geom)
 
 """
     centroid!(geom::AbstractGeometry, centroid::AbstractGeometry)
@@ -910,7 +1021,7 @@ function centroid!(
 end
 
 """
-    centroid(geom::AbstractGeometry)
+    centroid(geom)
 
 Compute the geometry centroid.
 
@@ -928,6 +1039,7 @@ function centroid(geom::AbstractGeometry)::IGeometry
     centroid!(geom, point)
     return point
 end
+centroid(geom) = centroid(to_gdal(geom))
 
 function unsafe_centroid(geom::AbstractGeometry)::Geometry
     point = unsafe_createpoint()
@@ -936,7 +1048,7 @@ function unsafe_centroid(geom::AbstractGeometry)::Geometry
 end
 
 """
-    pointalongline(geom::AbstractGeometry, distance::Real)
+    pointalongline(geom, distance::Real)
 
 Fetch point at given distance along curve.
 
@@ -950,6 +1062,7 @@ a point or NULL.
 """
 pointalongline(geom::AbstractGeometry, distance::Real)::IGeometry =
     IGeometry(GDAL.ogr_g_value(geom, distance))
+pointalongline(geom, distance::Real) = pointalongline(to_gdal(geom), distance)
 
 unsafe_pointalongline(geom::AbstractGeometry, distance::Real)::Geometry =
     Geometry(GDAL.ogr_g_value(geom, distance))
@@ -1113,6 +1226,7 @@ ismeasured(geom::_AbstractGeometryNoM) = false
 Returns `true` if the geometry has no points, otherwise `false`.
 """
 isempty(geom::AbstractGeometry)::Bool = GDAL.ogr_g_isempty(geom) != 0
+isempty(geom)::Bool = isempty(to_gdal(geom))
 
 """
     isvalid(geom::AbstractGeometry)
@@ -1120,6 +1234,7 @@ isempty(geom::AbstractGeometry)::Bool = GDAL.ogr_g_isempty(geom) != 0
 Returns `true` if the geometry is valid, otherwise `false`.
 """
 isvalid(geom::AbstractGeometry)::Bool = Bool(GDAL.ogr_g_isvalid(geom))
+isvalid(geom)::Bool = isvalid(to_gdal(geom))
 
 """
     issimple(geom::AbstractGeometry)
@@ -1127,6 +1242,7 @@ isvalid(geom::AbstractGeometry)::Bool = Bool(GDAL.ogr_g_isvalid(geom))
 Returns `true` if the geometry is simple, otherwise `false`.
 """
 issimple(geom::AbstractGeometry)::Bool = Bool(GDAL.ogr_g_issimple(geom))
+issimple(geom)::Bool = issimple(to_gdal(geom))
 
 """
     isring(geom::AbstractGeometry)
@@ -1134,6 +1250,7 @@ issimple(geom::AbstractGeometry)::Bool = Bool(GDAL.ogr_g_issimple(geom))
 Returns `true` if the geometry is a ring, otherwise `false`.
 """
 isring(geom::AbstractGeometry)::Bool = Bool(GDAL.ogr_g_isring(geom))
+isring(geom)::Bool = isring(to_gdal(geom))
 
 """
     polygonize(geom::AbstractGeometry)
@@ -1147,6 +1264,7 @@ impossible due to topological inconsistencies.
 """
 polygonize(geom::AbstractGeometry)::IGeometry =
     IGeometry(GDAL.ogr_g_polygonize(geom))
+polygonize(geom) = polygonize(to_gdal(geom))
 
 unsafe_polygonize(geom::AbstractGeometry)::Geometry =
     Geometry(GDAL.ogr_g_polygonize(geom))
@@ -1515,6 +1633,8 @@ function removegeom!(
     return geom
 end
 
+removegeom(geom, i::Integer) = removegeom!(clone(geom), i, true) 
+
 """
     removeallgeoms!(geom::AbstractGeometry, todelete::Bool = true)
 
@@ -1563,7 +1683,6 @@ MULTICURVE or MULTISURFACE in it, by approximating curve geometries.
 * `options`: options as a null-terminated list of strings or NULL.
     See OGRGeometryFactory::curveToLineString() for valid options.
 """
-
 function lineargeom(
     geom::AbstractGeometry,
     stepsize::Real = 0;
diff --git a/src/spatialref.jl b/src/spatialref.jl
index 087f26e4..14dab5a9 100644
--- a/src/spatialref.jl
+++ b/src/spatialref.jl
@@ -169,6 +169,21 @@ function reproject(
     end
 end
 
+function reproject(
+    geom,
+    sourcecrs::GFT.GeoFormat,
+    targetcrs::GFT.GeoFormat;
+    kwargs...,
+)
+    if GI.isgeometry(geom)
+        reproject(to_gdal(geom), sourcecrs, targetcrs; kwargs...)
+    elseif geom isa AbstractArray && (length(geom) > 0) && GI.isgeometry(first(geom))
+        reproject(to_gdal.(geom), Ref(sourcecrs), Ref(targetcrs); kwargs...)
+    else
+        throw(ArgumentError("geom is not a GeoInterface compatible geometry"))
+    end
+end
+
 """
     crs2transform(f::Function, sourcecrs::GeoFormat, targetcrs::GeoFormat;
         kwargs...)
diff --git a/test/test_geometry.jl b/test/test_geometry.jl
index fdb85c5a..3d5f1b53 100644
--- a/test/test_geometry.jl
+++ b/test/test_geometry.jl
@@ -37,6 +37,7 @@ import GeoFormatTypes as GFT
     @testset "Create a Point" begin
         # Method 1
         AG.createpoint(100, 70) do point
+            wrapped_point = GeoInterface.convert(GeoInterface, point)
             @test point isa AG.Geometry{AG.wkbPoint}
             @test isapprox(GI.coordinates(point), [100, 70], atol = 1e-6)
             @test AG.geomdim(point) == 0
@@ -115,14 +116,18 @@ import GeoFormatTypes as GFT
                 0x00,
             ]
             @test AG.toKML(point, "relativeToGround") ==
+                  AG.toKML(wrapped_point, "relativeToGround") ==
                   "<Point><altitudeMode>relativeToGround</altitudeMode>" *
                   "<coordinates>100,70,0</coordinates></Point>"
             @test AG.toKML(point, "clampToGround") ==
+                  AG.toKML(wrapped_point, "clampToGround") ==
                   "<Point><altitudeMode>clampToGround</altitudeMode>" *
                   "<coordinates>100,70,0</coordinates></Point>"
             @test AG.toKML(point) ==
+                  AG.toKML(wrapped_point) ==
                   "<Point><coordinates>100,70,0</coordinates></Point>"
             @test AG.toJSON(point) ==
+                  AG.toJSON(wrapped_point) ==
                   "{ \"type\": \"Point\", \"coordinates\": " *
                   "[ 100.0, 70.0, 0.0 ] }"
             @test startswith(
@@ -165,20 +170,31 @@ import GeoFormatTypes as GFT
 
         # Method 2
         point = AG.createpoint(100, 70)
+        tuple_point = (100, 70)
         @test AG.geomdim(point) == 0
+        @test AG.geomdim(tuple_point) == 0
         @test AG.getcoorddim(point) == 2
         AG.setcoorddim!(point, 3)
         @test AG.getcoorddim(point) == 3
         @test AG.isvalid(point) == true
+        @test AG.isvalid(tuple_point) == true
         @test AG.issimple(point) == true
+        @test AG.issimple(tuple_point) == true
         @test AG.isring(point) == false
+        @test AG.isring(tuple_point) == false
         @test AG.getz(point, 0) == 0
         @test typeof(point) == AG.IGeometry{AG.wkbPoint}
         @test sprint(print, AG.envelope(point)) ==
+              sprint(print, AG.envelope(tuple_point)) ==
+              "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
+        @test sprint(print, AG.envelope(tuple_point)) ==
               "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
         @test sprint(print, AG.envelope3d(point)) ==
               "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
-        @test AG.toISOWKB(point, AG.wkbNDR) == UInt8[
+        @test sprint(print, AG.envelope3d(tuple_point)) ==
+              "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
+        @test AG.toISOWKB(point, AG.wkbNDR) ==
+              AG.toISOWKB(point, tuple_point) == UInt8[
             0x01,
             0xe9,
             0x03,
@@ -209,7 +225,8 @@ import GeoFormatTypes as GFT
             0x00,
             0x00,
         ]
-        @test AG.toISOWKB(point, AG.wkbXDR) == UInt8[
+        @test AG.toISOWKB(point, AG.wkbXDR) ==
+              AG.toISOWKB(tuple_point, AG.wkbXDR) == UInt8[
             0x00,
             0x00,
             0x00,
@@ -241,17 +258,22 @@ import GeoFormatTypes as GFT
             0x00,
         ]
         @test AG.toKML(point, "relativeToGround") ==
+              AG.toKML(tuple_point, "relativeToGround") ==
               "<Point><altitudeMode>relativeToGround</altitudeMode>" *
               "<coordinates>100,70,0</coordinates></Point>"
         @test AG.toKML(point, "clampToGround") ==
+              AG.toKML(tuple_point, "clampToGround") ==
               "<Point><altitudeMode>clampToGround</altitudeMode>" *
               "<coordinates>100,70,0</coordinates></Point>"
         @test AG.toKML(point) ==
+              AG.toKML(tuple_point) ==
               "<Point><coordinates>100,70,0</coordinates></Point>"
         @test AG.toJSON(point) ==
+              AG.toJSON(tuple_point) ==
               "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0, 0.0 ] }"
         @test AG.equals(point, AG.createpoint(100, 70, 0)) == true
         @test AG.equals(point, AG.createpoint((100, 70, 0))) == true
+        @test AG.equals(tuple_point, AG.createpoint(100, 70, 0)) == true
         AG.flattento2d!(point)
         @test AG.getcoorddim(point) == 2
         @test AG.getnonlineargeomflag() == true
@@ -260,7 +282,9 @@ import GeoFormatTypes as GFT
         AG.setnonlineargeomflag!(true)
         @test AG.getnonlineargeomflag() == true
         AG.closerings!(point)
+        closed_tuple = AG.closerings(tuple_point)
         @test AG.toJSON(point) ==
+              AG.toJSON(closed_tuple) ==
               "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
     end
 
@@ -673,19 +697,27 @@ import GeoFormatTypes as GFT
             [4.0, 5.0, 6.0],
             [7.0, 8.0, 9.0],
         )
+        wrapped_geom1 = GeoInterface.convert(GeoInterface, geom1)
+        wrapped_geom2 = GeoInterface.convert(GeoInterface, geom2)
 
         AG.closerings!(geom1)
         @test AG.disjoint(geom1, geom2) == false
+        @test AG.disjoint(wrapped_geom1, wrapped_geom2) == false
         @test AG.touches(geom1, geom2) == true
+        @test AG.touches(wrapped_geom1, wrapped_geom2) == true
         @test AG.crosses(geom1, geom2) == false
+        @test AG.crosses(wrapped_geom1, wrapped_geom2) == false
         @test AG.overlaps(geom1, geom2) == false
+        @test AG.overlaps(wrapped_geom1, wrapped_geom2) == false
 
         @test AG.toWKT(AG.boundary(geom2)) == "GEOMETRYCOLLECTION EMPTY"
+        @test AG.toWKT(AG.boundary(wrapped_geom2)) == "GEOMETRYCOLLECTION EMPTY"
         AG.boundary(geom2) do result
             @test AG.toWKT(result) == "GEOMETRYCOLLECTION EMPTY"
         end
 
         @test AG.toWKT(AG.union(geom1, geom2)) ==
+              AG.toWKT(AG.union(wrapped_geom1, wrapped_geom2)) ==
               "GEOMETRYCOLLECTION (" *
               "POLYGON (" *
               "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
@@ -709,6 +741,7 @@ import GeoFormatTypes as GFT
         end
 
         @test AG.toWKT(AG.difference(geom1, geom2)) ==
+              AG.toWKT(AG.difference(wrapped_geom1, wrapped_geom2)) ==
               "MULTIPOLYGON (" *
               "((0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
               "(3 1 8,3 3 8,1 3 8,1 1 8,3 1 8))," *
@@ -769,6 +802,7 @@ import GeoFormatTypes as GFT
         end
 
         @test AG.toWKT(AG.symdifference(geom1, geom2)) ==
+              AG.toWKT(AG.symdifference(wrapped_geom1, wrapped_geom2)) ==
               "GEOMETRYCOLLECTION (" *
               "POLYGON (" *
               "(0 4 8,4 4 8,4 0 8,0 0 8,0 4 8)," *
@@ -813,6 +847,7 @@ import GeoFormatTypes as GFT
             "(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8))," *
             "POINT EMPTY)",
         )
+        wrapped_geom3 = GeoInterface.convert(GeoInterface, geom3)
         AG.clone(geom3) do geom4
             @test sprint(print, AG.clone(geom3)) ==
                   "Geometry: GEOMETRYCOLLECTION (" *
@@ -826,11 +861,25 @@ import GeoFormatTypes as GFT
                   " ... MPTY)"
             @test typeof(geom4) == AG.Geometry{AG.wkbGeometryCollection25D}
         end
+        AG.clone(wrapped_geom3) do geom4
+            @test sprint(print, AG.clone(geom3)) ==
+                  "Geometry: GEOMETRYCOLLECTION (" *
+                  "POINT (2 5 8)," *
+                  "POLYGON ((0 0 8," *
+                  " ... MPTY)"
+            @test sprint(print, AG.clone(geom4)) ==
+                  "Geometry: GEOMETRYCOLLECTION (" *
+                  "POINT (2 5 8)," *
+                  "POLYGON ((0 0 8," *
+                  " ... MPTY)"
+            @test typeof(geom4) == AG.Geometry{AG.wkbGeometryCollection25D}
+        end
         AG.clone(AG.getgeom(geom3, 3)) do geom4
             @test sprint(print, geom4) == "Geometry: POINT EMPTY"
         end
 
         @test AG.toISOWKT(geom3) ==
+              AG.toISOWKT(wrapped_geom3) ==
               "GEOMETRYCOLLECTION Z (" *
               "POINT Z (2 5 8)," *
               "POLYGON Z (" *
diff --git a/test/test_geos_operations.jl b/test/test_geos_operations.jl
index bcc2e4af..2ad6f9e6 100644
--- a/test/test_geos_operations.jl
+++ b/test/test_geos_operations.jl
@@ -1,6 +1,9 @@
 using Test
 import ArchGDAL as AG
+import GeoInterface
 using Extents
+const GI = GeoInterface
+
 
 @testset "test_geos_operations.jl" begin
     function equivalent_to_wkt(geom::AG.Geometry, wkt::String)
@@ -25,7 +28,8 @@ using Extents
                 AG.pointalongline(ls, dist) do pt1
                     AG.createpoint(dest) do pt2
                         @test AG.toWKT(pt1) == AG.toWKT(pt2)
-                    end
+                        @test AG.toWKT(GI.convert(GI, pt1)) == AG.toWKT(GI.convert(GI, pt2))
+                    end                
                 end
                 @test AG.toWKT(AG.pointalongline(ls, dist)) ==
                       AG.toWKT(AG.createpoint(dest))
@@ -38,6 +42,9 @@ using Extents
             AG.fromWKT("POLYGON EMPTY") do g2
                 @test AG.contains(g1, g2) == false
                 @test AG.contains(g2, g1) == false
+                # Empty geometries dont work so well accross packages yet
+                @test_broken AG.contains(GI.convert(GI, g1), GI.convert(GI, g2)) == false
+                @test_broken AG.contains(GI.convert(GI, g2), GI.convert(GI, g1)) == false
             end
         end
 
@@ -45,6 +52,8 @@ using Extents
             AG.fromWKT("POINT(2 2)") do g2
                 @test AG.contains(g1, g2) == true
                 @test AG.contains(g2, g1) == false
+                @test AG.contains(GI.convert(GI, g1), g2) == true
+                @test AG.contains(GI.convert(GI, g2), g1) == false
             end
         end
 
@@ -52,6 +61,8 @@ using Extents
             AG.fromWKT("POLYGON((1 1,1 2,2 2,2 1,1 1))") do g2
                 @test AG.contains(g1, g2) == true
                 @test AG.contains(g2, g1) == false
+                @test AG.contains(g1, GI.convert(GI, g2)) == true
+                @test AG.contains(g2, GI.convert(GI, g1)) == false
             end
         end
     end
@@ -65,7 +76,9 @@ using Extents
                     @test AG.isempty(output) == false
                     @test AG.toWKT(output) == AG.toWKT(expected)
                 end
-                @test AG.toWKT(AG.convexhull(input)) == AG.toWKT(expected)
+                @test AG.toWKT(AG.convexhull(input)) == 
+                      AG.toWKT(AG.convexhull(GI.convert(GI, input))) == 
+                      AG.toWKT(expected)
             end
         end
     end
@@ -77,29 +90,36 @@ using Extents
                 @test AG.isempty(g2) == true
                 @test AG.toWKT(g2) == "MULTILINESTRING EMPTY"
             end
-            @test AG.toWKT(AG.delaunaytriangulation(g1, 0, true)) ==
-                  "MULTILINESTRING EMPTY"
+            @test AG.toWKT(AG.delaunaytriangulation(g1, 0, true)) == "MULTILINESTRING EMPTY"
+            # Empty geometries dont work so well accross packages yet
+            @test_broken AG.toWKT(AG.delaunaytriangulation(GI.convert(GI, g1), 0, true)) == "MULTILINESTRING EMPTY"
         end
         AG.fromWKT("POINT(0 0)") do g1
             AG.delaunaytriangulation(g1, 0, false) do g2
-                @test AG.isempty(g2) == true
-                @test AG.toWKT(g2) == "GEOMETRYCOLLECTION EMPTY"
+                @test AG.isempty(g2) == AG.isempty(GI.convert(GI, g2)) == true
+                @test AG.toWKT(g2) == AG.toWKT(GI.convert(GI, g2)) == "GEOMETRYCOLLECTION EMPTY"
             end
             @test AG.toWKT(AG.delaunaytriangulation(g1, 0, false)) ==
+                  AG.toWKT(AG.delaunaytriangulation(GI.convert(GI, g1), 0, false)) ==
                   "GEOMETRYCOLLECTION EMPTY"
         end
         AG.fromWKT("MULTIPOINT(0 0, 5 0, 10 0)") do g1
             AG.delaunaytriangulation(g1, 0, false) do g2
-                @test AG.toWKT(g2) == "GEOMETRYCOLLECTION EMPTY"
+                @test AG.toWKT(g2) == 
+                      AG.toWKT(GI.convert(GI, g2)) == 
+                        "GEOMETRYCOLLECTION EMPTY"
             end
             AG.delaunaytriangulation(g1, 0, true) do g2
-                @test AG.toWKT(g2) == "MULTILINESTRING ((5 0,10 0),(0 0,5 0))"
+                @test AG.toWKT(g2) == 
+                      AG.toWKT(GI.convert(GI, g2)) == 
+                "MULTILINESTRING ((5 0,10 0),(0 0,5 0))"
             end
         end
         AG.fromWKT("MULTIPOINT(0 0, 10 0, 10 10, 11 10)") do g1
             AG.delaunaytriangulation(g1, 2.0, true) do g2
-                @test AG.toWKT(g2) ==
-                      "MULTILINESTRING ((0 0,10 10),(0 0,10 0),(10 0,10 10))"
+                @test AG.toWKT(g2) == 
+                      AG.toWKT(GI.convert(GI, g2)) == 
+                          "MULTILINESTRING ((0 0,10 10),(0 0,10 0),(10 0,10 10))"
             end
         end
     end
@@ -108,6 +128,7 @@ using Extents
         AG.fromWKT("POINT(10 10)") do g1
             AG.fromWKT("POINT(3 6)") do g2
                 @test AG.distance(g1, g2) ≈ 8.06225774829855 atol = 1e-12
+                @test AG.distance(GI.convert(GI, g1), GI.convert(GI, g2)) ≈ 8.06225774829855 atol = 1e-12
             end
         end
     end
@@ -122,6 +143,10 @@ using Extents
                 @test AG.toWKT(result) == wkt2
             end
             @test AG.toWKT(f(geom)) == wkt2
+            if parentmodule(f) != GeoInterface && GI.ngeom(geom) != 0
+                wrapped_geom = GI.convert(GI, geom)
+                @test AG.toWKT(f(wrapped_geom)) == wkt2
+            end
         end
     end
 
@@ -164,6 +189,11 @@ using Extents
                     @test AG.toWKT(result) == wkt3
                 end
                 @test AG.toWKT(f(geom1, geom2)) == wkt3
+                if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
+                    wrapped_geom1 = GI.convert(GI, geom1)
+                    wrapped_geom2 = GI.convert(GI, geom2)
+                    @test AG.toWKT(f(wrapped_geom1, wrapped_geom2)) == wkt3
+                end
             end
         end
     end
@@ -176,91 +206,145 @@ using Extents
         AG.fromWKT(wkt1) do geom1
             AG.fromWKT(wkt2) do geom2
                 @test AG.toWKT(f(geom1, geom2)) == wkt3
+                if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
+                    wrapped_geom1 = GI.convert(GI, geom1)
+                    wrapped_geom2 = GI.convert(GI, geom2)
+                    @test AG.toWKT(f(wrapped_geom1, wrapped_geom2)) == wkt3
+                end
             end
         end
     end
+    function test_method_simple(f1::Function, f2::Function, args...)
+        test_method_simple(f1, args...)
+        test_method_simple(f2, args...)
+    end
 
     function test_predicate(f::Function, wkt1, wkt2, result::Bool)
         AG.fromWKT(wkt1) do geom1
             AG.fromWKT(wkt2) do geom2
+                # Test GDAL geoms
                 @test f(geom1, geom2) == result
+                # Test GeoInterface geoms
+                if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
+                    wrapped_geom1 = GI.convert(GI, geom1)
+                    wrapped_geom2 = GI.convert(GI, geom2)
+                    @test f(wrapped_geom1, wrapped_geom2) == result
+                end
             end
         end
     end
+    function test_predicate(f1::Function, f2::Function, wkt1, wkt2, result::Bool)
+        test_predicate(f1, wkt1, wkt2, result::Bool)
+        test_predicate(f2, wkt1, wkt2, result::Bool)
+    end
 
-    @testset "GeoInterface" begin
+    @testset "Predicates" begin
         test_predicate(
-            GI.intersects,
+            GI.intersects, AG.intersects,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             true,
         )
         test_predicate(
-            GI.equals,
+            GI.equals, AG.equals, 
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             false,
         )
         test_predicate(
-            GI.disjoint,
+            GI.disjoint, AG.disjoint,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             false,
         )
         test_predicate(
-            GI.touches,
+            GI.touches, AG.touches,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(1 1)",
             true,
         )
         test_predicate(
-            GI.crosses,
+            GI.crosses, AG.crosses,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             false,
         )
         test_predicate(
-            GI.within,
+            GI.within, AG.within,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             false,
         )
         test_predicate(
-            GI.contains,
+            GI.contains, AG.contains,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             true,
         )
         test_predicate(
-            GI.overlaps,
+            GI.overlaps, AG.overlaps,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             false,
         )
+    end
+
+    @testset "Simple methods" begin
         test_method_simple(
-            GI.union,
+            GI.union, AG.union,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             "POLYGON ((1 5,5 5,5 1,1 1,1 5))",
         )
         test_method_simple(
-            GI.intersection,
+            GI.intersection, AG.intersection,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             "POINT (2 2)",
         )
         test_method_simple(
-            GI.difference,
+            GI.difference, AG.difference,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             "POLYGON ((1 5,5 5,5 1,1 1,1 5))",
         )
         test_method_simple(
-            GI.symdifference,
+            GI.symdifference, AG.symdifference,
             "POLYGON((1 1,1 5,5 5,5 1,1 1))",
             "POINT(2 2)",
             "POLYGON ((1 5,5 5,5 1,1 1,1 5))",
         )
+    end
+
+    @testset "ArchGDAL methods" begin
+        @test AG.distance(
+            AG.fromWKT("POLYGON((1 1,1 5,5 5,5 1,1 1))"),
+            AG.fromWKT("POINT(2 2)"),
+        ) == 0.0
+
+        # This could be `length` to match the GeoInterface method
+        @test AG.geomlength(AG.fromWKT("LINESTRING(0 0, 10 0)")) == 10
+
+        @test AG.area(AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")) == 100
+
+        @test GI.coordinates(
+            AG.buffer(AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"), 1),
+        )[1][1] == [-1.0, 0.0]
+
+        @test AG.toWKT(
+            AG.convexhull(AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")),
+        ) == "POLYGON ((0 0,0 10,10 10,10 0,0 0))"
+
+        # @test AG.asbinary(
+        #     AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"),
+        # )[1:10] ==
+        #       UInt8[0x01, 0x03, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x05]
+
+        # @test AG.astext(AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")) ==
+        #       "POLYGON ((0 0,10 0,10 10,0 10,0 0))"
+    end
+
+    @testset "GeoInterface methods" begin
         @test GI.distance(
             AG.fromWKT("POLYGON((1 1,1 5,5 5,5 1,1 1))"),
             AG.fromWKT("POINT(2 2)"),
@@ -398,6 +482,7 @@ using Extents
                 @test AG.toWKT(g2) == "POLYGON EMPTY"
             end
             @test AG.toWKT(AG.simplify(g1, 0.0)) == "POLYGON EMPTY"
+            @test AG.toWKT(AG.simplify(GI.convert(GI, g1), 0.0)) == "POLYGON EMPTY"
         end
 
         AG.fromWKT(
@@ -408,6 +493,7 @@ using Extents
                       "POLYGON ((56.5286666667 25.2101666667,56.529 25.2105,56.5288333333 25.2103333333,56.5286666667 25.2101666667))"
             end
             @test AG.toWKT(AG.simplifypreservetopology(g1, 43.2)) ==
+                  AG.toWKT(AG.simplifypreservetopology(GI.convert(GI, g1), 43.2)) ==
                   "POLYGON ((56.5286666667 25.2101666667,56.529 25.2105,56.5288333333 25.2103333333,56.5286666667 25.2101666667))"
         end
     end