fix crs on intersection

Author: asinghvi17

src/methods/clipping/clipping_processor.jl

@@ -937,11 +937,12 @@ A list of GeoInterface polygons is returned from this function.
 Note: `poly_a` and `poly_b` are temporary inputs used for debugging and can be removed
 eventually.
 =#
-function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b) where {T, M, A}
+function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b; crs = mutual_crs(poly_a, poly_b)) where {T, M, A}
     n_a_pts, n_b_pts = length(a_list), length(b_list)
     total_pts = n_a_pts + n_b_pts
     n_cross_pts = length(a_idx_list)
-    return_polys = Vector{_get_poly_type(T)}(undef, 0)
+
+    return_polys = Vector{_get_poly_type(T, crs)}(undef, 0)
     # Keep track of number of processed intersection points
     visited_pts = 0
     processed_pts = 0
@@ -996,15 +997,16 @@ function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b
             idx = curr.neighbor
             curr = curr_list[idx]
         end
-        push!(return_polys, GI.Polygon([pt_list]))
+        push!(return_polys, GI.Polygon([GI.LinearRing(pt_list; crs)]; crs))
     end
     return return_polys
 end
 
 # Get type of polygons that will be made
 # TODO: Increase type options
-_get_poly_type(::Type{T}) where T =
-    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing}
+function _get_poly_type(::Type{T}, crs::CRST = nothing) where {T, CRST}
+    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, CRST}}, Nothing, CRST}
+end
 
 #=
     _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)

src/methods/clipping/intersection.jl

@@ -80,15 +80,17 @@ function _intersection(
     # First we get the exteriors of 'poly_a' and 'poly_b'
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
+
+    crs = mutual_crs(poly_a, poly_b)
     # Then we find the intersection of the exteriors
     a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact)
     polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b)
     if isempty(polys) # 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([tuples(ext_a)]))
+            push!(polys, GI.Polygon([tuples(ext_a; crs)]; crs))
         elseif b_in_a
-            push!(polys, GI.Polygon([tuples(ext_b)]))
+            push!(polys, GI.Polygon([tuples(ext_b; crs)]; crs))
         end
     end
     remove_idx = falses(length(polys))
@@ -130,7 +132,8 @@ function _intersection(
     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)}()
+    crs = mutual_crs(poly_a, multipoly_b)
+    polys = Vector{_get_poly_type(T, crs)}()
     for poly_b in GI.getpolygon(multipoly_b)
         append!(polys, intersection(alg, poly_a, poly_b; target))
     end
@@ -163,7 +166,8 @@ function _intersection(
         multipoly_b = fix_multipoly(multipoly_b)
         fix_multipoly = nothing
     end
-    polys = Vector{_get_poly_type(T)}()
+    crs = mutual_crs(multipoly_a, multipoly_b)
+    polys = Vector{_get_poly_type(T, crs)}()
     for poly_a in GI.getpolygon(multipoly_a)
         append!(polys, intersection(alg, poly_a, multipoly_b; target, fix_multipoly))
     end

src/utils/utils.jl

@@ -296,3 +296,21 @@ end
 _geometry_or_error(g::Extents.Extent; kw...) = g
 
 
+"""
+    mutual_crs(a, b)
+
+Return the CRS of `a` if it is the same as `b`, or `nothing` if they are not the same.
+
+This is tolerant of `nothing` values in one of a or b, it assumes that they are the same crs.
+"""
+function mutual_crs(a, b)
+    if GI.crs(a) == GI.crs(b)
+        return GI.crs(a)
+    elseif GI.crs(a) == nothing && GI.crs(b) != nothing
+        return GI.crs(b)
+    elseif GI.crs(a) != nothing && GI.crs(b) == nothing
+        return GI.crs(a)
+    else
+        return nothing
+    end
+end