diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index 94bf4e859..1b376c6ff 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -630,11 +630,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 @@ -689,15 +690,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) diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl index 47d82077c..78dfb5a94 100644 --- a/src/methods/clipping/intersection.jl +++ b/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 diff --git a/src/utils/utils.jl b/src/utils/utils.jl index db27124c0..dc04fd100 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -289,5 +289,21 @@ end +""" + 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