Union between more than two polygons

[asinghvi17]

How would this be implemented? This is a pretty big blocker to seamlessly merging many polygons into, for example, a multipolygon.

This would also enable unions between multipolygons and polygons.

LibGEOS does this pretty seamlessly, so we probably should as well :D

[rafaqz]

reduce(union, polygons) ?

[asinghvi17]

It returns a list of polygons though - so you run into the same issue

[rafaqz]

Its also probably not optimal... like we probably want to decompose all polygons and calculated union simultaneously?

Im not accross the algorithms here but it seems like there is a chance to sort all the edges by Y axis in a single list and recombine them in one go.

[asinghvi17]

This can be done by reduce and MultiPolygonTarget now...

[tiemvanderdeure]

How can this be done? I am trying to figure it out but the documentation for union doesn't really tell me what I want to know.

Here is what I have been trying, everything errors except for the LibGEOS.union

using NaturalEarth, LibGEOS
import GeometryOps as GO, GeoInterface as GI

countries = naturalearth("ne_110m_admin_0_countries").geometry
world = reduce(LibGEOS.union, countries)
reduce((x,y) -> GO.union(x, y; target = GO.MultiPolygonTrait()), countries)
reduce((x,y) -> GO.union(x, y; target = GO.PolygonTrait()), countries)
reduce((x,y) -> GO.union(x, y; target = GO.PolygonTrait())[1], countries)

GI.MultiPolygon(countries)

Even converting to a MultiPolygon doesn't work because some countries are MultiPolygons to start with whereas others are Polygons.

If this is already possible we should probably take another look at the user-facing side of union.

[asinghvi17]

Ah I guess we just need to implement multipolygon union with the multipolygon trait, it looks like we can do it so I think that just got left behind.

This monkeypatch:

function _union(
    ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
    trait_a::GI.MultiPolygonTrait, geom_a, 
    trait_b::GI.PolygonTrait, geom_b;
    kwargs...
) where {T}
    return GI.MultiPolygon(
        _union(
            TraitTarget{GI.PolygonTrait}(), T, 
            trait_a, geom_a, 
            trait_b, geom_b; 
            kwargs...
        )
    )
end

at least gives a better error...

julia> reduce((x,y) -> GO.union(x, y; target = GO.MultiPolygonTrait()), countries)
ERROR: MethodError: Cannot `convert` an object of type
  GeoInterface.Wrappers.Polygon{false,false,Array{GeoInterface.Wrappers.LinearRing{false,false,Array{Tuple{Float64,Float64},1},Nothing,Nothing},1},Nothing,Nothing} to an object of type
  GeoInterface.Wrappers.Polygon{false,false,Array{GeoInterface.Wrappers.LinearRing{false,false,Array{Tuple{Float64,Float64},1},Nothing,EPSG{1}},1},Nothing,EPSG{1}}

(btw, this version reduce function should be the way to go)

[asinghvi17]

If you want easier syntax you can always do

reduce((x,y) -> GO.union(GO.GEOS(), x, y; target = GO.MultiPolygonTrait()), countries)

and once this is fixed just cut out the GEOS()...