fix crs on polygon clipping

GeometryOpsCore/Project.toml

@@ -1,7 +1,7 @@
 name = "GeometryOpsCore"
 uuid = "05efe853-fabf-41c8-927e-7063c8b9f013"
 authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"

GeometryOpsCore/src/geometry_utils.jl

@@ -1,3 +1,14 @@
 
+# Trait-based _linearring to handle any GeoInterface-compatible geometry
+_linearring(geom) = _linearring(GI.trait(geom), geom)
+
+# If it's already a LinearRing, return as-is
+_linearring(::GI.LinearRingTrait, geom) = geom
+
+# If it's a LineString (e.g., from ArchGDAL), convert to LinearRing preserving CRS
+_linearring(::GI.LineStringTrait, geom) =
+    GI.LinearRing(GI.getpoint(geom); crs=GI.crs(geom))
+
+# Concrete type specializations for GI wrappers (optimization)
 _linearring(geom::GI.LineString) = GI.LinearRing(parent(geom); extent=geom.extent, crs=geom.crs)
 _linearring(geom::GI.LinearRing) = geom

Project.toml

@@ -1,7 +1,7 @@
 name = "GeometryOps"
 uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
-version = "0.1.31"
+version = "0.1.32"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -21,6 +21,9 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
+[sources]
+GeometryOpsCore = {path = "GeometryOpsCore"}
+
 [weakdeps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
@@ -48,7 +51,7 @@ FlexiJoins = "0.1.30"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7, 0.5"
-GeometryOpsCore = "=0.1.7"
+GeometryOpsCore = "=0.1.8"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"

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/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/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/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)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b; crs)
     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([_linearring(tuples(ext_a; crs))]; crs))
         elseif b_in_a
-            push!(polys, GI.Polygon([tuples(ext_b)]))
+            push!(polys, GI.Polygon([_linearring(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/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

src/utils/utils.jl

@@ -296,3 +296,30 @@ 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 isnothing(GI.crs(a)) && !isnothing(GI.crs(b))
+        return GI.crs(b)
+    elseif !isnothing(GI.crs(a)) && isnothing(GI.crs(b))
+        return GI.crs(a)
+    else
+        error("""
+        Cannot perform clipping operation on two geometries with different CRS! 
+
+        Either pass an explicit `crs` keyword to the toplevel function you called (`intersection, difference, union`) or make sure the geometries you're clipping have the same crs.
+
+        Got the CRSes:
+        $(GI.crs(a))
+        and
+        $(GI.crs(b))
+        """)
+    end
+end