GeometryOps.jl

Appearance

Intersecting Polygons ​

julia
export UnionIntersectingPolygons

Multipolygon components must have disjoint interiors and may meet only at isolated points. UnionIntersectingPolygons merges overlaps while preserving the covered area.

Example ​

This multipolygon is invalid because it repeats the same polygon:

julia
import GeoInterface as GI
polygon = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]])
multipolygon = GI.MultiPolygon([polygon, polygon])
GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(0.0, 0.0), … (3) … , (0.0, 0.0)])]), GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(0.0, 0.0), … (3) … , (0.0, 0.0)])])])

Apply the correction:

julia
import GeometryOps as GO
GO.fix(multipolygon, corrections = [GO.UnionIntersectingPolygons()])
GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(0.0, 0.0), … (3) … , (0.0, 0.0)])])])

The corrected multipolygon contains one component.

Implementation ​

julia
"""
    UnionIntersectingPolygons([manifold = Planar()], [T = Float64]) <: GeometryCorrection
    UnionIntersectingPolygons(algorithm::FosterHormannClipping, [T = Float64])

Merge intersecting components with union, preserving the covered area. Result components are
disjoint except for possible point contacts.

Default clipping corrections inherit the caller's manifold, algorithm, and numeric type. An
explicitly supplied correction retains its own settings.

See also `GeometryCorrection`.
"""
struct UnionIntersectingPolygons{A <: FosterHormannClipping, T <: AbstractFloat} <: GeometryCorrection
    algorithm::A
end
UnionIntersectingPolygons(alg::FosterHormannClipping, ::Type{T} = Float64) where {T <: AbstractFloat} =
    UnionIntersectingPolygons{typeof(alg), T}(alg)
UnionIntersectingPolygons(m::Manifold = Planar(), ::Type{T} = Float64) where {T <: AbstractFloat} =
    UnionIntersectingPolygons(FosterHormannClipping(m), T)

application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait

function (correction::UnionIntersectingPolygons{A, T})(::GI.MultiPolygonTrait, multipoly) where {A, T}
    alg = correction.algorithm
    P = _fh_out_point_type(alg.manifold, multipoly, T)
    union_multipoly = _fh_multipolygon(_get_poly_type(T, P)[_fh_as_poly(P, poly, T) for poly in GI.getpolygon(multipoly)]; crs = GI.crs(multipoly))
    n_polys = GI.npolygon(multipoly)
    if n_polys > 1
        keep_idx = trues(n_polys)  # keep track of sub-polygons to remove

Combine any sub-polygons that intersect

julia
        for (curr_idx, _) in Iterators.filter(last, Iterators.enumerate(keep_idx))
            curr_poly = union_multipoly.geom[curr_idx]
            poly_disjoint = false
            while !poly_disjoint
                poly_disjoint = true  # assume current polygon is disjoint from others
                for (next_idx, _) in Iterators.filter(last, Iterators.drop(Iterators.enumerate(keep_idx), curr_idx))
                    next_poly = union_multipoly.geom[next_idx]
                    if intersects(alg.manifold, curr_poly, next_poly)  # if two polygons intersect
                        new_polys = union(alg, curr_poly, next_poly, T; target = GI.PolygonTrait())
                        n_new_polys = length(new_polys)
                        if n_new_polys == 1  # if polygons combined
                            poly_disjoint = false
                            union_multipoly.geom[curr_idx] = new_polys[1]
                            curr_poly = union_multipoly.geom[curr_idx]
                            keep_idx[next_idx] = false
                        end
                    end
                end
            end
        end
        keepat!(union_multipoly.geom, keep_idx)
    end
    return union_multipoly
end

"""
    DiffIntersectingPolygons([manifold = Planar()], [T = Float64]) <: GeometryCorrection
    DiffIntersectingPolygons(algorithm::FosterHormannClipping, [T = Float64])

Remove component overlaps with `difference`, preserving the covered area. Result
components are disjoint except for possible point contacts.

See also `GeometryCorrection`, `UnionIntersectingPolygons`.
"""
struct DiffIntersectingPolygons{A <: FosterHormannClipping, T <: AbstractFloat} <: GeometryCorrection
    algorithm::A
end
DiffIntersectingPolygons(alg::FosterHormannClipping, ::Type{T} = Float64) where {T <: AbstractFloat} =
    DiffIntersectingPolygons{typeof(alg), T}(alg)
DiffIntersectingPolygons(m::Manifold = Planar(), ::Type{T} = Float64) where {T <: AbstractFloat} =
    DiffIntersectingPolygons(FosterHormannClipping(m), T)

application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait

function (correction::DiffIntersectingPolygons{A, T})(::GI.MultiPolygonTrait, multipoly) where {A, T}
    alg = correction.algorithm
    P = _fh_out_point_type(alg.manifold, multipoly, T)
    diff_multipoly = _fh_multipolygon(_get_poly_type(T, P)[_fh_as_poly(P, poly, T) for poly in GI.getpolygon(multipoly)]; crs = GI.crs(multipoly))
    n_starting_polys = GI.npolygon(multipoly)
    n_polys = n_starting_polys
    if n_polys > 1
        keep_idx = trues(n_polys)  # keep track of sub-polygons to remove

Break apart any sub-polygons that intersect

julia
        for curr_idx in 1:n_starting_polys
            !keep_idx[curr_idx] && continue
            for next_idx in (curr_idx + 1):n_starting_polys
                !keep_idx[next_idx] && continue
                next_poly = diff_multipoly.geom[next_idx]
                n_new_polys = 0
                curr_pieces_added = (n_polys + 1):(n_polys + n_new_polys)
                for curr_piece_idx in Iterators.flatten((curr_idx:curr_idx, curr_pieces_added))
                    !keep_idx[curr_piece_idx] && continue
                    curr_poly = diff_multipoly.geom[curr_piece_idx]
                    if intersects(alg.manifold, curr_poly, next_poly)  # if two polygons intersect
                        new_polys = difference(alg, curr_poly, next_poly, T; target = GI.PolygonTrait())
                        n_new_pieces = length(new_polys) - 1
                        if n_new_pieces < 0  # current polygon is covered by next_polygon
                            keep_idx[curr_piece_idx] = false
                            break
                        elseif n_new_pieces ≥ 0
                            diff_multipoly.geom[curr_piece_idx] = new_polys[1]
                            curr_poly = diff_multipoly.geom[curr_piece_idx]
                            if n_new_pieces > 0 # current polygon breaks into several pieces
                                append!(diff_multipoly.geom, @view new_polys[2:end])
                                append!(keep_idx, trues(n_new_pieces))
                                n_new_polys += n_new_pieces
                            end
                        end
                    end
                end
                n_polys += n_new_polys
            end
        end
        keepat!(diff_multipoly.geom, keep_idx)
    end
    return diff_multipoly
end

This page was generated using Literate.jl.