Wrap all GEOS() results in GI wrappers (#216)

* Make all GEOS methods wrap the result in a GI wrapper

to preserve CRS and potentially calculate extent.

* Wrap the result of convexhull

* Convert directly to LibGEOS points

Co-authored-by: Rafael Schouten <[email protected]>

* less nesting

* Document `_wrap` and why we use it

---------

Co-authored-by: Rafael Schouten <[email protected]>

ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl

@@ -15,10 +15,23 @@ for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(Ge
     @eval using GeometryOps: $name
 end
 
+"""
+    _wrap(geom; crs, calc_extent)
+
+Wraps `geom` in a GI wrapper geometry of its geometry trait.  This allows us
+to attach CRS and extent info to geometry types which otherwise could not hold
+those, like LibGEOS and WKB geometries.
+
+Returns a GI wrapper geometry, for which `parent(result) == geom`.
+"""
+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")
 
 include("simple_overrides.jl")
 
-end
+end

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
@@ -66,15 +66,20 @@ end
 
 # ## Convex hull
 function GO.convex_hull(::GEOS, geoms)
-    return LG.convexhull(
+    chull = LG.convexhull(
         LG.MultiPoint(
             collect(
                 GO.flatten(
-                    x -> GI.convert(LG, x), 
+                    x -> GI.convert(LG.Point, x), 
                     GI.PointTrait, 
                     geoms
                 )
             )
         )
+    ); 
+    return _wrap(
+        chull;
+        crs = GI.crs(geoms), 
+        calc_extent = false
     )
-end
+end