Union and Intersection Degenerate Case

[skygering]

Sorry that these coordinates are complicated, I haven't been able to recreate with simpler coordinates.

using LibGEOS, Plots
c1 = [[[0.0, 65905.78568220709], [12540.144108116785, 66644.10887217366], [13639.90993528687, 64693.73062699103], [13323.160413385054, 61494.1194419435], [2375.0223287135154, 53673.4087205281], [0.0, 52759.58192468053], [0.0, 65905.78568220709]]]
c2 = [[[23303.577415035626, 55323.60484150198], [19851.52938808218, 49637.68131132904], [2375.022328713516, 53673.4087205281], [13323.160413385054, 61494.1194419435], [23303.577415035626, 55323.60484150198]]]
p1 = LibGEOS.Polygon(c1)
p2 = LibGEOS.Polygon(c2)
plot(p1) # top in image below
plot!(p2) # bottom in image below

When plotted, they look like follows:

They don't appear to be overlapping. However, the union is the lower of the two polygons, and the intersection is the bottom.

u = LibGEOS.union(p1, p2)
plot!(u)

i = LibGEOS.intersection(p1, p2)
plot!(i)

All polygons are valid if tested with ifValid. This doesn't happen if I set up two squares that share a border. Then the intersection is the line between them, and the union is a larger rectangle made of both.

Please let me know if I am missing something.

[jaakkor2]

Seems to work with p2 = LibGEOS.Polygon([reverse(c2[1])]). Therefore looks like https://github.com/libgeos/geos/issues/827 .

[skygering]

@jaakkor2 I am not convinced it is that problem. Both sets of coordinates above are oriented clockwise, so their winding order is the same in my original question, while the GEOS issue relates to opposite winding orders. This can be seen as follows:

cs1 = LibGEOS.createCoordSeq(c1[1])
cs2 = cs2 = LibGEOS.createCoordSeq(c2[1])
println(LibGEOS.isCCW(cs1))
println(LibGEOS.isCCW(cs2))

which will return true. I then decided to try flipping both sets of coordinates winding orders to see what would happen if they were both clockwise.

p1 = LibGEOS.Polygon([reverse(c1[1])])
p2 = LibGEOS.Polygon([reverse(c2[1])])
plot(p1) # top in image below
plot!(p2) # bottom in image below

u = LibGEOS.union(p1, p2)
plot!(u)

As you can see below, this actually creates a multipolygon!

Finally, we take the intersection:

i = LibGEOS.intersection(p1, p2)
plot!(i)

We get a point this time.

To me, these three totally different answers are suggesting that these functions don't support the degenerate case of two polygons being aligned directly along an edge.

This is supported by the function working when the polygons are slightly perturbed so that they don't line up directly. Clearly this is a problem with the algorithm used by the underlying code base, but it is worth noting here so that users are aware.

[jaakkor2]

Flipping winding of p2 gives the desired result. I agree the algorithm is faulty.

using LibGEOS, GLMakie, GeoInterfaceMakie
GeoInterfaceMakie.@enable(LibGEOS.AbstractGeometry)
c1 = [[[0.0, 65905.78568220709], [12540.144108116785, 66644.10887217366], [13639.90993528687, 64693.73062699103], [13323.160413385054, 61494.1194419435], [2375.0223287135154, 53673.4087205281], [0.0, 52759.58192468053], [0.0, 65905.78568220709]]]
c2 = [[[23303.577415035626, 55323.60484150198], [19851.52938808218, 49637.68131132904], [2375.022328713516, 53673.4087205281], [13323.160413385054, 61494.1194419435], [23303.577415035626, 55323.60484150198]]]
p1 = LibGEOS.Polygon(c1)
#p2 = LibGEOS.Polygon(c2)
p2 = LibGEOS.Polygon([reverse(c2[1])]) # flip winding of p2
union_p1p2 = LibGEOS.union(p1,p2)
isect_p1p2 = intersection(p1,p2)
area(p1), area(p2), area(isect_p1p2), area(union_p1p2)
pl = plot(union_p1p2, color = :green, axis = (aspect = DataAspect(),))
plot!(isect_p1p2, color = :red)
save("polys.png", pl.figure)
display(pl.figure)

gives

julia> union_p1p2 = LibGEOS.union(p1,p2)
POLYGON ((12540.144108116785 66644.10887217366, 13639.90993528687 64693.73062699103, 13323.160413385054 61494.1194419435, 23303.577415035626 55323.60484150198, 19851.52938808218 49637.68131132904, 2375.0223287135154 53673.4087205281, 0 52759.58192468053, 0 65905.78568220709, 12540.144108116785 66644.10887217366))

julia> isect_p1p2 = intersection(p1,p2)
LINESTRING (13323.160413385054 61494.1194419435, 2375.0223287135154 53673.4087205281)

julia> area(p1), area(p2), area(isect_p1p2), area(union_p1p2)
(1.2650748762841602e8, 1.2945560385123302e8, 0.0, 2.55963091479649e8)