Skip to content

Commit

Permalink
Multipolygon Clipping (#102)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@skygering
skygering authored Apr 11, 2024
1 parent 45315ee commit f5f9b0a
Show file tree
Hide file tree
Showing 8 changed files with 338 additions and 18 deletions.
1 change: 1 addition & 0 deletions src/GeometryOps.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
76 changes: 71 additions & 5 deletions src/methods/clipping/difference.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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},
Expand Down
76 changes: 69 additions & 7 deletions src/methods/clipping/intersection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand Down
82 changes: 76 additions & 6 deletions src/methods/clipping/union.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit f5f9b0a

Please sign in to comment.