Multipolygon Clipping (#102)

* Add minimal multipolygon correction

* Multipolygon clipping working

* Update multipolygon clipping documentation

* Explain fix_multipoly

* Make fix_multipoly more general

* Fix multipolygon type stability

* Remove GO ref

* Update fix_multipoly default comment

Author: skygering

src/GeometryOps.jl

@@ -54,6 +54,7 @@ include("transformations/tuples.jl")
 include("transformations/transform.jl")
 include("transformations/correction/geometry_correction.jl")
 include("transformations/correction/closed_ring.jl")
+include("transformations/correction/intersecting_polygons.jl")
 
 # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)

src/methods/clipping/difference.jl

@@ -3,13 +3,19 @@ export difference
 
 
 """
-    difference(geom_a, geom_b, [T::Type]; target::Type)
+    difference(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 Return the difference between two geometries as a list of geometries. Return an empty list
 if none are found. The type of the list will be constrained as much as possible given the
 input geometries. Furthermore, the user can provide a `taget` type as a keyword argument and
 a list of target geometries found in the difference will be returned. The user can also
-provide a float type that they would like the points of returned geometries to be. 
+provide a float type that they would like the points of returned geometries to be. If the
+user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to false if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
 
 ## Example 
 
@@ -27,9 +33,9 @@ GI.coordinates.(diff_poly)
 ```
 """
 function difference(
-    geom_a, geom_b, ::Type{T} = Float64; target=nothing,
+    geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
-    return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 #= The 'difference' function returns the difference of two polygons as a list of polygons.
@@ -38,7 +44,8 @@ polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.27
 function _difference(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...
 ) where T
     # Get the exterior of the polygons
     ext_a = GI.getexterior(poly_a)
@@ -99,6 +106,65 @@ tracing a_list, else step backwards, where x is the entry/exit status and y is a
 that is true if we are on a_list and false if we are on b_list. =#
 _diff_step(x, y) = (x ⊻ y) ? 1 : (-1)
 
+#= Polygon with multipolygon difference - note that all intersection regions between
+`poly_a` and any of the sub-polygons of `multipoly_b` are removed from `poly_a`. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    kwargs...,
+) where T
+    polys = [tuples(poly_a, T)]
+    for poly_b in GI.getpolygon(multipoly_b)
+        isempty(polys) && break
+        polys = mapreduce(p -> difference(p, poly_b; target), append!, polys)
+    end
+    return polys
+end
+
+#= Multipolygon with polygon difference - note that all intersection regions between
+sub-polygons of `multipoly_a` and `poly_b` will be removed from the corresponding
+sub-polygon. Unless specified with `fix_multipoly = nothing`, `multipolygon_a` will be
+validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
+        multipoly_a = fix_multipoly(multipoly_a)
+    end
+    polys = Vector{_get_poly_type(T)}()
+    sizehint!(polys, GI.npolygon(multipoly_a))
+    for poly_a in GI.getpolygon(multipoly_a)
+        append!(polys, difference(poly_a, poly_b; target))
+    end
+    return polys
+end
+
+#= Multipolygon with multipolygon difference - note that all intersection regions between
+sub-polygons of `multipoly_a` and sub-polygons of `multipoly_b` will be removed from the
+corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly = nothing`,
+`multipolygon_a` will be validated using the given (defauly is `UnionIntersectingPolygons()`)
+correction. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
+        multipoly_a = fix_multipoly(multipoly_a)
+        fix_multipoly = nothing
+    end
+    local polys
+    for (i, poly_b) in enumerate(GI.getpolygon(multipoly_b))
+        polys = difference(i == 1 ? multipoly_a : GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
+    end
+    return polys
+end
+
 # Many type and target combos aren't implemented
 function _difference(
     ::TraitTarget{Target}, ::Type{T},

src/methods/clipping/intersection.jl

@@ -10,13 +10,19 @@ Enum for the orientation of a line with respect to a curve. A line can be
 @enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4
 
 """
-    intersection(geom_a, geom_b, [T::Type]; target::Type)
+    intersection(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 Return the intersection between two geometries as a list of geometries. Return an empty list
 if none are found. The type of the list will be constrained as much as possible given the
 input geometries. Furthermore, the user can provide a `target` type as a keyword argument and
 a list of target geometries found in the intersection will be returned. The user can also
-provide a float type that they would like the points of returned geometries to be. 
+provide a float type that they would like the points of returned geometries to be. If the
+user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to nothing if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
 
 ## Example
 
@@ -34,16 +40,17 @@ GI.coordinates.(inter_points)
 ```
 """
 function intersection(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing,
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
-    return _intersection(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    return _intersection(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 # Curve-Curve Intersections with target Point
 _intersection(
     ::TraitTarget{GI.PointTrait}, ::Type{T},
     trait_a::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_a,
-    trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b,
+    trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b;
+    kwargs...,
 ) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b)
 
 #= Polygon-Polygon Intersections with target Polygon
@@ -53,7 +60,8 @@ DOI: https://doi.org/10.1145/274363.274364 =#
 function _intersection(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
 ) where {T}
     # First we get the exteriors of 'poly_a' and 'poly_b'
     ext_a = GI.getexterior(poly_a)
@@ -98,11 +106,65 @@ _inter_delay_bounce_f(x, _) = x
 point, else step backwards where x is the entry/exit status. =#
 _inter_step(x, _) =  x ? 1 : (-1)
 
+#= Polygon with multipolygon intersection - note that all intersection regions between
+`poly_a` and any of the sub-polygons of `multipoly_b` are counted as intersection polygons.
+Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be validated using
+the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_b to prevent duplicated intersection regions
+        multipoly_b = fix_multipoly(multipoly_b)
+    end
+    polys = Vector{_get_poly_type(T)}()
+    for poly_b in GI.getpolygon(multipoly_b)
+        append!(polys, intersection(poly_a, poly_b; target))
+    end
+    return polys
+end
+
+#= Multipolygon with polygon intersection is equivalent to taking the intersection of the
+poylgon with the multipolygon and thus simply switches the order of operations and calls the
+above method. =#
+_intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T = intersection(poly_b, multipoly_a; target , kwargs...)
+
+#= Multipolygon with multipolygon intersection - note that all intersection regions between
+any sub-polygons of `multipoly_a` and any of the sub-polygons of `multipoly_b` are counted
+as intersection polygons. Unless specified with `fix_multipoly = nothing`, both 
+`multipolygon_a` and `multipolygon_b` will be validated using the given (default is
+`UnionIntersectingPolygons()`) correction. =#
+function _intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix both multipolygons to prevent duplicated regions
+        multipoly_a = fix_multipoly(multipoly_a)
+        multipoly_b = fix_multipoly(multipoly_b)
+        fix_multipoly = nothing
+    end
+    polys = Vector{_get_poly_type(T)}()
+    for poly_a in GI.getpolygon(multipoly_a)
+        append!(polys, intersection(poly_a, multipoly_b; target, fix_multipoly))
+    end
+    return polys
+end
+
 # Many type and target combos aren't implemented
 function _intersection(
     ::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
-    trait_b::GI.AbstractTrait, geom_b,
+    trait_b::GI.AbstractTrait, geom_b;
+    kwargs...,
 ) where {Target, T}
     @assert(
         false,

src/methods/clipping/union.jl

@@ -2,13 +2,19 @@
 export union
 
 """
-    union(geom_a, geom_b, [::Type{T}]; target::Type)
+    union(geom_a, geom_b, [::Type{T}]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 Return the union between two geometries as a list of geometries. Return an empty list if
 none are found. The type of the list will be constrained as much as possible given the input
 geometries. Furthermore, the user can provide a `taget` type as a keyword argument and a
 list of target geometries found in the difference will be returned. The user can also
-provide a float type 'T' that they would like the points of returned geometries to be. 
+provide a float type 'T' that they would like the points of returned geometries to be. If
+the user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to false if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
     
 Calculates the union between two polygons.
 ## Example
@@ -27,9 +33,9 @@ GI.coordinates.(union_poly)
 ```
 """
 function union(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing,
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
 ) where {T<:AbstractFloat}
-    _union(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    _union(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 #= This 'union' implementation returns the union of two polygons. The algorithm to determine
@@ -38,7 +44,8 @@ Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
 function _union(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
 ) where T
     # First, I get the exteriors of the two polygons
     ext_a = GI.getexterior(poly_a)
@@ -196,12 +203,75 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
     return
 end
 
+#= Polygon with multipolygon union - note that all sub-polygons of `multipoly_b` will be
+included, unioning these sub-polygons with `poly_a` where they intersect. Unless specified
+with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default
+is `UnionIntersectingPolygons()`) correction. =#
+function _union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_b to prevent repeated regions in the output
+        multipoly_b = fix_multipoly(multipoly_b)
+    end
+    polys = [tuples(poly_a, T)]
+    for poly_b in GI.getpolygon(multipoly_b)
+        if intersects(polys[1], poly_b)
+            # If polygons intersect and form a new polygon, swap out polygon
+            new_polys = union(polys[1], poly_b; target)
+            if length(new_polys) > 1 # case where they intersect by just one point
+                push!(polys, tuples(poly_b, T))  # add poly_b to list
+            else
+                polys[1] = new_polys[1]
+            end
+        else
+            # If they don't intersect, poly_b is now a part of the union as its own polygon
+            push!(polys, tuples(poly_b, T))
+        end
+    end
+    return polys
+end
+
+#= Multipolygon with polygon union is equivalent to taking the union of the poylgon with the
+multipolygon and thus simply switches the order of operations and calls the above method. =#
+_union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T = union(poly_b, multipoly_a; target, kwargs...)
+
+#= Multipolygon with multipolygon union - note that all of the sub-polygons of `multipoly_a`
+and the sub-polygons of `multipoly_b` are included and combined together where there are
+intersections. Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be
+validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_b to prevent repeated regions in the output
+        multipoly_b = fix_multipoly(multipoly_b)
+        fix_multipoly = nothing
+    end
+    multipolys = multipoly_b
+    local polys
+    for poly_a in GI.getpolygon(multipoly_a)
+        polys = union(poly_a, multipolys; target, fix_multipoly)
+        multipolys = GI.MultiPolygon(polys)
+    end
+    return polys
+end
 
 # Many type and target combos aren't implemented
 function _union(
     ::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
-    trait_b::GI.AbstractTrait, geom_b,
+    trait_b::GI.AbstractTrait, geom_b;
+    kwargs...
 ) where {Target,T}
     throw(ArgumentError("Union between $trait_a and $trait_b with target $Target isn't implemented yet."))
     return nothing