Fix CRS propagation in union and difference clipping operations

Extend CRS handling from intersection to union and difference operations.
Use `mutual_crs` to preserve coordinate reference system information when
creating result polygons in non-crossing cases, hole handling, and
multipolygon operations. Fixes type mismatches when clipping geometries
with different CRS configurations.

Author: asinghvi17

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 AccurateArithmetic = "22286c92-06ac-501d-9306-4abd417d9753"
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
+ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -33,17 +34,17 @@ LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MakieThemes = "e296ed71-da82-5faf-88ab-0034a9761098"
+Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
-Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 
 [sources]
@@ -52,4 +53,4 @@ GeometryOpsCore = {path = "../GeometryOpsCore"}
 
 [compat]
 DocumenterVitepress = ">=0.2.1"
-GeoInterface = "1.6"
+GeoInterface = "1.6"

src/methods/clipping/difference.jl

@@ -57,19 +57,21 @@ function _difference(
     # Get the exterior of the polygons
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
+
+    crs = mutual_crs(poly_a, poly_b)
     # Find the difference of the exterior of the polygons
     a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact)
-    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b; crs)
     # if no crossing points, determine if either poly is inside of the other
     if isempty(polys)
         a_in_b, b_in_a = _find_non_cross_orientation(alg.manifold, a_list, b_list, ext_a, ext_b; exact)
         # add case for if they polygons are the same (all intersection points!)
         # add a find_first check to find first non-inter poly!
         if b_in_a && !a_in_b  # b in a and can't be the same polygon
-            poly_a_b_hole = GI.Polygon([tuples(ext_a), tuples(ext_b)])
+            poly_a_b_hole = GI.Polygon([tuples(ext_a; crs), tuples(ext_b; crs)]; crs)
             push!(polys, poly_a_b_hole)
         elseif !b_in_a && !a_in_b # polygons don't intersect
-            push!(polys, tuples(poly_a))
+            push!(polys, tuples(poly_a; crs))
             return polys
         end
     end
@@ -119,7 +121,8 @@ function _difference(
     ::GI.MultiPolygonTrait, multipoly_b;
     kwargs...,
 ) where T
-    polys = [tuples(poly_a, T)]
+    crs = mutual_crs(poly_a, multipoly_b)
+    polys = [tuples(poly_a, T; crs)]
     for poly_b in GI.getpolygon(multipoly_b)
         isempty(polys) && break
         polys = mapreduce(p -> difference(alg, p, poly_b; target), append!, polys)
@@ -140,7 +143,8 @@ function _difference(
     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)}()
+    crs = mutual_crs(multipoly_a, poly_b)
+    polys = Vector{_get_poly_type(T, crs)}()
     sizehint!(polys, GI.npolygon(multipoly_a))
     for poly_a in GI.getpolygon(multipoly_a)
         append!(polys, difference(alg, poly_a, poly_b; target))
@@ -163,6 +167,7 @@ function _difference(
         multipoly_a = fix_multipoly(multipoly_a)
         fix_multipoly = nothing
     end
+    crs = mutual_crs(multipoly_a, multipoly_b)
     local polys
     for (i, poly_b) in enumerate(GI.getpolygon(multipoly_b))
         #= Removing intersections of `multipoly_a`` with pieces of `multipoly_b`` - as
@@ -171,7 +176,7 @@ function _difference(
         polys = if i == 1
             difference(alg, multipoly_a, poly_b; target, fix_multipoly)
         else
-            difference(alg, GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
+            difference(alg, GI.MultiPolygon(polys; crs), poly_b; target, fix_multipoly)
         end
         #= One multipoly_a has been completely covered (and thus removed) there is no need to
         continue taking the difference =#

src/methods/clipping/union.jl

@@ -66,21 +66,23 @@ function _union(
     # First, I get the exteriors of the two polygons
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
+
+    crs = mutual_crs(poly_a, poly_b)
     # Then, I get the union of the exteriors
     a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact)
-    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b; crs)
     n_pieces = length(polys)
     # Check if one polygon totally within other and if so, return the larger polygon
     a_in_b, b_in_a = false, false
     if n_pieces == 0 # no crossing points, determine if either poly is inside the other
         a_in_b, b_in_a = _find_non_cross_orientation(alg, a_list, b_list, ext_a, ext_b; exact)
         if a_in_b
-            push!(polys, GI.Polygon([_linearring(tuples(ext_b))]))
+            push!(polys, GI.Polygon([_linearring(tuples(ext_b; crs))]; crs))
         elseif b_in_a
-            push!(polys,  GI.Polygon([_linearring(tuples(ext_a))]))
+            push!(polys,  GI.Polygon([_linearring(tuples(ext_a; crs))]; crs))
         else
-            push!(polys, tuples(poly_a))
-            push!(polys, tuples(poly_b))
+            push!(polys, tuples(poly_a; crs))
+            push!(polys, tuples(poly_b; crs))
             return polys
         end
     elseif n_pieces > 1
@@ -94,7 +96,7 @@ function _union(
     end
     # Add in holes
     if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
-        _add_union_holes!(alg, polys, a_in_b, b_in_a, poly_a, poly_b; exact)
+        _add_union_holes!(alg, polys, a_in_b, b_in_a, poly_a, poly_b; exact, crs)
     end
     # Remove unneeded collinear points on same edge
     _remove_collinear_points!(alg, polys, [false], poly_a, poly_b)
@@ -121,11 +123,11 @@ _union_step(x, _) = x ? (-1) : 1
 #= Add holes from two polygons to the exterior polygon formed by their union. If adding the
 the holes reveals that the polygons aren't actually intersecting, return the original
 polygons. =#
-function _add_union_holes!(alg::FosterHormannClipping, polys, a_in_b, b_in_a, poly_a, poly_b; exact)
+function _add_union_holes!(alg::FosterHormannClipping, polys, a_in_b, b_in_a, poly_a, poly_b; exact, crs)
     if a_in_b
-        _add_union_holes_contained_polys!(alg, polys, poly_a, poly_b; exact)
+        _add_union_holes_contained_polys!(alg, polys, poly_a, poly_b; exact, crs)
     elseif b_in_a
-        _add_union_holes_contained_polys!(alg, polys, poly_b, poly_a; exact)
+        _add_union_holes_contained_polys!(alg, polys, poly_b, poly_a; exact, crs)
     else  # Polygons intersect, but neither is contained in the other
         n_a_holes = GI.nhole(poly_a)
         ext_poly_a = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_a)))
@@ -166,7 +168,7 @@ end
 #= Add holes holes to the union of two polygons where one of the original polygons was
 inside of the other. If adding the the holes reveal that the polygons aren't actually
 intersecting, return the original polygons.=#
-function _add_union_holes_contained_polys!(alg::FosterHormannClipping, polys, interior_poly, exterior_poly; exact)
+function _add_union_holes_contained_polys!(alg::FosterHormannClipping, polys, interior_poly, exterior_poly; exact, crs)
     union_poly = polys[1]
     ext_int_ring = GI.getexterior(interior_poly)
     for (i, ih) in enumerate(GI.gethole(exterior_poly))
@@ -176,8 +178,8 @@ function _add_union_holes_contained_polys!(alg::FosterHormannClipping, polys, in
             if !on_ih && !out_ih
                 #= interior polygon is completely within the ith hole - polygons aren't
                 touching and do not actually form a union =#
-                polys[1] = tuples(interior_poly)
-                push!(polys, tuples(exterior_poly))
+                polys[1] = tuples(interior_poly; crs)
+                push!(polys, tuples(exterior_poly; crs))
                 return polys
             else
                 #= interior polygon is partially within the ith hole - area of interior
@@ -231,19 +233,20 @@ function _union(
     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)]
+    crs = mutual_crs(poly_a, multipoly_b)
+    polys = [tuples(poly_a, T; crs)]
     for poly_b in GI.getpolygon(multipoly_b)
         if intersects(#=TODO: alg.manifold, =#polys[1], poly_b)
             # If polygons intersect and form a new polygon, swap out polygon
             new_polys = union(alg, 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
+                push!(polys, tuples(poly_b, T; crs))  # 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))
+            push!(polys, tuples(poly_b, T; crs))
         end
     end
     return polys