Intersection points (#164)

* Start updating intersection point function * Improve intersection point calculation and add tests * Add test and comments * Add check for missing data in primitive data tests * Bump patch version * Remove Proj from Project.toml Co-authored-by: Anshul Singhvi <anshulsinghvi@gmail.com> --------- Co-authored-by: Anshul Singhvi <anshulsinghvi@gmail.com>
JuliaGeo · Jun 24, 2024 · 1e46749 · 1e46749 · skygering · Jun 24, 2024
1 parent 696326c
commit 1e46749
Show file tree

Hide file tree

Showing 5 changed files with 78 additions and 46 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "GeometryOps"
 uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 authors = ["Anshul Singhvi <anshulsinghvi@gmail.com> and contributors"]
-version = "0.1.9"
+version = "0.1.10"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
@@ -15,13 +15,13 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
-LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
+LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 
 [extensions]
-GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsFlexiJoinsExt = "FlexiJoins"
+GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsProjExt = "Proj"
 
 [compat]
@@ -50,11 +50,11 @@ GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
-Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
@@ -175,57 +175,53 @@ function _intersection(
 end
 
 """
-    intersection_points(
-        geom_a,
-        geom_b,
-    )::Union{
-        ::Vector{::Tuple{::Real, ::Real}},
-        ::Nothing,
-    }
-
-Return a list of intersection points between two geometries of type GI.Point.
-If no intersection point was possible given geometry extents, returns an empty
-list.
+    intersection_points(geom_a, geom_b, [T::Type])
+
+Return a list of intersection tuple points between two geometries. If no intersection points
+exist, returns an empty list.
+
+## Example
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
+line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
+inter_points = GO.intersection_points(line1, line2)
+
+# output
+1-element Vector{Tuple{Float64, Float64}}:
+ (125.58375366067548, -14.83572303404496)
 """
 intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat =
     _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
 
 
 #= Calculates the list of intersection points between two geometries, inlcuding line
-segments, line strings, linear rings, polygons, and multipolygons. If no intersection points
-were possible given geometry extents or if none are found, return an empty list of
-GI.Points. =#
-function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = _False()) where T
+segments, line strings, linear rings, polygons, and multipolygons. =#
+function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = _True()) where T
     # Initialize an empty list of points
-    result = GI.Point[]
+    result = Tuple{T, T}[]
     # Check if the geometries extents even overlap
     Extents.intersects(GI.extent(a), GI.extent(b)) || return result
     # Create a list of edges from the two input geometries
     edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
-    npoints_a, npoints_b  = length(edges_a), length(edges_b)
-    a_closed = npoints_a > 1 && edges_a[1][1] == edges_a[end][1]
-    b_closed = npoints_b > 1 && edges_b[1][1] == edges_b[end][1]
-    if npoints_a > 0 && npoints_b > 0
-        # Loop over pairs of edges and add any intersection points to results
-        for i in eachindex(edges_a), j in eachindex(edges_b)
-            line_orient, intr1, _ = _intersection_point(T, edges_a[i], edges_b[j]; exact)
-            # TODO: Add in degenerate intersection points when line_over
-            if line_orient == line_cross || line_orient == line_hinge
-                #=
-                Determine if point is on edge (all edge endpoints excluded
-                except for the last edge for an open geometry)
-                =#
-                point, (α, β) = intr1
-                on_a_edge = (!a_closed && i == npoints_a && 0 <= α <= 1) ||
-                    (0 <= α < 1)
-                on_b_edge = (!b_closed && j == npoints_b && 0 <= β <= 1) ||
-                    (0 <= β < 1)
-                if on_a_edge && on_b_edge
-                    push!(result, GI.Point(point))
-                end
-            end
+    # Loop over pairs of edges and add any unique intersection points to results
+    for a_edge in edges_a, b_edge in edges_b
+        line_orient, intr1, intr2 = _intersection_point(T, a_edge, b_edge; exact)
+        line_orient == line_out && continue  # no intersection points
+        pt1, _ = intr1
+        push!(result, pt1)  # if not line_out, there is at least one intersection point
+        if line_orient == line_over # if line_over, there are two intersection points
+            pt2, _ = intr2
+            push!(result, pt2)
         end
     end
+    #= TODO: We might be able to just add unique points with checks on the α and β values
+    returned from `_intersection_point`, but this would be different for curves vs polygons
+    vs multipolygons depending on if the shape is closed. This then wouldn't allow using the
+    `to_edges` functionality.  =# 
+    unique!(sort!(result))
     return result
 end
 

diff --git a/test/methods/clipping/intersection_points.jl b/test/methods/clipping/intersection_points.jl
@@ -0,0 +1,34 @@
+using Test
+import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+
+l1 = GI.LineString([(90000.0, 1000.0), (90000.0, 22500.0), (95000.0, 22500.0), (95000.0, 1000.0), (90000.0, 1000.0)])
+l2 = GI.LineString([(90000.0, 7500.0), (107500.0, 27500.0), (112500.0, 27500.0), (95000.0, 7500.0), (90000.0, 7500.0)])
+l3 = GI.LineString([(90000.0, 90000.0), (90000.0, 105000.0), (105000.0, 105000.0), (105000.0, 90000.0), (90000.0, 90000.0)])
+l4 = GI.LineString([(-98000.0, 90000.0), (-98000.0, 105000.0), (98000.0, 105000.0), (98000.0, 90000.0), (-98000.0, 90000.0)])
+l5 = GI.LineString([(19999.999, 25000.0), (19999.999, 29000.0), (39999.998999999996, 29000.0), (39999.998999999996, 25000.0), (19999.999, 25000.0)])
+l6 = GI.LineString([(0.0, 25000.0), (0.0, 29000.0), (20000.0, 29000.0), (20000.0, 25000.0), (0.0, 25000.0)])
+
+p1, p2 = GI.Polygon([l1]), GI.Polygon([l2])
+
+# Three intersection points
+GO_l1_l2_mp = GI.MultiPoint(GO.intersection_points(l1, l2))
+LG_l1_l2_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l1, l2))))
+@test GO.equals(GO_l1_l2_mp, LG_l1_l2_mp)
+
+# Four intersection points with large intersection
+GO_l3_l4_mp = GI.MultiPoint(GO.intersection_points(l3, l4))
+LG_l3_l4_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l3, l4))))
+@test GO.equals(GO_l3_l4_mp, LG_l3_l4_mp)
+
+# Four intersection points with very small intersection
+GO_l5_l6_mp = GI.MultiPoint(GO.intersection_points(l5, l6))
+LG_l5_l6_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l5, l6))))
+@test GO.equals(GO_l5_l6_mp, LG_l5_l6_mp)
+
+# Test that intersection points between lines and polygons is equivalent
+GO_p1_p2_mp = GI.MultiPoint(GO.intersection_points(p1, p2))
+@test GO.equals(GO_p1_p2_mp, GO_l1_l2_mp)
+
+# No intersection points between polygon and line
+GO_p1_l6 = GO.intersection_points(p1, l6)
+@test isempty(GO_p1_l6)
diff --git a/test/primitives.jl b/test/primitives.jl
@@ -21,6 +21,9 @@ poly = GI.Polygon([lr1, lr2])
                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
 
     @testset "Tables.jl support" begin
+        # check to account for missing data
+        missing_or_equal(x, y) = (ismissing(x) && ismissing(y)) || (x == y)
+        # file setup
         mktempdir() do dir
         cd(dir) do
 
@@ -37,7 +40,7 @@ poly = GI.Polygon([lr1, lr2])
                 @test all(centroid_geometry .== GO.centroid.(countries_table.geometry))
                 @testset "Columns are preserved" begin  
                     for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_table))
-                        @test all(GO.Tables.getcolumn(centroid_table, column) .== GO.Tables.getcolumn(countries_table, column))
+                        @test all(missing_or_equal.(GO.Tables.getcolumn(centroid_table, column), GO.Tables.getcolumn(countries_table, column)))
                     end
                 end
             end
@@ -51,7 +54,7 @@ poly = GI.Polygon([lr1, lr2])
                 @test all(centroid_geometry .== GO.centroid.(countries_df.geometry))
                 @testset "Columns are preserved" begin  
                     for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_df))
-                        @test all(centroid_df[!, column] .== countries_df[!, column])
+                        @test all(missing_or_equal.(centroid_df[!, column], countries_df[!, column]))
                     end
                 end
             end
@@ -60,8 +63,6 @@ poly = GI.Polygon([lr1, lr2])
     end
 end
 
-
-
 @testset "unwrap" begin
     flipped_vectors = GO.unwrap(GI.PointTrait, poly) do p
         (GI.y(p), GI.x(p))

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -28,6 +28,7 @@ const GO = GeometryOps
     # Clipping
     @testset "Coverage" begin include("methods/clipping/coverage.jl") end
     @testset "Cut" begin include("methods/clipping/cut.jl") end
+    @testset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end
     @testset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end
     # Transformations
     @testset "Embed Extent" begin include("transformations/extent.jl") end