Make all GEOS methods wrap the result in a GI wrapper

to preserve CRS and potentially calculate extent.

ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl

@@ -15,6 +15,10 @@ for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(Ge
     @eval using GeometryOps: $name
 end
 
+function _wrap(geom; crs=GI.crs(geom), calc_extent = true)
+    return GI.geointerface_geomtype(GI.geomtrait(geom))(geom; crs, extent = GI.extent(geom, calc_extent))
+end
+
 include("buffer.jl")
 include("segmentize.jl")
 include("simplify.jl")

ext/GeometryOpsLibGEOSExt/buffer.jl

@@ -18,16 +18,17 @@ to_join_style(style::Symbol) = _GEOS_JOINSTYLE_LOOKUP[style]
 to_join_style(style::LG.GEOSBufJoinStyles) = style
 to_join_style(num::Integer) = num
 
-function GO.buffer(alg::GEOS, geometry, distance)
+function GO.buffer(alg::GEOS, geometry, distance; calc_extent = true, kwargs...)
     # The reason we use apply here is so that this also works with featurecollections,
     # tables, vectors of geometries, etc!
-    return apply(TraitTarget{GI.AbstractGeometryTrait}(), geometry) do geom
-        LG.bufferWithStyle(
+    return apply(TraitTarget{GI.AbstractGeometryTrait}(), geometry; kwargs...) do geom
+        newgeom = LG.bufferWithStyle(
             GI.convert(LG, geom), distance; 
             quadsegs = get(alg, :quadsegs, 8),
             endCapStyle = to_cap_style(get(alg, :endCapStyle, :round)),
             joinStyle = to_join_style(get(alg, :joinStyle, :round)),
             mitreLimit = get(alg, :mitreLimit, 5.0),
         )
+        return _wrap(newgeom; crs = GI.crs(geom), calc_extent)
     end
 end

ext/GeometryOpsLibGEOSExt/segmentize.jl

@@ -16,13 +16,20 @@ end
 
 _segmentize_geos(geom, max_distance) = _segmentize_geos(GI.convert(LG, geom), max_distance)
 
+function _wrap_and_segmentize_geos(geom, max_distance)
+    _wrap(_segmentize_geos(geom, max_distance); crs = GI.crs(geom), calc_extent = false)
+end
+
 # 2 behaviours:
 # - enforce: enforce the presence of a kwargs
 # - fetch: fetch the value of a kwargs, or return a default value
 @inline function GO.segmentize(alg::GEOS, geom; threaded::Union{Bool, GO.BoolsAsTypes} = _False())
     max_distance = enforce(alg, :max_distance, GO.segmentize)
     return GO.apply(
-        Base.Fix2(_segmentize_geos, max_distance), 
+        Base.Fix2(_wrap_and_segmentize_geos, max_distance), 
+        # TODO: should this just be a target on GI.AbstractGeometryTrait()?
+        # But Geos doesn't support eg RectangleTrait
+        # Maybe we need an abstract trait `GI.AbstractWKBGeomTrait`?
         GO.TraitTarget(GI.GeometryCollectionTrait(), GI.MultiPolygonTrait(), GI.PolygonTrait(), GI.MultiLineStringTrait(), GI.LineStringTrait(), GI.LinearRingTrait(), GI.MultiPointTrait(), GI.PointTrait()),
         geom; 
         threaded

ext/GeometryOpsLibGEOSExt/simple_overrides.jl

@@ -6,20 +6,20 @@ require conversion before calling.
 =#
 # ## Polygon set operations
 # ### Difference
-function GO.difference(::GEOS, geom_a, geom_b; target=nothing)
-    return LG.difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+function GO.difference(::GEOS, geom_a, geom_b; target=nothing, calc_extent = false)
+    return _wrap(LG.difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b)); crs = GI.crs(geom_a), calc_extent)
 end
 # ### Union
-function GO.union(::GEOS, geom_a, geom_b; target=nothing)
-    return LG.union(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+function GO.union(::GEOS, geom_a, geom_b; target=nothing, calc_extent = false)
+    return _wrap(LG.union(GI.convert(LG, geom_a), GI.convert(LG, geom_b)); crs = GI.crs(geom_a), calc_extent)
 end
 # ### Intersection
-function GO.intersection(::GEOS, geom_a, geom_b; target=nothing)
-    return LG.intersection(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+function GO.intersection(::GEOS, geom_a, geom_b; target=nothing, calc_extent = false)
+    return _wrap(LG.intersection(GI.convert(LG, geom_a), GI.convert(LG, geom_b)); crs = GI.crs(geom_a), calc_extent)
 end
 # ### Symmetric difference
-function GO.symdifference(::GEOS, geom_a, geom_b; target=nothing)
-    return LG.symmetric_difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+function GO.symdifference(::GEOS, geom_a, geom_b; target=nothing, calc_extent = false)
+    return _wrap(LG.symmetric_difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b)); crs = GI.crs(geom_a), calc_extent)
 end
 
 # ## DE-9IM boolean methods