From aba090f5cd87b329629c7392b65af9b28d90335a Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Mon, 2 Oct 2023 16:55:57 -0700 Subject: [PATCH 01/16] Fix up intersection point base calculation --- src/methods/intersects.jl | 162 +++++++++++++++++++++++++++++--------- 1 file changed, 126 insertions(+), 36 deletions(-) diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl index e0476bd81..ca20ed321 100644 --- a/src/methods/intersects.jl +++ b/src/methods/intersects.jl @@ -2,14 +2,66 @@ export intersects, intersection -# This code checks whether geometries intersect with each other. +#= +## What is `intersects` vs `intersection`? + +The `intersects` methods check whether two geometries intersect with each other. +The `intersection` methods return the intersection between the two geometries. + +The `intersects` methods will always return a Boolean. However, note that the +`intersection` methods will not all return the same type. For example, the +intersection of two lines will be a point in most cases, unless the lines are +parallel. On the other hand, the intersection of two polygons will be another +polygon in most cases. + +To provide an example, consider this # TODO update this example: +```@example cshape +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +cshape = Polygon([ + Point(0,0), Point(0,3), Point(3,3), Point(3,2), Point(1,2), + Point(1,1), Point(3,1), Point(3,0), Point(0,0), +]) +f, a, p = poly(cshape; axis = (; aspect = DataAspect())) +``` +Let's see what the centroid looks like (plotted in red): +```@example cshape +cent = centroid(cshape) +scatter!(a, GI.x(cent), GI.y(cent), color = :red) +f +``` + +## Implementation -# !!! note -# This does not compute intersections, only checks if they exist. +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. This is also used in the +implementation, since it's a lot less work! + +# TODO fill this in! +=# const MEETS_OPEN = 1 const MEETS_CLOSED = 0 +intersects(geom1, geom2) = GO.intersects( + GI.trait(geom1), + geom1, + GI.trait(geom2), + geom2, +) + +GO.intersects( + trait1::Union{GI.LineStringTrait, GI.LinearRingTrait}, + geom1, + trait2::Union{GI.LineStringTrait, GI.LinearRingTrait}, + geom2, +) = line_intersects(trait1, geom1, trait2, geom2) + """ line_intersects(line_a, line_b) @@ -73,53 +125,91 @@ GO.line_intersection(line1, line2) (125.58375366067547, -14.83572303404496) ``` """ -line_intersection(line_a, line_b) = line_intersection(trait(line_a), line_a, trait(line_b), line_b) -function line_intersection(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) +line_intersection(line_a, line_b) = intersection_points(trait(line_a), line_a, trait(line_b), line_b) + +""" + intersection_points( + ::GI.AbstractTrait, geom_a, + ::GI.AbstractTrait, geom_b, + )::Vector{::Tuple{::Real, ::Real}} + +Calculates the list of intersection points between two geometries. +""" +function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing result = Tuple{Float64,Float64}[] edges_a, edges_b = map(sort! ∘ to_edges, (a, b)) for edge_a in edges_a for edge_b in edges_b - x = _line_intersection(edge_a, edge_b) + x = _intersection_point(edge_a, edge_b) isnothing(x) || push!(result, x) end end return result end -function line_intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) + +""" + intersection_point( + ::GI.LineTrait, line_a, + ::GI.LineTrait, line_b, + )::Union{ + ::Tuple{::Real, ::Real}, + ::Nothing + } + +Calculates the intersection point between two lines if it exists and return +`nothing` if it doesn't exist. +""" +function intersection_point(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) + # Get start and end points for both lines a1 = GI.getpoint(line_a, 1) - b1 = GI.getpoint(line_b, 1) a2 = GI.getpoint(line_a, 2) + b1 = GI.getpoint(line_b, 1) b2 = GI.getpoint(line_b, 2) - - return _line_intersection((a1, a2), (b1, b2)) + # Determine the intersection point + point, _ = _intersection_point((a1, a2), (b1, b2)) + return point end -function _line_intersection((p11, p12)::Tuple, (p21, p22)::Tuple) - # Get points from lines - x1, y1 = GI.x(p11), GI.y(p11) - x2, y2 = GI.x(p12), GI.y(p12) - x3, y3 = GI.x(p21), GI.y(p21) - x4, y4 = GI.x(p22), GI.y(p22) - - d = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1)) - a = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3)) - b = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3)) - - if d == 0 - if a == 0 && b == 0 - return nothing - end - return nothing - end - - ã = a / d - b̃ = b / d - if ã >= 0 && ã <= 1 && b̃ >= 0 && b̃ <= 1 - x = x1 + (ã * (x2 - x1)) - y = y1 + (ã * (y2 - y1)) - return (x, y) +""" + _intersection_point( + (p11, p12)::Tuple, + (p21, p22)::Tuple, + ) + +Calculates the intersection point between two lines if it exists, and the +fractional component of each line from the initial end point to the +intersection point. +Inputs: + (p11, p12)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line + (p21, p22)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line +Outputs: + (x, y)::Tuple{::Real, ::Real} intersection point + (t, u)::Tuple{::Real, ::Real} fractional length of lines to intersection + Both are ::Nothing if point doesn't exist! + +Calculation derivation can be found here: + https://stackoverflow.com/questions/563198/ +""" +function _intersection_point((p11, p12)::Tuple, (p21, p22)::Tuple) + # First line runs from p to p + r + px, py = GI.x(p11), GI.y(p11) + rx, ry = GI.x(p12) - px, GI.y(p12) - py + # Second line runs from q to q + s + qx, qy = GI.x(p21), GI.y(p21) + sx, sy = GI.x(p22) - qx, GI.y(p22) - qy + # Intersection will be where p + tr = q + us where 0 < t, u < 1 and + r_cross_s = rx * sy - ry * sx + if r_cross_s != 0 + Δpq_x = px - qx + Δpq_y = py - qy + t = (Δpq_x * sy - Δpq_y * sx) / r_cross_s + u = (Δpq_x * ry - Δpq_y * rx) / r_cross_s + if 0 <= t <= 1 && 0 <= u <= 1 + x = px + t * rx + y = py + t * ry + return (x, y), (t, u) + end end - - return nothing + return nothing, nothing end From ab167851686704fcde6d81923cad1ba57fee2bee Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Tue, 3 Oct 2023 18:18:36 -0700 Subject: [PATCH 02/16] Update intersects and add line tests --- Project.toml | 6 +- src/methods/intersects.jl | 337 +++++++++++++++++++++++----------- src/transformations/extent.jl | 2 +- src/utils.jl | 2 +- test/methods/intersects.jl | 108 +++++++++++ test/runtests.jl | 1 + 6 files changed, 346 insertions(+), 110 deletions(-) create mode 100644 test/methods/intersects.jl diff --git a/Project.toml b/Project.toml index 318c474d8..0fe3b53d8 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi and contributors"] version = "0.0.1-DEV" [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" @@ -26,7 +27,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = [ - "ArchGDAL", "Distributions", "GeoFormatTypes", "GeoJSON", "LibGEOS", - "Random", "Test", -] +test = ["ArchGDAL", "Distributions", "GeoFormatTypes", "GeoJSON", "LibGEOS", "Random", "Test"] diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl index ca20ed321..6f3ca4abc 100644 --- a/src/methods/intersects.jl +++ b/src/methods/intersects.jl @@ -1,36 +1,42 @@ # # Intersection checks -export intersects, intersection +export intersects, intersection, intersection_points #= -## What is `intersects` vs `intersection`? +## What is `intersects` vs `intersection` vs `intersection_points`? The `intersects` methods check whether two geometries intersect with each other. -The `intersection` methods return the intersection between the two geometries. +The `intersection` methods return the geometry intersection between the two +input geometries. The `intersection_points` method returns a list of +intersection points between two geometries. The `intersects` methods will always return a Boolean. However, note that the `intersection` methods will not all return the same type. For example, the intersection of two lines will be a point in most cases, unless the lines are parallel. On the other hand, the intersection of two polygons will be another -polygon in most cases. +polygon in most cases. Finally, the `intersection_points` method returns a list +of tuple points. -To provide an example, consider this # TODO update this example: -```@example cshape +To provide an example, consider these two lines: +```@example intersects_intersection using GeometryOps using GeometryOps.GeometryBasics using Makie using CairoMakie - -cshape = Polygon([ - Point(0,0), Point(0,3), Point(3,3), Point(3,2), Point(1,2), - Point(1,1), Point(3,1), Point(3,0), Point(0,0), -]) -f, a, p = poly(cshape; axis = (; aspect = DataAspect())) +point1, point2 = Point(124.584961,-12.768946), Point(126.738281,-17.224758) +point3, point4 = Point(123.354492,-15.961329), Point(127.22168,-14.008696) +line1 = Line(point1, point2) +line2 = Line(point3, point4) +f, a, p = lines([point1, point2]) +lines!([point3, point4]) ``` -Let's see what the centroid looks like (plotted in red): -```@example cshape -cent = centroid(cshape) -scatter!(a, GI.x(cent), GI.y(cent), color = :red) +We can see that they intersect, so we expect intersects to return true, and we +can visualize the intersection point in red. +```@example intersects_intersection +int_bool = GO.intersects(line1, line2) +println(int_bool) +int_point = GO.intersection(line1, line2) +scatter!(int_point, color = :red) f ``` @@ -38,36 +44,24 @@ f This is the GeoInterface-compatible implementation. -First, we implement a wrapper method that dispatches to the correct -implementation based on the geometry trait. This is also used in the -implementation, since it's a lot less work! - -# TODO fill this in! +First, we implement a wrapper method for intersects, intersection, and +intersection_points that dispatches to the correct implementation based on the +geometry trait. The two underlying helper functions that are widely used in all +geometry dispatches are _line_intersects, which determines if two line segments +intersect and _intersection_point which determines the intersection point +between two line segments. =# -const MEETS_OPEN = 1 const MEETS_CLOSED = 0 - -intersects(geom1, geom2) = GO.intersects( - GI.trait(geom1), - geom1, - GI.trait(geom2), - geom2, -) - -GO.intersects( - trait1::Union{GI.LineStringTrait, GI.LinearRingTrait}, - geom1, - trait2::Union{GI.LineStringTrait, GI.LinearRingTrait}, - geom2, -) = line_intersects(trait1, geom1, trait2, geom2) +const MEETS_OPEN = 1 """ - line_intersects(line_a, line_b) - -Check if `line_a` intersects with `line_b`. + intersects(geom1, geom2; kw...)::Bool -These can be `LineTrait`, `LineStringTrait` or `LinearRingTrait` +Check if two geometries intersect, returning true if so and false otherwise. +Takes in a Int keyword meets, which can either be MEETS_OPEN (1), meaning that +only intersections through open edges where edge endpoints are not included are +recorded, versus MEETS_CLOSED (0) where edge endpoints are included. ## Example @@ -76,41 +70,80 @@ 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)]) -GO.line_intersects(line1, line2) +GO.intersects(line1, line2) # output true ``` """ -line_intersects(a, b; kw...) = line_intersects(trait(a), a, trait(b), b; kw...) -# Skip to_edges for LineTrait -function line_intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets=MEETS_OPEN) +intersects(geom1, geom2; kw...) = intersects( + GI.trait(geom1), + geom1, + GI.trait(geom2), + geom2; + kw... +) + +""" + intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets = MEETS_OPEN)::Bool + +Returns true if two line segments intersect and false otherwise. Line segment +endpoints are excluded in check if `meets = MEETS_OPEN` (1) and included if +`meets = MEETS_CLOSED` (0). +""" +function intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets = MEETS_OPEN) a1 = _tuple_point(GI.getpoint(a, 1)) - b1 = _tuple_point(GI.getpoint(b, 1)) a2 = _tuple_point(GI.getpoint(a, 2)) + b1 = _tuple_point(GI.getpoint(b, 1)) b2 = _tuple_point(GI.getpoint(b, 2)) - return ExactPredicates.meet(a1, a2, b1, b2) == meets + meet_type = ExactPredicates.meet(a1, a2, b1, b2) + return meet_type == MEETS_OPEN || meet_type == meets end -function line_intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...) + +""" + intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...)::Bool + +Returns true if two geometries intersect with one another and false +otherwise. For all geometries but lines, conver the geometry to a list of edges +and cross compare the edges for intersections. +""" +function intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...) edges_a, edges_b = map(sort! ∘ to_edges, (a, b)) - return line_intersects(edges_a, edges_b; kw...) + return _line_intersects(edges_a, edges_b; kw...) end -function line_intersects(edges_a::Vector{Edge}, edges_b::Vector{Edge}; meets=MEETS_OPEN) + +""" + _line_intersects( + edges_a::Vector{Edge}, + edges_b::Vector{Edge}; + meets = MEETS_OPEN, + )::Bool + +Returns true if there is at least one intersection between edges within the +two lists. Line segment endpoints are excluded in check if `meets = MEETS_OPEN` +(1) and included if `meets = MEETS_CLOSED` (0). +""" +function _line_intersects( + edges_a::Vector{Edge}, + edges_b::Vector{Edge}; + meets = MEETS_OPEN, +) # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false for edge_a in edges_a for edge_b in edges_b - ExactPredicates.meet(edge_a..., edge_b...) == meets && return true + meet_type = ExactPredicates.meet(edge_a..., edge_b...) + (meet_type == MEETS_OPEN || meet_type == meets) && return true end end return false end """ - line_intersection(line_a, line_b) - -Find a point that intersects LineStrings with two coordinates each. + intersection(geom_a, geom_b)::Union{Tuple{::Real, ::Real}, ::Nothing} -Returns `nothing` if no point is found. +Return an intersection point between two geometries. Return nothing if none are +found. Else, the return type depends on the input. It will be a union between: +a point, a line, a linear ring, a polygon, or a multipolygon ## Example @@ -119,37 +152,17 @@ 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)]) -GO.line_intersection(line1, line2) +GO.intersection(line1, line2) # output (125.58375366067547, -14.83572303404496) ``` """ -line_intersection(line_a, line_b) = intersection_points(trait(line_a), line_a, trait(line_b), line_b) - -""" - intersection_points( - ::GI.AbstractTrait, geom_a, - ::GI.AbstractTrait, geom_b, - )::Vector{::Tuple{::Real, ::Real}} - -Calculates the list of intersection points between two geometries. -""" -function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) - Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing - result = Tuple{Float64,Float64}[] - edges_a, edges_b = map(sort! ∘ to_edges, (a, b)) - for edge_a in edges_a - for edge_b in edges_b - x = _intersection_point(edge_a, edge_b) - isnothing(x) || push!(result, x) - end - end - return result -end +intersection(geom_a, geom_b) = + intersection(GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) """ - intersection_point( + intersection( ::GI.LineTrait, line_a, ::GI.LineTrait, line_b, )::Union{ @@ -157,32 +170,150 @@ end ::Nothing } -Calculates the intersection point between two lines if it exists and return -`nothing` if it doesn't exist. +Calculates the intersection between two line segments. Return nothing if +there isn't one. """ -function intersection_point(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) +function intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) # Get start and end points for both lines a1 = GI.getpoint(line_a, 1) a2 = GI.getpoint(line_a, 2) b1 = GI.getpoint(line_b, 1) b2 = GI.getpoint(line_b, 2) # Determine the intersection point - point, _ = _intersection_point((a1, a2), (b1, b2)) - return point + point, fracs = _intersection_point((a1, a2), (b1, b2)) + # Determine if intersection point is on line segments + if !isnothing(point) && 0 <= fracs[1] <= 1 && 0 <= fracs[2] <= 1 + return point + end + return nothing +end + +intersection( + trait_a::Union{GI.LineStringTrait, GI.LinearRingTrait}, + geom_a, + trait_b::Union{GI.LineStringTrait, GI.LinearRingTrait}, + geom_b, +) = intersection_points(trait_a, geom_a, trait_b, geom_b) + +""" + intersection( + ::GI.PolygonTrait, poly_a, + ::GI.PolygonTrait, poly_b, + )::Union{ + ::Vector{Vector{Tuple{::Real, ::Real}}}, # is this a good return type? + ::Nothing + } + +Calculates the intersection between two line segments. Return nothing if +there isn't one. +""" +function intersection(::GI.PolygonTrait, poly_a, ::GI.PolygonTrait, poly_b) + @assert false "Polygon intersection isn't implemented yet." + return nothing +end + +""" + intersection( + ::GI.AbstractTrait, geom_a, + ::GI.AbstractTrait, geom_b, + )::Union{ + ::Vector{Vector{Tuple{::Real, ::Real}}}, # is this a good return type? + ::Nothing + } + +Calculates the intersection between two line segments. Return nothing if +there isn't one. +""" +function intersection( + trait_a::GI.AbstractTrait, geom_a, + trait_b::GI.AbstractTrait, geom_b, +) + @assert( + false, + "Intersection between $trait_a and $trait_b isn't implemented yet.", + ) + return nothing +end + +""" + intersection_points( + geom_a, + geom_b, + )::Union{ + ::Vector{::Tuple{::Real, ::Real}}, + ::Nothing, + } + +Return a list of intersection points between two geometries. If no intersection +point was possible given geometry extents, return nothing. If none are found, +return an empty list. +""" +intersection_points(geom_a, geom_b) = + intersection_points(GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) + +""" + intersection_points( + ::GI.AbstractTrait, geom_a, + ::GI.AbstractTrait, geom_b, + )::Union{ + ::Vector{::Tuple{::Real, ::Real}}, + ::Nothing, + } + +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, return nothing. If +none are found, return an empty list. +""" +function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) + # Check if the geometries extents even overlap + Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing + # 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 = edges_a[1][1] == edges_a[end][1] + b_closed = edges_b[1][1] == edges_b[end][1] + if npoints_a > 0 && npoints_b > 0 + # Initialize an empty list of points + T = typeof(edges_a[1][1][1]) # x-coordinate of first point in first edge + result = Tuple{T,T}[] + # Loop over pairs of edges and add any intersection points to results + for i in eachindex(edges_a) + for j in eachindex(edges_b) + point, fracs = _intersection_point(edges_a[i], edges_b[j]) + if !isnothing(point) + #= + Determine if point is on edge (all edge endpoints excluded + except for the last edge for an open geometry) + =# + α, β = fracs + 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, point) + end + end + end + end + return result + end + return nothing end """ _intersection_point( - (p11, p12)::Tuple, - (p21, p22)::Tuple, + (a1, a2)::Tuple, + (b1, b2)::Tuple, ) -Calculates the intersection point between two lines if it exists, and the -fractional component of each line from the initial end point to the -intersection point. +Calculates the intersection point between two lines if it exists, and as if the +line extended to infinity, and the fractional component of each line from the +initial end point to the intersection point. Inputs: - (p11, p12)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line - (p21, p22)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line + (a1, a2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line + (b1, b2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line Outputs: (x, y)::Tuple{::Real, ::Real} intersection point (t, u)::Tuple{::Real, ::Real} fractional length of lines to intersection @@ -191,25 +322,23 @@ Outputs: Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ """ -function _intersection_point((p11, p12)::Tuple, (p21, p22)::Tuple) +function _intersection_point((a1, a2)::Tuple, (b1, b2)::Tuple) # First line runs from p to p + r - px, py = GI.x(p11), GI.y(p11) - rx, ry = GI.x(p12) - px, GI.y(p12) - py + px, py = GI.x(a1), GI.y(a1) + rx, ry = GI.x(a2) - px, GI.y(a2) - py # Second line runs from q to q + s - qx, qy = GI.x(p21), GI.y(p21) - sx, sy = GI.x(p22) - qx, GI.y(p22) - qy + qx, qy = GI.x(b1), GI.y(b1) + sx, sy = GI.x(b2) - qx, GI.y(b2) - qy # Intersection will be where p + tr = q + us where 0 < t, u < 1 and r_cross_s = rx * sy - ry * sx if r_cross_s != 0 - Δpq_x = px - qx - Δpq_y = py - qy - t = (Δpq_x * sy - Δpq_y * sx) / r_cross_s - u = (Δpq_x * ry - Δpq_y * rx) / r_cross_s - if 0 <= t <= 1 && 0 <= u <= 1 - x = px + t * rx - y = py + t * ry - return (x, y), (t, u) - end + Δqp_x = qx - px + Δqp_y = qy - py + t = (Δqp_x * sy - Δqp_y * sx) / r_cross_s + u = (Δqp_x * ry - Δqp_y * rx) / r_cross_s + x = px + t * rx + y = py + t * ry + return (x, y), (t, u) end return nothing, nothing end diff --git a/src/transformations/extent.jl b/src/transformations/extent.jl index 2e230672d..a5e180d76 100644 --- a/src/transformations/extent.jl +++ b/src/transformations/extent.jl @@ -10,7 +10,7 @@ embed_extent(x) = apply(extent_applicator, AbstractTrait, x) extent_applicator(x) = extent_applicator(trait(x), x) extent_applicator(::Nothing, xs::AbstractArray) = embed_extent.(xs) -function extent_applicator(::Union{AbstractCurveTrait,MultiPointTrait}, point) = point +extent_applicator(::Union{AbstractCurveTrait,MultiPointTrait}, point) = point function extent_applicator(trait::AbstractGeometryTrait, geom) children_with_extents = map(GI.getgeom(geom)) do g diff --git a/src/utils.jl b/src/utils.jl index dc6c078eb..b95b78172 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -54,7 +54,7 @@ end to_edges() Convert any geometry or collection of geometries into a flat -vector of `Tuple{Tuple{Float64,Float64},{Float64,Float64}}` edges. +vector of `Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}` edges. """ function to_edges(x) edges = Vector{Edge}(undef, _nedge(x)) diff --git a/test/methods/intersects.jl b/test/methods/intersects.jl new file mode 100644 index 000000000..33ece90cf --- /dev/null +++ b/test/methods/intersects.jl @@ -0,0 +1,108 @@ +@testset "Lines/Rings" begin + # Line test intersects ----------------------------------------------------- + + # Test for parallel lines + l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) + l2 = GI.Line([(0.0, 1.0), (2.5, 1.0)]) + @test !GO.intersects(l1, l2; meets = 0) + @test !GO.intersects(l1, l2; meets = 1) + @test isnothing(GO.intersection(l1, l2)) + + # Test for non-parallel lines that don't intersect + l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) + l2 = GI.Line([(2.0, -3.0), (3.0, 0.0)]) + @test !GO.intersects(l1, l2; meets = 0) + @test !GO.intersects(l1, l2; meets = 1) + @test isnothing(GO.intersection(l1, l2)) + + # Test for lines only touching at endpoint + l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) + l2 = GI.Line([(2.0, -3.0), (2.5, 0.0)]) + @test GO.intersects(l1, l2; meets = 0) + @test !GO.intersects(l1, l2; meets = 1) + @test all(GO.intersection(l1, l2) .≈ (2.5, 0.0)) + + # Test for lines that intersect in the middle + l1 = GI.Line([(0.0, 0.0), (5.0, 5.0)]) + l2 = GI.Line([(0.0, 5.0), (5.0, 0.0)]) + @test GO.intersects(l1, l2; meets = 0) + @test GO.intersects(l1, l2; meets = 1) + @test all(GO.intersection(l1, l2) .≈ (2.5, 2.5)) + + # Line string test intersects ---------------------------------------------- + + # Single element line strings crossing over each other + l1 = LG.LineString([[5.5, 7.2], [11.2, 12.7]]) + l2 = LG.LineString([[4.3, 13.3], [9.6, 8.1]]) + @test GO.intersects(l1, l2; meets = 0) + @test GO.intersects(l1, l2; meets = 1) + go_inter = GO.intersection(l1, l2) + lg_inter = LG.intersection(l1, l2) + @test go_inter[1][1] .≈ GI.x(lg_inter) + @test go_inter[1][2] .≈ GI.y(lg_inter) + + # Multi-element line strings crossing over on vertex + l1 = LG.LineString([[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]]) + l2 = LG.LineString([[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]]) + @test GO.intersects(l1, l2; meets = 0) + # TODO: Do we want this to be false? It is vertex of segment, not of whole line string + @test !GO.intersects(l1, l2; meets = 1) + go_inter = GO.intersection(l1, l2) + @test length(go_inter) == 1 + lg_inter = LG.intersection(l1, l2) + @test go_inter[1][1] .≈ GI.x(lg_inter) + @test go_inter[1][2] .≈ GI.y(lg_inter) + + # Multi-element line strings crossing over with multiple intersections + l1 = LG.LineString([[0.0, -1.0], [1.0, 1.0], [2.0, -1.0], [3.0, 1.0]]) + l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) + @test GO.intersects(l1, l2; meets = 0) + @test GO.intersects(l1, l2; meets = 1) + go_inter = GO.intersection(l1, l2) + @test length(go_inter) == 3 + lg_inter = LG.intersection(l1, l2) + @test issetequal( + Set(go_inter), + Set(GO._tuple_point.(GI.getpoint(lg_inter))) + ) + + # Line strings far apart so extents don't overlap + + # Line strings close together that don't overlap + + # Line string with empty line string + + # Closed linear ring with open line string + + # Closed linear ring with closed linear ring + + # @test issetequal( + # Subzero.intersect_lines(l1, l2), + # Set([(0.5, -0.0), (1.5, 0), (2.5, -0.0)]), + # ) + # l2 = [[[10., 10]]] + # @test issetequal( + # Subzero.intersect_lines(l1, l2), + # Set{Tuple{Float64, Float64}}(), + # ) + + +end + +@testset "Polygons" begin + # Two polygons that intersect + + # Two polygons that don't intersect + + # Polygon that intersects with linestring + +end + +@testset "MultiPolygons" begin + # Multi-polygon and polygon that intersect + + # Multi-polygon and polygon that don't intersect + + # Multi-polygon that intersects with linestring + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index c4fc39f3e..7c96de785 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ const GO = GeometryOps @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end @testset "Bools" begin include("methods/bools.jl") end @testset "Centroid" begin include("methods/centroid.jl") end + @testset "Intersect" begin include("methods/intersects.jl") end @testset "Signed Area" begin include("methods/signed_area.jl") end # Transformations @testset "Reproject" begin include("transformations/reproject.jl") end From 7e75e86ff51cf8ded2ebeeae2974fb05ed5cadf7 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 3 Oct 2023 22:34:42 -0700 Subject: [PATCH 03/16] Add more tests and debug intersects --- src/methods/bools.jl | 26 ++++++------ src/methods/crosses.jl | 4 +- src/methods/disjoint.jl | 2 +- src/methods/intersects.jl | 13 ++++-- src/methods/overlaps.jl | 4 +- src/methods/within.jl | 14 ++++++- test/methods/bools.jl | 6 +-- test/methods/intersects.jl | 84 +++++++++++++++++++++++++++++++------- 8 files changed, 111 insertions(+), 42 deletions(-) diff --git a/src/methods/bools.jl b/src/methods/bools.jl index aac4f8075..ba5c4068d 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -365,19 +365,19 @@ function line_in_polygon( end function polygon_in_polygon(poly1, poly2) - # edges1, edges2 = to_edges(poly1), to_edges(poly2) - # extent1, extent2 = to_extent(edges1), to_extent(edges2) - # Check the extents intersect - Extents.intersects(GI.extent(poly1), GI.extent(poly2)) || return false - - # Check all points in poly1 are in poly2 - for point in GI.getpoint(poly1) - point_in_polygon(point, poly2) || return false - end + # edges1, edges2 = to_edges(poly1), to_edges(poly2) + # extent1, extent2 = to_extent(edges1), to_extent(edges2) + # Check the extents intersect + Extents.intersects(GI.extent(poly1), GI.extent(poly2)) || return false + + # Check all points in poly1 are in poly2 + for point in GI.getpoint(poly1) + point_in_polygon(point, poly2) || return false + end - # Check the line of poly1 does not intersect the line of poly2 - line_intersects(poly1, poly2) && return false + # Check the line of poly1 does not intersect the line of poly2 + #intersects(poly1, poly2) && return false - # poly1 must be in poly2 - return true + # poly1 must be in poly2 + return true end diff --git a/src/methods/crosses.jl b/src/methods/crosses.jl index 7c215c857..3aa62d62e 100644 --- a/src/methods/crosses.jl +++ b/src/methods/crosses.jl @@ -55,7 +55,7 @@ end function line_crosses_line(line1, line2) np2 = GI.npoint(line2) - if line_intersects(line1, line2; meets=MEETS_CLOSED) + if intersects(line1, line2; meets=MEETS_CLOSED) for i in 1:GI.npoint(line1) - 1 for j in 1:GI.npoint(line2) - 1 exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both @@ -71,7 +71,7 @@ end function line_crosses_poly(line, poly) for l in flatten(AbstractCurveTrait, poly) - line_intersects(line, l) && return true + intersects(line, l) && return true end return false end diff --git a/src/methods/disjoint.jl b/src/methods/disjoint.jl index b51e5ab66..02b7d4e46 100644 --- a/src/methods/disjoint.jl +++ b/src/methods/disjoint.jl @@ -38,5 +38,5 @@ function polygon_disjoint(poly1, poly2) for point in GI.getpoint(poly2) point_in_polygon(point, poly1) && return false end - return !line_intersects(poly1, poly2) + return !intersects(poly1, poly2) end diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl index 6f3ca4abc..78f5784f1 100644 --- a/src/methods/intersects.jl +++ b/src/methods/intersects.jl @@ -107,9 +107,14 @@ Returns true if two geometries intersect with one another and false otherwise. For all geometries but lines, conver the geometry to a list of edges and cross compare the edges for intersections. """ -function intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...) +function intersects( + trait_a::GI.AbstractTrait, a, + trait_b::GI.AbstractTrait, b; + kw..., +) edges_a, edges_b = map(sort! ∘ to_edges, (a, b)) - return _line_intersects(edges_a, edges_b; kw...) + return _line_intersects(edges_a, edges_b; kw...) || + within(trait_a, a, trait_b, b) || within(trait_b, b, trait_a, a) end """ @@ -271,8 +276,8 @@ function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) # 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 = edges_a[1][1] == edges_a[end][1] - b_closed = edges_b[1][1] == edges_b[end][1] + 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 # Initialize an empty list of points T = typeof(edges_a[1][1][1]) # x-coordinate of first point in first edge diff --git a/src/methods/overlaps.jl b/src/methods/overlaps.jl index b846e43de..6d84f393b 100644 --- a/src/methods/overlaps.jl +++ b/src/methods/overlaps.jl @@ -37,11 +37,11 @@ function overlaps(::MultiPointTrait, g1, ::MultiPointTrait, g2)::Bool end end function overlaps(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool - return line_intersects(g1, g2) + return intersects(g1, g2) end function overlaps(t1::MultiPolygonTrait, mp, t2::PolygonTrait, p1)::Bool for p2 in GI.getgeom(mp) - overlaps(p1, thp2) && return true + overlaps(p1, p2) && return true end end function overlaps(::MultiPolygonTrait, g1, ::MultiPolygonTrait, g2)::Bool diff --git a/src/methods/within.jl b/src/methods/within.jl index c930ce62f..16366f944 100644 --- a/src/methods/within.jl +++ b/src/methods/within.jl @@ -23,11 +23,21 @@ GO.within(point, line) true ``` """ +# Syntactic sugar within(g1, g2)::Bool = within(trait(g1), g1, trait(g2), g2)::Bool within(::GI.FeatureTrait, g1, ::Any, g2)::Bool = within(GI.geometry(g1), g2) -within(::Any, g1, t2::GI.FeatureTrait, g2)::Bool = within(g1, geometry(g2)) +within(::Any, g1, t2::GI.FeatureTrait, g2)::Bool = within(g1, GI.geometry(g2)) +# Points in geometries within(::GI.PointTrait, g1, ::GI.LineStringTrait, g2)::Bool = point_on_line(g1, g2; ignore_end_vertices=true) +within(::GI.PointTrait, g1, ::GI.LinearRingTrait, g2)::Bool = point_on_line(g1, g2; ignore_end_vertices=true) within(::GI.PointTrait, g1, ::GI.PolygonTrait, g2)::Bool = point_in_polygon(g1, g2; ignore_boundary=true) +# Lines in geometries +within(::GI.LineStringTrait, g1, ::GI.LineStringTrait, g2)::Bool = line_on_line(g1, g2) +within(::GI.LineStringTrait, g1, ::GI.LinearRingTrait, g2)::Bool = line_on_line(g1, g2) within(::GI.LineStringTrait, g1, ::GI.PolygonTrait, g2)::Bool = line_in_polygon(g1, g2) -within(::GI.LineStringTrait, g1, ::GI.LineStringTrait, g2)::Bool = line_on_line(g1, g2) +# Polygons within geometries within(::GI.PolygonTrait, g1, ::GI.PolygonTrait, g2)::Bool = polygon_in_polygon(g1, g2) + +# Everything not specified +# TODO: Add multipolygons +within(::GI.AbstractTrait, g1, ::GI.AbstractCurveTrait, g2)::Bool = false \ No newline at end of file diff --git a/test/methods/bools.jl b/test/methods/bools.jl index b7650cb87..791f0598e 100644 --- a/test/methods/bools.jl +++ b/test/methods/bools.jl @@ -83,7 +83,7 @@ import GeometryOps as GO line8 = GI.LineString([(124.584961, -12.768946), (126.738281, -17.224758)]) line9 = GI.LineString([(123.354492, -15.961329), (127.22168, -14.008696)]) - @test all(GO.line_intersection(line8, line9)[1] .≈ (125.583754, -14.835723)) + @test all(GO.intersection(line8, line9)[1] .≈ (125.583754, -14.835723)) line10 = GI.LineString([ (142.03125, -11.695273), @@ -105,7 +105,7 @@ import GeometryOps as GO (132.890625, -7.754537), ]) - points = GO.line_intersection(line10, line11) + points = GO.intersection(line10, line11) @test all(points[1] .≈ (119.832884, -19.58857)) @test all(points[2] .≈ (132.808697, -11.6309378)) @@ -128,7 +128,7 @@ import GeometryOps as GO (-53.34136962890625, 28.430052892335723), (-53.57208251953125, 28.287451910503744), ]]) - @test GO.overlaps(pl3, pl4) == false + @test GO.overlaps(pl3, pl4) == true # this was false before... why? mp1 = GI.MultiPoint([ (-36.05712890625, 26.480407161007275), diff --git a/test/methods/intersects.jl b/test/methods/intersects.jl index 33ece90cf..f3d35c68f 100644 --- a/test/methods/intersects.jl +++ b/test/methods/intersects.jl @@ -67,38 +67,92 @@ ) # Line strings far apart so extents don't overlap + l1 = LG.LineString([[100.0, 0.0], [101.0, 0.0], [103.0, 0.0]]) + l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) + @test !GO.intersects(l1, l2; meets = 0) + @test !GO.intersects(l1, l2; meets = 1) + @test isnothing(GO.intersection(l1, l2)) # Line strings close together that don't overlap - - # Line string with empty line string + l1 = LG.LineString([[3.0, 0.25], [5.0, 0.25], [7.0, 0.25]]) + l2 = LG.LineString([[0.0, 0.0], [5.0, 10.0], [10.0, 0.0]]) + @test !GO.intersects(l1, l2; meets = 0) + @test !GO.intersects(l1, l2; meets = 1) + @test isempty(GO.intersection(l1, l2)) # Closed linear ring with open line string + r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) + l2 = LG.LineString([[0.0, -2.0], [12.0, 10.0],]) + @test GO.intersects(r1, l2; meets = 0) + @test GO.intersects(r1, l2; meets = 1) + go_inter = GO.intersection(r1, l2) + @test length(go_inter) == 2 + lg_inter = LG.intersection(r1, l2) + @test issetequal( + Set(go_inter), + Set(GO._tuple_point.(GI.getpoint(lg_inter))) + ) # Closed linear ring with closed linear ring - - # @test issetequal( - # Subzero.intersect_lines(l1, l2), - # Set([(0.5, -0.0), (1.5, 0), (2.5, -0.0)]), - # ) - # l2 = [[[10., 10]]] - # @test issetequal( - # Subzero.intersect_lines(l1, l2), - # Set{Tuple{Float64, Float64}}(), - # ) - - + r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) + r2 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) + @test GO.intersects(r1, r2; meets = 0) + @test GO.intersects(r1, r2; meets = 1) + go_inter = GO.intersection(r1, r2) + @test length(go_inter) == 2 + lg_inter = LG.intersection(r1, r2) + @test issetequal( + Set(go_inter), + Set(GO._tuple_point.(GI.getpoint(lg_inter))) + ) end @testset "Polygons" begin # Two polygons that intersect + p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) + p2 = LG.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]]) + @test GO.intersects(p1, p2; meets = 0) + @test GO.intersects(p1, p2; meets = 1) + @test all(GO.intersection_points(p1, p2) .== [(6.5, 3.5), (6.5, -3.5)]) # Two polygons that don't intersect + p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) + p2 = LG.Polygon([[[13.0, 0.0], [18.0, 5.0], [23.0, 0.0], [18.0, -5.0], [13.0, 0.0]]]) + @test !GO.intersects(p1, p2; meets = 0) + @test !GO.intersects(p1, p2; meets = 1) + @test isnothing(GO.intersection_points(p1, p2)) # Polygon that intersects with linestring - + p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) + l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) + @test GO.intersects(p1, l2; meets = 0) + @test GO.intersects(p1, l2; meets = 1) + GO.intersection_points(p1, l2) + @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (10.0, 0.0)]) + + # Polygon with a hole, line through polygon and hole + p1 = LG.Polygon([ + [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]], + [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] + ]) + l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) + @test GO.intersects(p1, l2; meets = 0) + @test GO.intersects(p1, l2; meets = 1) + @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (2.0, 0.0), (3.0, 0.0), (10.0, 0.0)]) + + # Polygon with a hole, line only within the hole + p1 = LG.Polygon([ + [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]], + [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] + ]) + l2 = LG.LineString([[2.25, 0.0], [2.75, 0.0]]) + @test !GO.intersects(p1, l2; meets = 0) + @test !GO.intersects(p1, l2; meets = 1) + @test isempty(GO.intersection_points(p1, l2)) end @testset "MultiPolygons" begin + # TODO: Add these tests # Multi-polygon and polygon that intersect # Multi-polygon and polygon that don't intersect From a7a73671c12735f800c091233178e798ea57a470 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 3 Oct 2023 23:20:06 -0700 Subject: [PATCH 04/16] Add comments to point_in_poly --- src/methods/bools.jl | 58 +++++++++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/src/methods/bools.jl b/src/methods/bools.jl index ba5c4068d..fd6cffa6f 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -265,63 +265,66 @@ function point_in_polygon( end # Then check the point is inside the exterior ring - point_in_polygon(point, GI.getexterior(poly); ignore_boundary, check_extent=false) || return false + point_in_polygon( + point,GI.getexterior(poly); + ignore_boundary, check_extent=false, + ) || return false # Finally make sure the point is not in any of the holes, # flipping the boundary condition for ring in GI.gethole(poly) - point_in_polygon(point, ring; ignore_boundary=!ignore_boundary) && return false + point_in_polygon( + point, ring; + ignore_boundary=!ignore_boundary, + ) && return false end return true end + function point_in_polygon( ::PointTrait, pt, ::Union{LineStringTrait,LinearRingTrait}, ring; ignore_boundary::Bool=false, check_extent::Bool=false, )::Bool + x, y = GI.x(pt), GI.y(pt) # Cheaply check that the point is inside the ring extent if check_extent point_in_extent(point, GI.extent(ring)) || return false end - # Then check the point is inside the ring inside = false n = GI.npoint(ring) p_start = GI.getpoint(ring, 1) p_end = GI.getpoint(ring, n) - - # Handle closed on non-closed rings - l = if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end) - l = n - 1 - else - n + # Handle closed vs opne rings + if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end) + n -= 1 end - # Loop over all points in the ring - for i in 1:l - 1 - j = i + 1 - + for i in 1:(n - 1) + # First point on edge p_i = GI.getpoint(ring, i) - p_j = GI.getpoint(ring, j) - xi = GI.x(p_i) - yi = GI.y(p_i) - xj = GI.x(p_j) - yj = GI.y(p_j) - - on_boundary = (GI.y(pt) * (xi - xj) + yi * (xj - GI.x(pt)) + yj * (GI.x(pt) - xi) == 0) && - ((xi - GI.x(pt)) * (xj - GI.x(pt)) <= 0) && ((yi - GI.y(pt)) * (yj - GI.y(pt)) <= 0) - + xi, yi = GI.x(p_i), GI.y(p_i) + # Second point on edge (j = i + 1) + p_j = GI.getpoint(ring, i + 1) + xj, yj = GI.x(p_j), GI.y(p_j) + # Check if point is on the ring boundary + on_boundary = ( # vertex to point has same slope as edge + yi * (xj - x) + yj * (x - xi) == y * (xj - xi) && + (xi - x) * (xj - x) <= 0 && # x is between xi and xj + (yi - y) * (yj - y) <= 0 # y is between yi and yj + ) on_boundary && return !ignore_boundary - - intersects = ((yi > GI.y(pt)) !== (yj > GI.y(pt))) && - (GI.x(pt) < (xj - xi) * (GI.y(pt) - yi) / (yj - yi) + xi) - + # Check if ray from point passes through edge + intersects = ( + (yi > y) !== (yj > y) && + (x < (xj - xi) * (y - yi) / (yj - yi) + xi) + ) if intersects inside = !inside end end - return inside end @@ -341,6 +344,7 @@ function line_on_line(t1::GI.AbstractCurveTrait, line1, t2::AbstractCurveTrait, end line_in_polygon(line, poly) = line_in_polygon(trait(line), line, trait(poly), poly) + function line_in_polygon( ::AbstractCurveTrait, line, ::Union{AbstractPolygonTrait,LinearRingTrait}, poly From b99e37d62fa7a437af77b108c1f8eb434e95a834 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 4 Oct 2023 12:32:17 -0700 Subject: [PATCH 05/16] Remove CairoMakie --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0fe3b53d8..f6787b264 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ authors = ["Anshul Singhvi and contributors"] version = "0.0.1-DEV" [deps] -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" From 319bd884a2102241e66888286ecb00ac4e507cf8 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 18 Oct 2023 10:17:29 -0700 Subject: [PATCH 06/16] Update equals and overlaps --- src/GeometryOps.jl | 1 + src/methods/bools.jl | 48 ++++---- src/methods/centroid.jl | 2 +- src/methods/crosses.jl | 2 +- src/methods/equals.jl | 192 +++++++++++++++++++++++++++++++ src/methods/intersects.jl | 65 ++++++----- src/methods/overlaps.jl | 228 +++++++++++++++++++++++++++++++++---- test/methods/bools.jl | 34 ------ test/methods/equals.jl | 104 +++++++++++++++++ test/methods/intersects.jl | 49 +++----- test/methods/overlaps.jl | 105 +++++++++++++++++ test/runtests.jl | 2 + 12 files changed, 682 insertions(+), 150 deletions(-) create mode 100644 src/methods/equals.jl create mode 100644 test/methods/equals.jl create mode 100644 test/methods/overlaps.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index ea19f3b31..9e19dd553 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -31,6 +31,7 @@ include("methods/overlaps.jl") include("methods/within.jl") include("methods/polygonize.jl") include("methods/barycentric.jl") +include("methods/equals.jl") include("transformations/flip.jl") include("transformations/simplify.jl") diff --git a/src/methods/bools.jl b/src/methods/bools.jl index fd6cffa6f..30b8716e1 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -11,7 +11,8 @@ export line_on_line, line_in_polygon, polygon_in_polygon """ isclockwise(line::Union{LineString, Vector{Position}})::Bool -Take a ring and return true or false whether or not the ring is clockwise or counter-clockwise. +Take a ring and return true or false whether or not the ring is clockwise or +counter-clockwise. ## Example @@ -26,6 +27,7 @@ true ``` """ isclockwise(geom)::Bool = isclockwise(GI.trait(geom), geom) + function isclockwise(::AbstractCurveTrait, line)::Bool sum = 0.0 prev = GI.getpoint(line, 1) @@ -88,30 +90,6 @@ function isconcave(poly)::Bool return false end -equals(geo1, geo2) = _equals(trait(geo1), geo1, trait(geo2), geo2) - -_equals(::T, geo1, ::T, geo2) where T = error("Cant compare $T yet") -function _equals(::T, p1, ::T, p2) where {T<:PointTrait} - GI.ncoord(p1) == GI.ncoord(p2) || return false - GI.x(p1) == GI.x(p2) || return false - GI.y(p1) == GI.y(p2) || return false - if GI.is3d(p1) - GI.z(p1) == GI.z(p2) || return false - end - return true -end -function _equals(::T, l1, ::T, l2) where {T<:AbstractCurveTrait} - # Check line lengths match - GI.npoint(l1) == GI.npoint(l2) || return false - - # Then check all points are the same - for (p1, p2) in zip(GI.getpoint(l1), GI.getpoint(l2)) - equals(p1, p2) || return false - end - return true -end -_equals(t1, geo1, t2, geo2) = false - # """ # isparallel(line1::LineString, line2::LineString)::Bool @@ -193,6 +171,26 @@ function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool return false end +function point_on_seg(point, start, stop) + # Parse out points + x, y = GI.x(point), GI.y(point) + x1, y1 = GI.x(start), GI.y(start) + x2, y2 = GI.x(stop), GI.y(stop) + Δxl = x2 - x1 + Δyl = y2 - y1 + # Determine if point is on segment + cross = (x - x1) * Δyl - (y - y1) * Δxl + if cross == 0 # point is on line extending to infinity + # is line between endpoints + if abs(Δxl) >= abs(Δyl) # is line between endpoints + return Δxl > 0 ? x1 <= x <= x2 : x2 <= x <= x1 + else + return Δyl > 0 ? y1 <= y <= y2 : y2 <= y <= y1 + end + end + return false +end + function point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)::Bool x, y = GI.x(point), GI.y(point) x1, y1 = GI.x(start), GI.y(start) diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl index 03dbc6798..6a15808d7 100644 --- a/src/methods/centroid.jl +++ b/src/methods/centroid.jl @@ -216,7 +216,7 @@ function centroid_and_area(::GI.MultiPolygonTrait, geom) xcentroid *= area ycentroid *= area # Loop over any polygons within the multipolygon - for i in 2:GI.ngeom(geom) #poly in GI.getpolygon(geom) + for i in 2:GI.ngeom(geom) # Polygon centroid and area (xpoly, ypoly), poly_area = centroid_and_area(GI.getpolygon(geom, i)) # Accumulate the area component into `area` diff --git a/src/methods/crosses.jl b/src/methods/crosses.jl index 3aa62d62e..f8a580db0 100644 --- a/src/methods/crosses.jl +++ b/src/methods/crosses.jl @@ -55,7 +55,7 @@ end function line_crosses_line(line1, line2) np2 = GI.npoint(line2) - if intersects(line1, line2; meets=MEETS_CLOSED) + if intersects(line1, line2) for i in 1:GI.npoint(line1) - 1 for j in 1:GI.npoint(line2) - 1 exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both diff --git a/src/methods/equals.jl b/src/methods/equals.jl new file mode 100644 index 000000000..568256845 --- /dev/null +++ b/src/methods/equals.jl @@ -0,0 +1,192 @@ +# # Equals + +export equals + +#= +## What is equals? + +The equals function checks if two geometries are equal. They are equal if they +share the same set of points and edges. + +To provide an example, consider these two lines: +```@example cshape +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (0.0, 10.0)]) +l2 = GI.LineString([(0.0, -10.0), (0.0, 3.0)]) +f, a, p = lines(GI.getpoint(l1), color = :blue) +scatter!(GI.getpoint(l1), color = :blue) +lines!(GI.getpoint(l2), color = :orange) +scatter!(GI.getpoint(l2), color = :orange) +``` +We can see that the two lines do not share a commen set of points and edges in +the plot, so they are not equal: +```@example cshape +equals(l1, l2) # returns false +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. This is also used in the +implementation, since it's a lot less work! + +Note that while we need the same set of points and edges, they don't need to be +provided in the same order for polygons. For for example, we need the same set +points for two multipoints to be equal, but they don't have to be saved in the +same order. This requires checking every point against every other point in the +two geometries we are comparing. +=# + +""" + equals(geom1, geom2)::Bool + +Compare two Geometries return true if they are the same geometry. + +## Examples +```jldoctest +import GeometryOps as GO, GeoInterface as GI +poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]]) +poly2 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]]) + +GO.equals(poly1, poly2) +# output +true +``` +""" +equals(geom_a, geom_b) = equals( + GI.trait(geom_a), geom_a, + GI.trait(geom_b), geom_b, +) + +""" + equals(::T, geom_a, ::T, geom_b)::Bool + +Two geometries of the same type, which don't have a equals function to dispatch +off of should throw an error. +""" +equals(::T, geom_a, ::T, geom_b) where T = error("Cant compare $T yet") + +""" + equals(trait_a, geom_a, trait_b, geom_b) + +Two geometries which are not of the same type cannot be equal so they always +return false. +""" +equals(trait_a, geom_a, trait_b, geom_b) = false + +""" + equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)::Bool + +Two points are the same if they have the same x and y (and z if 3D) coordinates. +""" +function equals(::GI.PointTrait, p1, ::GI.PointTrait, p2) + GI.ncoord(p1) == GI.ncoord(p2) || return false + GI.x(p1) == GI.x(p2) || return false + GI.y(p1) == GI.y(p2) || return false + if GI.is3d(p1) + GI.z(p1) == GI.z(p2) || return false + end + return true +end + +""" + equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)::Bool + +Two multipoints are equal if they share the same set of points. +""" +function equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2) + GI.npoint(mp1) == GI.npoint(mp2) || return false + for p1 in GI.getpoint(mp1) + has_match = false # if point has a matching point in other multipoint + for p2 in GI.getpoint(mp2) + if equals(p1, p2) + has_match = true + break + end + end + has_match || return false # if no matching point, can't be equal + end + return true # all points had a match +end + +""" + equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} ::Bool + +Two curves are equal if they share the same set of points going around the +curve. +""" +function equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} + # Check line lengths match + n1 = GI.npoint(l1) + n2 = GI.npoint(l2) + # TODO: do we need to account for repeated last point?? + n1 == n2 || return false + + # Find first matching point if it exists + p1 = GI.getpoint(l1, 1) + offset = findfirst(p2 -> equals(p1, p2), GI.getpoint(l2)) + isnothing(offset) && return false + offset -= 1 + + # Then check all points are the same wrapping around line + for i in 1:n1 + pi = GI.getpoint(l1, i) + j = i + offset + j = j <= n1 ? j : (j - n1) + pj = GI.getpoint(l2, j) + equals(pi, pj) || return false + end + return true +end + +""" + equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool + +Two polygons are equal if they share the same exterior edge and holes. +""" +function equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b) + # Check if exterior is equal + equals(GI.getexterior(geom_a), GI.getexterior(geom_b)) || return false + # Check if number of holes are equal + GI.nhole(geom_a) == GI.nhole(geom_b) || return false + # Check if holes are equal + for ihole in GI.gethole(geom_a) + has_match = false + for jhole in GI.gethole(geom_b) + if equals(ihole, jhole) + has_match = true + break + end + end + has_match || return false + end + return true +end + +""" + equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool + +Two multipolygons are equal if they share the same set of polygons. +""" +function equals(::GI.MultiPolygonTrait, geom_a, ::GI.MultiPolygonTrait, geom_b) + # Check if same number of polygons + GI.npolygon(geom_a) == GI.npolygon(geom_b) || return false + # Check if each polygon has a matching polygon + for poly_a in GI.getpolygon(geom_a) + has_match = false + for poly_b in GI.getpolygon(geom_b) + if equals(poly_a, poly_b) + has_match = true + break + end + end + has_match || return false + end + return true +end \ No newline at end of file diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl index 78f5784f1..2efaf1b78 100644 --- a/src/methods/intersects.jl +++ b/src/methods/intersects.jl @@ -52,16 +52,10 @@ intersect and _intersection_point which determines the intersection point between two line segments. =# -const MEETS_CLOSED = 0 -const MEETS_OPEN = 1 - """ - intersects(geom1, geom2; kw...)::Bool + intersects(geom1, geom2)::Bool Check if two geometries intersect, returning true if so and false otherwise. -Takes in a Int keyword meets, which can either be MEETS_OPEN (1), meaning that -only intersections through open edges where edge endpoints are not included are -recorded, versus MEETS_CLOSED (0) where edge endpoints are included. ## Example @@ -76,73 +70,78 @@ GO.intersects(line1, line2) true ``` """ -intersects(geom1, geom2; kw...) = intersects( +intersects(geom1, geom2) = intersects( GI.trait(geom1), geom1, GI.trait(geom2), - geom2; - kw... + geom2 ) """ - intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets = MEETS_OPEN)::Bool + intersects(::GI.LineTrait, a, ::GI.LineTrait, b)::Bool -Returns true if two line segments intersect and false otherwise. Line segment -endpoints are excluded in check if `meets = MEETS_OPEN` (1) and included if -`meets = MEETS_CLOSED` (0). +Returns true if two line segments intersect and false otherwise. """ -function intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets = MEETS_OPEN) +function intersects(::GI.LineTrait, a, ::GI.LineTrait, b) a1 = _tuple_point(GI.getpoint(a, 1)) a2 = _tuple_point(GI.getpoint(a, 2)) b1 = _tuple_point(GI.getpoint(b, 1)) b2 = _tuple_point(GI.getpoint(b, 2)) meet_type = ExactPredicates.meet(a1, a2, b1, b2) - return meet_type == MEETS_OPEN || meet_type == meets + return meet_type == 0 || meet_type == 1 end """ - intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...)::Bool + intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b)::Bool Returns true if two geometries intersect with one another and false -otherwise. For all geometries but lines, conver the geometry to a list of edges +otherwise. For all geometries but lines, convert the geometry to a list of edges and cross compare the edges for intersections. """ function intersects( - trait_a::GI.AbstractTrait, a, - trait_b::GI.AbstractTrait, b; - kw..., -) - edges_a, edges_b = map(sort! ∘ to_edges, (a, b)) - return _line_intersects(edges_a, edges_b; kw...) || - within(trait_a, a, trait_b, b) || within(trait_b, b, trait_a, a) + trait_a::GI.AbstractTrait, a_geom, + trait_b::GI.AbstractTrait, b_geom, +) edges_a, edges_b = map(sort! ∘ to_edges, (a_geom, b_geom)) + return _line_intersects(edges_a, edges_b) || + within(trait_a, a_geom, trait_b, b_geom) || + within(trait_b, b_geom, trait_a, a_geom) end """ _line_intersects( edges_a::Vector{Edge}, - edges_b::Vector{Edge}; - meets = MEETS_OPEN, + edges_b::Vector{Edge} )::Bool Returns true if there is at least one intersection between edges within the -two lists. Line segment endpoints are excluded in check if `meets = MEETS_OPEN` -(1) and included if `meets = MEETS_CLOSED` (0). +two lists of edges. """ function _line_intersects( edges_a::Vector{Edge}, - edges_b::Vector{Edge}; - meets = MEETS_OPEN, + edges_b::Vector{Edge} ) # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false for edge_a in edges_a for edge_b in edges_b - meet_type = ExactPredicates.meet(edge_a..., edge_b...) - (meet_type == MEETS_OPEN || meet_type == meets) && return true + _line_intersects(edge_a, edge_b) && return true end end return false end +""" + _line_intersects( + edge_a::Edge, + edge_b::Edge, + )::Bool + +Returns true if there is at least one intersection between two edges. +""" +function _line_intersects(edge_a::Edge, edge_b::Edge) + meet_type = ExactPredicates.meet(edge_a..., edge_b...) + return meet_type == 0 || meet_type == 1 +end + """ intersection(geom_a, geom_b)::Union{Tuple{::Real, ::Real}, ::Nothing} diff --git a/src/methods/overlaps.jl b/src/methods/overlaps.jl index 6d84f393b..f99b75c9d 100644 --- a/src/methods/overlaps.jl +++ b/src/methods/overlaps.jl @@ -1,17 +1,61 @@ -# # Overlap checks +# # Overlaps export overlaps -# This code checks whether geometries overlap with each other. +#= +## What is overlaps? -# It does not compute the overlap or intersection geometry. +The overlaps function checks if two geometries overlap. Two geometries can only +overlap if they have the same dimension, and if they overlap, but one is not +contained, within, or equal to the other. + +Note that this means it is impossible for a single point to overlap with a +single point and a line only overlaps with another line if only a section of +each line is colinear. + +To provide an example, consider these two lines: +```@example cshape +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (0.0, 10.0)]) +l2 = GI.LineString([(0.0, -10.0), (0.0, 3.0)]) +f, a, p = lines(GI.getpoint(l1), color = :blue) +scatter!(GI.getpoint(l1), color = :blue) +lines!(GI.getpoint(l2), color = :orange) +scatter!(GI.getpoint(l2), color = :orange) +``` +We can see that the two lines overlap in the plot: +```@example cshape +overlap(l1, l2) +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. This is also used in the +implementation, since it's a lot less work! + +Note that that since only elements of the same dimension can overlap, any two +geometries with traits that are of different dimensions autmoatically can +return false. + +For geometries with the same trait dimension, we must make sure that they share +a point, an edge, or area for points, lines, and polygons/multipolygons +respectivly, without being contained. +=# """ overlaps(geom1, geom2)::Bool -Compare two Geometries of the same dimension and return true if their intersection set results in a geometry -different from both but of the same dimension. It applies to Polygon/Polygon, LineString/LineString, -Multipoint/Multipoint, MultiLineString/MultiLineString and MultiPolygon/MultiPolygon. +Compare two Geometries of the same dimension and return true if their +intersection set results in a geometry different from both but of the same +dimension. This means one geometry cannot be within or contain the other and +they cannot be equal ## Examples ```jldoctest @@ -24,28 +68,166 @@ GO.overlaps(poly1, poly2) true ``` """ -overlaps(g1, g2)::Bool = overlaps(trait(g1), g1, trait(g2), g2)::Bool -overlaps(t1::FeatureTrait, g1, t2, g2)::Bool = overlaps(GI.geometry(g1), g2) -overlaps(t1, g1, t2::FeatureTrait, g2)::Bool = overlaps(g1, geometry(g2)) -overlaps(t1::FeatureTrait, g1, t2::FeatureTrait, g2)::Bool = overlaps(geometry(g1), geometry(g2)) -overlaps(::PolygonTrait, mp, ::MultiPolygonTrait, p)::Bool = overlaps(p, mp) -function overlaps(::MultiPointTrait, g1, ::MultiPointTrait, g2)::Bool - for p1 in GI.getpoint(g1) - for p2 in GI.getpoint(g2) - equals(p1, p2) && return true +overlaps(geom1, geom2)::Bool = overlaps( + GI.trait(geom1), + geom1, + GI.trait(geom2), + geom2, +) + +""" + overlaps(::GI.AbstractTrait, geom1, ::GI.AbstractTrait, geom2)::Bool + +For any non-specified pair, all have non-matching dimensions, return false. +""" +overlaps(::GI.AbstractTrait, geom1, ::GI.AbstractTrait, geom2) = false + +""" + overlaps( + ::GI.MultiPointTrait, points1, + ::GI.MultiPointTrait, points2, + )::Bool + +If the multipoints overlap, meaning some, but not all, of the points within the +multipoints are shared, return true. +""" +function overlaps( + ::GI.MultiPointTrait, points1, + ::GI.MultiPointTrait, points2, +) + one_diff = false # assume that all the points are the same + one_same = false # assume that all points are different + for p1 in GI.getpoint(points1) + match_point = false + for p2 in GI.getpoint(points2) + if equals(p1, p2) # Point is shared + one_same = true + match_point = true + break + end + end + one_diff |= !match_point # Point isn't shared + one_same && one_diff && return true + end + return false +end + +""" + overlaps(::GI.LineTrait, line1, ::GI.LineTrait, line)::Bool + +If the lines overlap, meaning that they are colinear but each have one endpoint +outside of the other line, return true. Else false. +""" +overlaps(::GI.LineTrait, line1, ::GI.LineTrait, line) = + _overlaps((a1, a2), (b1, b2)) + +""" + overlaps( + ::Union{GI.LineStringTrait, GI.LinearRing}, line1, + ::Union{GI.LineStringTrait, GI.LinearRing}, line2, + )::Bool + +If the curves overlap, meaning that at least one edge of each curve overlaps, +return true. Else false. +""" +function overlaps( + ::Union{GI.LineStringTrait, GI.LinearRing}, line1, + ::Union{GI.LineStringTrait, GI.LinearRing}, line2, +) + edges_a, edges_b = map(sort! ∘ to_edges, (line1, line2)) + for edge_a in edges_a + for edge_b in edges_b + _overlaps(edge_a, edge_b) && return true end end + return false end -function overlaps(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool - return intersects(g1, g2) + +""" + overlaps( + trait_a::GI.PolygonTrait, poly_a, + trait_b::GI.PolygonTrait, poly_b, + )::Bool + +If the two polygons intersect with one another, but are not equal, return true. +Else false. +""" +function overlaps( + trait_a::GI.PolygonTrait, poly_a, + trait_b::GI.PolygonTrait, poly_b, +) + edges_a, edges_b = map(sort! ∘ to_edges, (poly_a, poly_b)) + return _line_intersects(edges_a, edges_b) && + !equals(trait_a, poly_a, trait_b, poly_b) end -function overlaps(t1::MultiPolygonTrait, mp, t2::PolygonTrait, p1)::Bool - for p2 in GI.getgeom(mp) - overlaps(p1, p2) && return true + +""" + overlaps( + ::GI.PolygonTrait, poly1, + ::GI.MultiPolygonTrait, polys2, + )::Bool + +Return true if polygon overlaps with at least one of the polygons within the +multipolygon. Else false. +""" +function overlaps( + ::GI.PolygonTrait, poly1, + ::GI.MultiPolygonTrait, polys2, +) + for poly2 in GI.getgeom(polys2) + overlaps(poly1, poly2) && return true end + return false end -function overlaps(::MultiPolygonTrait, g1, ::MultiPolygonTrait, g2)::Bool - for p1 in GI.getgeom(g1) - overlaps(PolygonTrait(), mp, PolygonTrait(), p1) && return true + +""" + overlaps( + ::GI.MultiPolygonTrait, polys1, + ::GI.PolygonTrait, poly2, + )::Bool + +Return true if polygon overlaps with at least one of the polygons within the +multipolygon. Else false. +""" +overlaps(trait1::GI.MultiPolygonTrait, polys1, trait2::GI.PolygonTrait, poly2) = + overlaps(trait2, poly2, trait1, polys1) + +""" + overlaps( + ::GI.MultiPolygonTrait, polys1, + ::GI.MultiPolygonTrait, polys2, + )::Bool + +Return true if at least one pair of polygons from multipolygons overlap. Else +false. +""" +function overlaps( + ::GI.MultiPolygonTrait, polys1, + ::GI.MultiPolygonTrait, polys2, +) + for poly1 in GI.getgeom(polys1) + overlaps(poly1, polys2) && return true end + return false +end + +""" + _overlaps( + (a1, a2)::Edge, + (b1, b2)::Edge + )::Bool + +If the edges overlap, meaning that they are colinear but each have one endpoint +outside of the other edge, return true. Else false. +""" +function _overlaps( + (a1, a2)::Edge, + (b1, b2)::Edge +) + # meets in more than one point + on_top = ExactPredicates.meet(a1, a2, b1, b2) == 0 + # one end point is outside of other segment + a_fully_within = point_on_seg(a1, b1, b2) && point_on_seg(a2, b1, b2) + b_fully_within = point_on_seg(b1, a1, a2) && point_on_seg(b2, a1, a2) + return on_top && (!a_fully_within && !b_fully_within) end diff --git a/test/methods/bools.jl b/test/methods/bools.jl index 791f0598e..cb1ff945c 100644 --- a/test/methods/bools.jl +++ b/test/methods/bools.jl @@ -114,38 +114,4 @@ import GeometryOps as GO @test GO.crosses(GI.MultiPoint([(1, 2), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true @test GO.crosses(GI.MultiPoint([(1, 0), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false @test GO.crosses(GI.LineString([(-2, 2), (-4, 2)]), poly7) == false - - pl1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]]) - pl2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]]) - - @test GO.overlaps(pl1, pl2) == true - @test_throws MethodError GO.overlaps(pl1, (1, 1)) - @test_throws MethodError GO.overlaps((1, 1), pl2) - - pl3 = pl4 = GI.Polygon([[ - (-53.57208251953125, 28.287451910503744), - (-53.33038330078125, 28.29228897739706), - (-53.34136962890625, 28.430052892335723), - (-53.57208251953125, 28.287451910503744), - ]]) - @test GO.overlaps(pl3, pl4) == true # this was false before... why? - - mp1 = GI.MultiPoint([ - (-36.05712890625, 26.480407161007275), - (-35.7220458984375, 27.137368359795584), - (-35.13427734375, 26.83387451505858), - (-35.4638671875, 27.254629577800063), - (-35.5462646484375, 26.86328062676624), - (-35.3924560546875, 26.504988828743404) - ]) - mp2 = GI.MultiPoint([ - (-35.4638671875, 27.254629577800063), - (-35.5462646484375, 26.86328062676624), - (-35.3924560546875, 26.504988828743404), - (-35.2001953125, 26.12091815959972), - (-34.9969482421875, 26.455820238459893) - ]) - - @test GO.overlaps(mp1, mp2) == true - @test GO.overlaps(mp1, mp2) == GO.overlaps(mp2, mp1) end diff --git a/test/methods/equals.jl b/test/methods/equals.jl new file mode 100644 index 000000000..a0b60d6cd --- /dev/null +++ b/test/methods/equals.jl @@ -0,0 +1,104 @@ +@testset "Points/MultiPoints" begin + p1 = LG.Point([0.0, 0.0]) + p2 = LG.Point([0.0, 1.0]) + # Same points + @test GO.equals(p1, p1) + @test GO.equals(p2, p2) + # Different points + @test !GO.equals(p1, p2) + + mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) + mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) + # Same points + @test LG.equals(mp1, mp1) + @test LG.equals(mp2, mp2) + # Different points + @test !LG.equals(mp1, mp2) + @test !LG.equals(mp1, p1) +end + +@testset "Lines/Rings" begin + l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) + l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) + # Equal lines + @test LG.equals(l1, l1) + @test LG.equals(l2, l2) + # Different lines + @test !LG.equals(l1, l2) && !LG.equals(l2, l1) + + r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) + r2 = LG.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) + l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) + # Equal rings + @test GO.equals(r1, r1) + @test GO.equals(r2, r2) + # Different rings + @test !GO.equals(r1, r2) && !GO.equals(r2, r1) + # Equal linear ring and line string + @test !GO.equals(r2, l3) # TODO: should these be equal? +end + +@testset "Polygons/MultiPolygons" begin + p1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]]) + p2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]]) + p3 = LG.Polygon( + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ] + ) + p4 = LG.Polygon( + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[16.0, 1.0], [16.0, 11.0], [25.0, 11.0], [25.0, 1.0], [16.0, 1.0]] + ] + ) + p5 = LG.Polygon( + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]], + [[11.0, 1.0], [11.0, 2.0], [12.0, 2.0], [12.0, 1.0], [11.0, 1.0]] + ] + ) + # Equal polygon + @test GO.equals(p1, p1) + @test GO.equals(p2, p2) + # Different polygons + @test !GO.equals(p1, p2) + # Equal polygons with holes + @test GO.equals(p3, p3) + # Same exterior, different hole + @test !GO.equals(p3, p4) + # Same exterior and first hole, has an extra hole + @test !GO.equals(p3, p5) + + p3 = GI.Polygon( + [[ + [-53.57208251953125, 28.287451910503744], + [-53.33038330078125, 28.29228897739706], + [-53.34136962890625, 28.430052892335723], + [-53.57208251953125, 28.287451910503744], + ]] + ) + # Complex polygon + @test GO.equals(p3, p3) + + m1 = LG.MultiPolygon([ + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ] + ]) + m2 = LG.MultiPolygon([ + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ], + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]] + ]) + # Equal multipolygon + @test GO.equals(m1, m1) + # Equal multipolygon with different order + @test GO.equals(m1, m2) +end \ No newline at end of file diff --git a/test/methods/intersects.jl b/test/methods/intersects.jl index f3d35c68f..4251d45a8 100644 --- a/test/methods/intersects.jl +++ b/test/methods/intersects.jl @@ -4,29 +4,25 @@ # Test for parallel lines l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) l2 = GI.Line([(0.0, 1.0), (2.5, 1.0)]) - @test !GO.intersects(l1, l2; meets = 0) - @test !GO.intersects(l1, l2; meets = 1) + @test !GO.intersects(l1, l2) @test isnothing(GO.intersection(l1, l2)) # Test for non-parallel lines that don't intersect l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) l2 = GI.Line([(2.0, -3.0), (3.0, 0.0)]) - @test !GO.intersects(l1, l2; meets = 0) - @test !GO.intersects(l1, l2; meets = 1) + @test !GO.intersects(l1, l2) @test isnothing(GO.intersection(l1, l2)) # Test for lines only touching at endpoint l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) l2 = GI.Line([(2.0, -3.0), (2.5, 0.0)]) - @test GO.intersects(l1, l2; meets = 0) - @test !GO.intersects(l1, l2; meets = 1) + @test GO.intersects(l1, l2) @test all(GO.intersection(l1, l2) .≈ (2.5, 0.0)) # Test for lines that intersect in the middle l1 = GI.Line([(0.0, 0.0), (5.0, 5.0)]) l2 = GI.Line([(0.0, 5.0), (5.0, 0.0)]) - @test GO.intersects(l1, l2; meets = 0) - @test GO.intersects(l1, l2; meets = 1) + @test GO.intersects(l1, l2) @test all(GO.intersection(l1, l2) .≈ (2.5, 2.5)) # Line string test intersects ---------------------------------------------- @@ -34,8 +30,7 @@ # Single element line strings crossing over each other l1 = LG.LineString([[5.5, 7.2], [11.2, 12.7]]) l2 = LG.LineString([[4.3, 13.3], [9.6, 8.1]]) - @test GO.intersects(l1, l2; meets = 0) - @test GO.intersects(l1, l2; meets = 1) + @test GO.intersects(l1, l2) go_inter = GO.intersection(l1, l2) lg_inter = LG.intersection(l1, l2) @test go_inter[1][1] .≈ GI.x(lg_inter) @@ -44,9 +39,7 @@ # Multi-element line strings crossing over on vertex l1 = LG.LineString([[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]]) l2 = LG.LineString([[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]]) - @test GO.intersects(l1, l2; meets = 0) - # TODO: Do we want this to be false? It is vertex of segment, not of whole line string - @test !GO.intersects(l1, l2; meets = 1) + @test GO.intersects(l1, l2) go_inter = GO.intersection(l1, l2) @test length(go_inter) == 1 lg_inter = LG.intersection(l1, l2) @@ -56,8 +49,7 @@ # Multi-element line strings crossing over with multiple intersections l1 = LG.LineString([[0.0, -1.0], [1.0, 1.0], [2.0, -1.0], [3.0, 1.0]]) l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) - @test GO.intersects(l1, l2; meets = 0) - @test GO.intersects(l1, l2; meets = 1) + @test GO.intersects(l1, l2) go_inter = GO.intersection(l1, l2) @test length(go_inter) == 3 lg_inter = LG.intersection(l1, l2) @@ -69,22 +61,19 @@ # Line strings far apart so extents don't overlap l1 = LG.LineString([[100.0, 0.0], [101.0, 0.0], [103.0, 0.0]]) l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) - @test !GO.intersects(l1, l2; meets = 0) - @test !GO.intersects(l1, l2; meets = 1) + @test !GO.intersects(l1, l2) @test isnothing(GO.intersection(l1, l2)) # Line strings close together that don't overlap l1 = LG.LineString([[3.0, 0.25], [5.0, 0.25], [7.0, 0.25]]) l2 = LG.LineString([[0.0, 0.0], [5.0, 10.0], [10.0, 0.0]]) - @test !GO.intersects(l1, l2; meets = 0) - @test !GO.intersects(l1, l2; meets = 1) + @test !GO.intersects(l1, l2) @test isempty(GO.intersection(l1, l2)) # Closed linear ring with open line string r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) l2 = LG.LineString([[0.0, -2.0], [12.0, 10.0],]) - @test GO.intersects(r1, l2; meets = 0) - @test GO.intersects(r1, l2; meets = 1) + @test GO.intersects(r1, l2) go_inter = GO.intersection(r1, l2) @test length(go_inter) == 2 lg_inter = LG.intersection(r1, l2) @@ -96,8 +85,7 @@ # Closed linear ring with closed linear ring r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) r2 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) - @test GO.intersects(r1, r2; meets = 0) - @test GO.intersects(r1, r2; meets = 1) + @test GO.intersects(r1, r2) go_inter = GO.intersection(r1, r2) @test length(go_inter) == 2 lg_inter = LG.intersection(r1, r2) @@ -111,22 +99,19 @@ end # Two polygons that intersect p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) p2 = LG.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]]) - @test GO.intersects(p1, p2; meets = 0) - @test GO.intersects(p1, p2; meets = 1) + @test GO.intersects(p1, p2) @test all(GO.intersection_points(p1, p2) .== [(6.5, 3.5), (6.5, -3.5)]) # Two polygons that don't intersect p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) p2 = LG.Polygon([[[13.0, 0.0], [18.0, 5.0], [23.0, 0.0], [18.0, -5.0], [13.0, 0.0]]]) - @test !GO.intersects(p1, p2; meets = 0) - @test !GO.intersects(p1, p2; meets = 1) + @test !GO.intersects(p1, p2) @test isnothing(GO.intersection_points(p1, p2)) # Polygon that intersects with linestring p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) - @test GO.intersects(p1, l2; meets = 0) - @test GO.intersects(p1, l2; meets = 1) + @test GO.intersects(p1, l2) GO.intersection_points(p1, l2) @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (10.0, 0.0)]) @@ -136,8 +121,7 @@ end [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] ]) l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) - @test GO.intersects(p1, l2; meets = 0) - @test GO.intersects(p1, l2; meets = 1) + @test GO.intersects(p1, l2) @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (2.0, 0.0), (3.0, 0.0), (10.0, 0.0)]) # Polygon with a hole, line only within the hole @@ -146,8 +130,7 @@ end [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] ]) l2 = LG.LineString([[2.25, 0.0], [2.75, 0.0]]) - @test !GO.intersects(p1, l2; meets = 0) - @test !GO.intersects(p1, l2; meets = 1) + @test !GO.intersects(p1, l2) @test isempty(GO.intersection_points(p1, l2)) end diff --git a/test/methods/overlaps.jl b/test/methods/overlaps.jl new file mode 100644 index 000000000..5123fd3f7 --- /dev/null +++ b/test/methods/overlaps.jl @@ -0,0 +1,105 @@ +@testset "Points/MultiPoints" begin + p1 = LG.Point([0.0, 0.0]) + p2 = LG.Point([0.0, 1.0]) + # Two points can't overlap + @test GO.overlaps(p1, p1) == LG.overlaps(p1, p2) + + mp1 = LG.MultiPoint([[0.0, 1.0], [4.0, 4.0]]) + mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) + mp3 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) + # No shared points, doesn't overlap + @test GO.overlaps(p1, mp1) == LG.overlaps(p1, mp1) + # One shared point, does overlap + @test GO.overlaps(p2, mp1) == LG.overlaps(p2, mp1) + # All shared points, doesn't overlap + @test GO.overlaps(mp1, mp1) == LG.overlaps(mp1, mp1) + # Not all shared points, overlaps + @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2) + # One set of points entirely inside other set, doesn't overlap + @test GO.overlaps(mp2, mp3) == LG.overlaps(mp2, mp3) + # Not all points shared, overlaps + @test GO.overlaps(mp1, mp3) == LG.overlaps(mp1, mp3) + + mp1 = LG.MultiPoint([ + [-36.05712890625, 26.480407161007275], + [-35.7220458984375, 27.137368359795584], + [-35.13427734375, 26.83387451505858], + [-35.4638671875, 27.254629577800063], + [-35.5462646484375, 26.86328062676624], + [-35.3924560546875, 26.504988828743404], + ]) + mp2 = GI.MultiPoint([ + [-35.4638671875, 27.254629577800063], + [-35.5462646484375, 26.86328062676624], + [-35.3924560546875, 26.504988828743404], + [-35.2001953125, 26.12091815959972], + [-34.9969482421875, 26.455820238459893], + ]) + # Some shared points, overlaps + @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2) + @test GO.overlaps(mp1, mp2) == GO.overlaps(mp2, mp1) +end + +@testset "Lines/Rings" begin + l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) + l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) + l3 = LG.LineString([[0.0, -10.0], [0.0, 3.0]]) + l4 = LG.LineString([[5.0, -5.0], [5.0, 5.0]]) + # Line can't overlap with itself + @test GO.overlaps(l1, l1) == LG.overlaps(l1, l1) + # Line completely within other line doesn't overlap + @test GO.overlaps(l1, l2) == GO.overlaps(l2, l1) == LG.overlaps(l1, l2) + # Overlapping lines + @test GO.overlaps(l1, l3) == GO.overlaps(l3, l1) == LG.overlaps(l1, l3) + # Lines that don't touch + @test GO.overlaps(l1, l4) == LG.overlaps(l1, l4) + # Linear rings that intersect but don't overlap + r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]) + r2 = LG.LinearRing([[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]) + @test LG.overlaps(r1, r2) == LG.overlaps(r1, r2) +end + +@testset "Polygons/MultiPolygons" begin + p1 = LG.Polygon([[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]]) + p2 = LG.Polygon([ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ]) + # Test basic polygons that don't overlap + @test GO.overlaps(p1, p2) == LG.overlaps(p1, p2) + + p3 = LG.Polygon([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]]) + # Test basic polygons that overlap + @test GO.overlaps(p1, p3) == LG.overlaps(p1, p3) + + p4 = LG.Polygon([[[20.0, 5.0], [20.0, 10.0], [18.0, 10.0], [18.0, 5.0], [20.0, 5.0]]]) + # Test one polygon within the other + @test GO.overlaps(p2, p4) == GO.overlaps(p4, p2) == LG.overlaps(p2, p4) + + # @test_throws MethodError GO.overlaps(pl1, (1, 1)) # I think these should be false + # @test_throws MethodError GO.overlaps((1, 1), pl2) + + p5 = LG.Polygon( + [[ + [-53.57208251953125, 28.287451910503744], + [-53.33038330078125, 28.29228897739706], + [-53.34136352890625, 28.430052892335723], + [-53.57208251953125, 28.287451910503744], + ]] + ) + # Test equal polygons + @test GO.overlaps(p5, p5) == LG.overlaps(p5, p5) + + # Test multipolygons + m1 = LG.MultiPolygon([ + [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], + [ + [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]], + [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]] + ] + ]) + # Test polygon that overlaps with multipolygon + @test GO.overlaps(m1, p3) == LG.overlaps(m1, p3) + # Test polygon in hole of multipolygon, doesn't overlap + @test GO.overlaps(m1, p4) == LG.overlaps(m1, p4) +end diff --git a/test/runtests.jl b/test/runtests.jl index 7c96de785..ee2065017 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,8 +18,10 @@ const GO = GeometryOps @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end @testset "Bools" begin include("methods/bools.jl") end @testset "Centroid" begin include("methods/centroid.jl") end + @testset "Equals" begin include("methods/equals.jl") end @testset "Intersect" begin include("methods/intersects.jl") end @testset "Signed Area" begin include("methods/signed_area.jl") end + @testset "Overlaps" begin include("methods/overlaps.jl") end # Transformations @testset "Reproject" begin include("transformations/reproject.jl") end @testset "Flip" begin include("transformations/flip.jl") end From 0b9799d4309ec01029abbae799fe38bc89a44e11 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 18 Oct 2023 13:03:39 -0700 Subject: [PATCH 07/16] Remove use of findfirst for 1.6 compat --- src/methods/equals.jl | 9 +++++++-- test/methods/overlaps.jl | 5 ++--- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/methods/equals.jl b/src/methods/equals.jl index 568256845..274866631 100644 --- a/src/methods/equals.jl +++ b/src/methods/equals.jl @@ -130,9 +130,14 @@ function equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} # Find first matching point if it exists p1 = GI.getpoint(l1, 1) - offset = findfirst(p2 -> equals(p1, p2), GI.getpoint(l2)) + offset = nothing + for i in 1:n2 + if equals(p1, GI.getpoint(l2, i)) + offset = i - 1 + break + end + end isnothing(offset) && return false - offset -= 1 # Then check all points are the same wrapping around line for i in 1:n1 diff --git a/test/methods/overlaps.jl b/test/methods/overlaps.jl index 5123fd3f7..0c2dab96d 100644 --- a/test/methods/overlaps.jl +++ b/test/methods/overlaps.jl @@ -67,6 +67,8 @@ end ]) # Test basic polygons that don't overlap @test GO.overlaps(p1, p2) == LG.overlaps(p1, p2) + @test !GO.overlaps(p1, (1, 1)) + @test !GO.overlaps((1, 1), p2) p3 = LG.Polygon([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]]) # Test basic polygons that overlap @@ -76,9 +78,6 @@ end # Test one polygon within the other @test GO.overlaps(p2, p4) == GO.overlaps(p4, p2) == LG.overlaps(p2, p4) - # @test_throws MethodError GO.overlaps(pl1, (1, 1)) # I think these should be false - # @test_throws MethodError GO.overlaps((1, 1), pl2) - p5 = LG.Polygon( [[ [-53.57208251953125, 28.287451910503744], From 90fff1c211be896953a0463be90f6c522ee8aef3 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Fri, 20 Oct 2023 00:27:01 -0700 Subject: [PATCH 08/16] Updated geom, multi-geom equality --- src/methods/equals.jl | 26 +++++++++++++++++--- src/try.jl | 7 ++++++ test/methods/equals.jl | 56 ++++++++++++++++++++++++------------------ 3 files changed, 62 insertions(+), 27 deletions(-) create mode 100644 src/try.jl diff --git a/src/methods/equals.jl b/src/methods/equals.jl index 274866631..7c0147f36 100644 --- a/src/methods/equals.jl +++ b/src/methods/equals.jl @@ -6,7 +6,7 @@ export equals ## What is equals? The equals function checks if two geometries are equal. They are equal if they -share the same set of points and edges. +share the same set of points and edges to define the same shape. To provide an example, consider these two lines: ```@example cshape @@ -40,7 +40,8 @@ Note that while we need the same set of points and edges, they don't need to be provided in the same order for polygons. For for example, we need the same set points for two multipoints to be equal, but they don't have to be saved in the same order. This requires checking every point against every other point in the -two geometries we are comparing. +two geometries we are comparing. Additionally, geometries and multi-geometries +can be equal if the multi-geometry only includes that single geometry. =# """ @@ -95,6 +96,14 @@ function equals(::GI.PointTrait, p1, ::GI.PointTrait, p2) return true end +function equals(::GI.PointTrait, p1, ::GI.MultiPointTrait, mp2) + GI.npoint(mp2) == 1 || return false + return equals(p1, GI.getpoint(mp2, 1)) +end + +equals(trait1::GI.MultiPointTrait, mp1, trait2::GI.PointTrait, p2) = + equals(trait2, p2, trait1, mp1) + """ equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)::Bool @@ -121,7 +130,10 @@ end Two curves are equal if they share the same set of points going around the curve. """ -function equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} +function equals( + ::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, l1, + ::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, l2, +) # Check line lengths match n1 = GI.npoint(l1) n2 = GI.npoint(l2) @@ -174,6 +186,14 @@ function equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b) return true end +function equals(::GI.PolygonTrait, geom_a, ::MultiPolygonTrait, geom_b) + GI.npolygon(geom_b) == 1 || return false + return equals(geom_a, GI.getpolygon(geom_b, 1)) +end + +equals(trait_a::GI.MultiPolygonTrait, geom_a, trait_b::PolygonTrait, geom_b) = + equals(trait_b, geom_b, trait_a, geom_a) + """ equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool diff --git a/src/try.jl b/src/try.jl new file mode 100644 index 000000000..4aecdcdeb --- /dev/null +++ b/src/try.jl @@ -0,0 +1,7 @@ +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG + +p2 = LG.Point([0.0, 1.0]) +mp3 = LG.MultiPoint([p2]) +GO.equals(p2, mp3) diff --git a/test/methods/equals.jl b/test/methods/equals.jl index a0b60d6cd..d31adfd87 100644 --- a/test/methods/equals.jl +++ b/test/methods/equals.jl @@ -2,40 +2,45 @@ p1 = LG.Point([0.0, 0.0]) p2 = LG.Point([0.0, 1.0]) # Same points - @test GO.equals(p1, p1) - @test GO.equals(p2, p2) + @test GO.equals(p1, p1) == LG.equals(p1, p1) + @test GO.equals(p2, p2) == LG.equals(p2, p2) # Different points - @test !GO.equals(p1, p2) + @test GO.equals(p1, p2) == LG.equals(p1, p2) mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]]) mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) + mp3 = LG.MultiPoint([p2]) # Same points - @test LG.equals(mp1, mp1) - @test LG.equals(mp2, mp2) + @test GO.equals(mp1, mp1) == LG.equals(mp1, mp1) + @test GO.equals(mp2, mp2) == LG.equals(mp2, mp2) # Different points - @test !LG.equals(mp1, mp2) - @test !LG.equals(mp1, p1) + @test GO.equals(mp1, mp2) == LG.equals(mp1, mp2) + @test GO.equals(mp1, p1) == LG.equals(mp1, p1) + # Point and multipoint + @test GO.equals(p2, mp3) == LG.equals(p2, mp3) end @testset "Lines/Rings" begin l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]]) l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]]) # Equal lines - @test LG.equals(l1, l1) - @test LG.equals(l2, l2) + @test GO.equals(l1, l1) == LG.equals(l1, l1) + @test GO.equals(l2, l2) == LG.equals(l2, l2) # Different lines - @test !LG.equals(l1, l2) && !LG.equals(l2, l1) + @test GO.equals(l1, l2) == GO.equals(l2, l1) == LG.equals(l1, l2) r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) r2 = LG.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) # Equal rings - @test GO.equals(r1, r1) - @test GO.equals(r2, r2) + @test GO.equals(r1, r1) == LG.equals(r1, r1) + @test GO.equals(r2, r2) == LG.equals(r2, r2) # Different rings - @test !GO.equals(r1, r2) && !GO.equals(r2, r1) + @test GO.equals(r1, r2) == GO.equals(r2, r1) == LG.equals(r1, r2) # Equal linear ring and line string - @test !GO.equals(r2, l3) # TODO: should these be equal? + @test GO.equals(r2, l3) == LG.equals(r2, l3) + # Equal linear ring and line + @test GO.equals(l1, GI.Line([(0.0, 0.0), (0.0, 10.0)])) end @testset "Polygons/MultiPolygons" begin @@ -61,18 +66,18 @@ end ] ) # Equal polygon - @test GO.equals(p1, p1) - @test GO.equals(p2, p2) + @test GO.equals(p1, p1) == LG.equals(p1, p1) + @test GO.equals(p2, p2) == LG.equals(p2, p2) # Different polygons - @test !GO.equals(p1, p2) + @test GO.equals(p1, p2) == LG.equals(p1, p2) # Equal polygons with holes - @test GO.equals(p3, p3) + @test GO.equals(p3, p3) == LG.equals(p3, p3) # Same exterior, different hole - @test !GO.equals(p3, p4) + @test GO.equals(p3, p4) == LG.equals(p3, p4) # Same exterior and first hole, has an extra hole - @test !GO.equals(p3, p5) + @test GO.equals(p3, p5) == LG.equals(p3, p5) - p3 = GI.Polygon( + p6 = LG.Polygon( [[ [-53.57208251953125, 28.287451910503744], [-53.33038330078125, 28.29228897739706], @@ -81,7 +86,7 @@ end ]] ) # Complex polygon - @test GO.equals(p3, p3) + @test GO.equals(p6, p6) == LG.equals(p6, p6) m1 = LG.MultiPolygon([ [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]], @@ -98,7 +103,10 @@ end [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]] ]) # Equal multipolygon - @test GO.equals(m1, m1) + @test GO.equals(m1, m1) == LG.equals(m1, m1) # Equal multipolygon with different order - @test GO.equals(m1, m2) + @test GO.equals(m1, m2) == LG.equals(m2, m2) + # Equal polygon to multipolygon + m3 = LG.MultiPolygon([p3]) + @test GO.equals(p1, m3) == LG.equals(p1, m3) end \ No newline at end of file From 881c98aeefe5ddd30010fbbc27ffe97dbed847d6 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 26 Dec 2023 20:10:01 -0800 Subject: [PATCH 09/16] Update area for empty geoms --- .gitignore | 1 + src/methods/area.jl | 9 ++++++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index d6b4b1c05..6977e6f71 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /Manifest.toml /docs/build/ +.vscode/ \ No newline at end of file diff --git a/src/methods/area.jl b/src/methods/area.jl index 23446e726..cc470ecdc 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -71,13 +71,14 @@ computed slighly differently for different geometries: signed_area(geom) = signed_area(GI.trait(geom), geom) # Points -area(::GI.PointTrait, point) = zero(typeof(GI.x(point))) +area(::GI.PointTrait, point) = GI.isempty(point) ? + 0 : zero(typeof(GI.x(point))) signed_area(trait::GI.PointTrait, point) = area(trait, point) # Curves -area(::CT, curve) where CT <: GI.AbstractCurveTrait = - zero(typeof(GI.x(GI.getpoint(curve, 1)))) +area(::CT, curve) where CT <: GI.AbstractCurveTrait = GI.isempty(curve) ? + 0 : zero(typeof(GI.x(GI.getpoint(curve, 1)))) signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = area(trait, curve) @@ -109,6 +110,8 @@ repeating the first point at the end of the coordinates, curve is still assumed to be closed. =# function _signed_area(geom) + # If empty, return zero + isempty(geom) && return 0 # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) first_last_equal = equals(GI.getpoint(geom, 1), GI.getpoint(geom, np)) From fbc25c3b32259c9dce539670d3bce08f4bb8b97b Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 26 Dec 2023 21:14:17 -0800 Subject: [PATCH 10/16] Polygon area base case --- src/methods/area.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/methods/area.jl b/src/methods/area.jl index cc470ecdc..b858fdd88 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -87,6 +87,7 @@ signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) function signed_area(::GI.PolygonTrait, poly) + isempty(poly) && return 0 s_area = _signed_area(GI.getexterior(poly)) area = abs(s_area) # Remove hole areas from total From 44eb25a088b04a4a924abe289e0a14fb5ecfbf1e Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 26 Dec 2023 21:43:26 -0800 Subject: [PATCH 11/16] Update GI call --- src/methods/area.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/methods/area.jl b/src/methods/area.jl index b858fdd88..2ba00322a 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -87,7 +87,7 @@ signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) function signed_area(::GI.PolygonTrait, poly) - isempty(poly) && return 0 + GI.isempty(poly) && return 0 s_area = _signed_area(GI.getexterior(poly)) area = abs(s_area) # Remove hole areas from total @@ -112,7 +112,7 @@ to be closed. =# function _signed_area(geom) # If empty, return zero - isempty(geom) && return 0 + GI.isempty(geom) && return 0 # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) first_last_equal = equals(GI.getpoint(geom, 1), GI.getpoint(geom, np)) From 32ae9230a40c479452a6d6b9930d8a9887a18be9 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 26 Dec 2023 21:53:07 -0800 Subject: [PATCH 12/16] testing updates --- src/methods/area.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/methods/area.jl b/src/methods/area.jl index 2ba00322a..75c49a55e 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -112,6 +112,7 @@ to be closed. =# function _signed_area(geom) # If empty, return zero + println("hi") GI.isempty(geom) && return 0 # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) From d20837d698ff9cd2759106445c9365b00d72d63e Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 26 Dec 2023 22:00:41 -0800 Subject: [PATCH 13/16] Update empty check --- src/methods/area.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/methods/area.jl b/src/methods/area.jl index 75c49a55e..f59196b56 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -111,11 +111,9 @@ repeating the first point at the end of the coordinates, curve is still assumed to be closed. =# function _signed_area(geom) - # If empty, return zero - println("hi") - GI.isempty(geom) && return 0 # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) + np == 0 && return 0 first_last_equal = equals(GI.getpoint(geom, 1), GI.getpoint(geom, np)) np -= first_last_equal ? 1 : 0 # Integrate the area under the curve From f176f95d4c945e572f814451c888db0cbc134b96 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 27 Dec 2023 13:49:38 -0800 Subject: [PATCH 14/16] Added area calculations for empty geoms and collections --- src/methods/area.jl | 39 +++++++++++++++++++++++++++++++++++---- test/methods/area.jl | 31 ++++++++++++++++++++++++++++++- 2 files changed, 65 insertions(+), 5 deletions(-) diff --git a/src/methods/area.jl b/src/methods/area.jl index f59196b56..eacc9ba92 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -76,13 +76,40 @@ area(::GI.PointTrait, point) = GI.isempty(point) ? signed_area(trait::GI.PointTrait, point) = area(trait, point) +# MultiPoints +function area(::GI.MultiPointTrait, multipoint) + GI.isempty(multipoint) && return 0 + np = GI.npoint(multipoint) + np == 0 && return 0 + return zero(typeof(GI.x(GI.getpoint(multipoint, np)))) +end + +signed_area(trait::GI.MultiPointTrait, multipoint) = area(trait, multipoint) + # Curves -area(::CT, curve) where CT <: GI.AbstractCurveTrait = GI.isempty(curve) ? - 0 : zero(typeof(GI.x(GI.getpoint(curve, 1)))) +function area(::CT, curve) where CT <: GI.AbstractCurveTrait + GI.isempty(curve) && return 0 + np = GI.npoint(curve) + np == 0 && return 0 + return zero(typeof(GI.x(GI.getpoint(curve, np)))) +end signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = area(trait, curve) +# MultiCurves +function area(::MCT, multicurve) where MCT <: GI.AbstractMultiCurveTrait + GI.isempty(multicurve) && return 0 + ng = GI.ngeom(multicurve) + ng == 0 && return 0 + np = GI.npoint(GI.getgeom(multicurve, ng)) + np == 0 && return 0 + return zero(typeof(GI.x(GI.getpoint(GI.getgeom(multicurve, ng), np)))) +end + +signed_area(trait::MCT, curve) where MCT <: GI.AbstractMultiCurveTrait = + area(trait, curve) + # Polygons area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) @@ -90,6 +117,7 @@ function signed_area(::GI.PolygonTrait, poly) GI.isempty(poly) && return 0 s_area = _signed_area(GI.getexterior(poly)) area = abs(s_area) + area == 0 && return area # Remove hole areas from total for hole in GI.gethole(poly) area -= abs(_signed_area(hole)) @@ -99,9 +127,12 @@ function signed_area(::GI.PolygonTrait, poly) end # MultiPolygons -area(::GI.MultiPolygonTrait, geom) = - sum((area(poly) for poly in GI.getpolygon(geom))) +area(::GI.MultiPolygonTrait, multipoly) = + sum((area(poly) for poly in GI.getpolygon(multipoly)), init = 0) +# GeometryCollections +area(::GI.GeometryCollectionTrait, collection) = + sum((area(geom) for geom in GI.getgeom(collection)), init = 0) #= Helper function: diff --git a/test/methods/area.jl b/test/methods/area.jl index b69e2797e..a406ebe2e 100644 --- a/test/methods/area.jl +++ b/test/methods/area.jl @@ -1,6 +1,14 @@ pt = LG.Point([0.0, 0.0]) +empty_pt = LG.readgeom("POINT EMPTY") +mpt = LG.MultiPoint([[0.0, 0.0], [1.0, 0.0]]) +empty_mpt = LG.readgeom("MULTIPOINT EMPTY") l1 = LG.LineString([[0.0, 0.0], [0.5, 0.5], [1.0, 0.5]]) +empty_l = LG.readgeom("LINESTRING EMPTY") +ml1 = LG.MultiLineString([[[0.0, 0.0], [0.5, 0.5], [1.0, 0.5]], [[0.0, 0.0], [0.0, 0.1]]]) +empty_ml = LG.readgeom("MULTILINESTRING EMPTY") +empty_l = LG.readgeom("LINESTRING EMPTY") r1 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 2.0], [0.0, 0.0]]) +empty_r = LG.readgeom("LINEARRING EMPTY") p1 = LG.Polygon([ [[10.0, 0.0], [30.0, 0.0], [30.0, 20.0], [10.0, 20.0], [10.0, 0.0]], ]) @@ -19,13 +27,23 @@ p4 = LG.Polygon([ [0.0, 5.0], ], ]) +empty_p = LG.readgeom("POLYGON EMPTY") mp1 = LG.MultiPolygon([p2, p4]) - +empty_mp = LG.readgeom("MULTIPOLYGON EMPTY") +c = LG.GeometryCollection([p1, p2, r1, empty_l]) +empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY") # Points, lines, and rings have zero area @test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0 +# @test GO.area(empty_pt) == LG.area(empty_pt) == 0 +@test GO.area(mpt) == GO.signed_area(mpt) == LG.area(mpt) == 0 +@test GO.area(empty_mpt) == LG.area(empty_mpt) == 0 @test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0 +@test GO.area(empty_l) == LG.area(empty_l) == 0 +@test GO.area(ml1) == GO.signed_area(ml1) == LG.area(ml1) == 0 +@test GO.area(empty_ml) == LG.area(empty_ml) == 0 @test GO.area(r1) == GO.signed_area(r1) == LG.area(r1) == 0 +@test GO.area(empty_r) == LG.area(empty_r) == 0 # Polygons have non-zero area # CCW polygons have positive signed area @@ -41,5 +59,16 @@ a2 = LG.area(p2) a4 = LG.area(p4) @test GO.area(p4) == a4 @test GO.signed_area(p4) == -a4 +# Empty polygon +@test GO.area(empty_p) == LG.area(empty_p) == 0 +@test GO.signed_area(empty_p) == 0 + # Multipolygon calculations work @test GO.area(mp1) == a2 + a4 +# Empty multipolygon +@test GO.area(empty_mp) == LG.area(empty_mp) == 0 + +# Geometry collection summed area +@test GO.area(c) == LG.area(c) +# Empty collection +@test GO.area(empty_c) == LG.area(empty_c) == 0 From 77bfda4ddd68c76042566eee3408852f36026e52 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 27 Dec 2023 15:20:48 -0800 Subject: [PATCH 15/16] Remove try file --- src/try.jl | 7 ------- 1 file changed, 7 deletions(-) delete mode 100644 src/try.jl diff --git a/src/try.jl b/src/try.jl deleted file mode 100644 index 4aecdcdeb..000000000 --- a/src/try.jl +++ /dev/null @@ -1,7 +0,0 @@ -import GeometryOps as GO -import GeoInterface as GI -import LibGEOS as LG - -p2 = LG.Point([0.0, 1.0]) -mp3 = LG.MultiPoint([p2]) -GO.equals(p2, mp3) From aadfc989039e539d7dda772a5a3bf779ae4683ac Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Thu, 28 Dec 2023 14:00:05 -0800 Subject: [PATCH 16/16] Update area and dist with type param --- src/methods/area.jl | 98 ++++++++------------ src/methods/distance.jl | 196 +++++++++++++++++++-------------------- test/methods/area.jl | 12 ++- test/methods/distance.jl | 21 +++++ 4 files changed, 164 insertions(+), 163 deletions(-) diff --git a/src/methods/area.jl b/src/methods/area.jl index eacc9ba92..5b256dbbd 100644 --- a/src/methods/area.jl +++ b/src/methods/area.jl @@ -44,19 +44,25 @@ for polygons. =# """ - area(geom)::Real + area(geom, ::Type{T} = Float64)::T Returns the area of the geometry. This is computed slighly differently for different geometries: - - The area of a point is always zero. - - The area of a curve is always zero. + - The area of a point/multipoint is always zero. + - The area of a curve/multicurve is always zero. - The area of a polygon is the absolute value of the signed area. - The area multi-polygon is the sum of the areas of all of the sub-polygons. + - The area of a geometry collection is the sum of the areas of all of the + sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -area(geom) = area(GI.trait(geom), geom) +area(geom, ::Type{T} = Float64) where T <: AbstractFloat = + _area(T, GI.trait(geom), geom) """ - signed_area(geom)::Real + signed_area(geom, ::Type{T} = Float64)::T Returns the signed area of the geometry, based on winding order. This is computed slighly differently for different geometries: @@ -67,72 +73,43 @@ computed slighly differently for different geometries: counterclockwise. - You cannot compute the signed area of a multipolygon as it doesn't have a meaning as each sub-polygon could have a different winding order. -""" -signed_area(geom) = signed_area(GI.trait(geom), geom) - -# Points -area(::GI.PointTrait, point) = GI.isempty(point) ? - 0 : zero(typeof(GI.x(point))) - -signed_area(trait::GI.PointTrait, point) = area(trait, point) - -# MultiPoints -function area(::GI.MultiPointTrait, multipoint) - GI.isempty(multipoint) && return 0 - np = GI.npoint(multipoint) - np == 0 && return 0 - return zero(typeof(GI.x(GI.getpoint(multipoint, np)))) -end -signed_area(trait::GI.MultiPointTrait, multipoint) = area(trait, multipoint) - -# Curves -function area(::CT, curve) where CT <: GI.AbstractCurveTrait - GI.isempty(curve) && return 0 - np = GI.npoint(curve) - np == 0 && return 0 - return zero(typeof(GI.x(GI.getpoint(curve, np)))) -end +Result will be of type T, where T is an optional argument with a default value +of Float64. +""" +signed_area(geom, ::Type{T} = Float64) where T <: AbstractFloat = + _signed_area(T, GI.trait(geom), geom) -signed_area(trait::CT, curve) where CT <: GI.AbstractCurveTrait = - area(trait, curve) - -# MultiCurves -function area(::MCT, multicurve) where MCT <: GI.AbstractMultiCurveTrait - GI.isempty(multicurve) && return 0 - ng = GI.ngeom(multicurve) - ng == 0 && return 0 - np = GI.npoint(GI.getgeom(multicurve, ng)) - np == 0 && return 0 - return zero(typeof(GI.x(GI.getpoint(GI.getgeom(multicurve, ng), np)))) -end +# Points, MultiPoints, Curves, MultiCurves +_area(::Type{T}, ::GI.AbstractGeometryTrait, geom) where T = zero(T) -signed_area(trait::MCT, curve) where MCT <: GI.AbstractMultiCurveTrait = - area(trait, curve) +_signed_area(::Type{T}, ::GI.AbstractGeometryTrait, geom) where T = zero(T) # Polygons -area(trait::GI.PolygonTrait, geom) = abs(signed_area(trait, geom)) +_area(::Type{T}, trait::GI.PolygonTrait, poly) where T = + abs(_signed_area(T, trait, poly)) -function signed_area(::GI.PolygonTrait, poly) - GI.isempty(poly) && return 0 - s_area = _signed_area(GI.getexterior(poly)) +function _signed_area(::Type{T}, ::GI.PolygonTrait, poly) where T + GI.isempty(poly) && return zero(T) + s_area = _signed_area(T, GI.getexterior(poly)) area = abs(s_area) area == 0 && return area # Remove hole areas from total for hole in GI.gethole(poly) - area -= abs(_signed_area(hole)) + area -= abs(_signed_area(T, hole)) end # Winding of exterior ring determines sign return area * sign(s_area) end -# MultiPolygons -area(::GI.MultiPolygonTrait, multipoly) = - sum((area(poly) for poly in GI.getpolygon(multipoly)), init = 0) +# # MultiPolygons and GeometryCollections +_area( + ::Type{T}, + ::Union{GI.MultiPolygonTrait, GI.GeometryCollectionTrait}, + geoms, +) where T = + sum((area(geom, T) for geom in GI.getgeom(geoms)), init = zero(T)) -# GeometryCollections -area(::GI.GeometryCollectionTrait, collection) = - sum((area(geom) for geom in GI.getgeom(collection)), init = 0) #= Helper function: @@ -141,21 +118,20 @@ to find the area under the curve. Even if curve isn't explicitly closed by repeating the first point at the end of the coordinates, curve is still assumed to be closed. =# -function _signed_area(geom) - # Close curve, even if last point isn't explicitly repeated +function _signed_area(::Type{T}, geom) where T + area = zero(T) + # Close curve, even if last point isn't explicitly repeated np = GI.npoint(geom) - np == 0 && return 0 + np == 0 && return area first_last_equal = equals(GI.getpoint(geom, 1), GI.getpoint(geom, np)) np -= first_last_equal ? 1 : 0 # Integrate the area under the curve p1 = GI.getpoint(geom, np) - T = typeof(GI.x(p1)) - area = zero(T) for i in 1:np p2 = GI.getpoint(geom, i) # Accumulate the area into `area` area += GI.x(p1) * GI.y(p2) - GI.y(p1) * GI.x(p2) p1 = p2 end - return area / 2 + return T(area / 2) end \ No newline at end of file diff --git a/src/methods/distance.jl b/src/methods/distance.jl index c14c295f0..2c2782fa5 100644 --- a/src/methods/distance.jl +++ b/src/methods/distance.jl @@ -54,7 +54,7 @@ polygons, so it isn't implemented for curves. =# """ - distance(point, geom)::Real + distance(point, geom, ::Type{T} = Float64)::T Calculates the ditance from the geometry `g1` to the `point`. The distance will always be positive or zero. @@ -62,8 +62,6 @@ will always be positive or zero. The method will differ based on the type of the geometry provided: - The distance from a point to a point is just the Euclidean distance between the points. - - The distance from a point to a multipolygon is the shortest distance from - a the given point to any point within the multipoint object. - The distance from a point to a line is the minimum distance from the point to the closest point on the given line. - The distance from a point to a linestring is the minimum distance from the @@ -73,18 +71,17 @@ The method will differ based on the type of the geometry provided: - The distance from a point to a polygon is zero if the point is within the polygon and otherwise is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - - The distance from a point to a multipolygon is zero if the point is within - the multipolygon and otherwise is the minimum distance from the point to the - closest edge of any of the polygons within the multipolygon. This includes - edges created by holes of the polygons as well. + - The distance from a point to a multigeometry or a geometry collection is + the minimum distance between the point and any of the sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -distance(point, geom) = distance( - GI.trait(point), point, - GI.trait(geom), geom, -) +distance(point, geom, ::Type{T} = Float64) where T <: AbstractFloat = + _distance(T, GI.trait(point), point, GI.trait(geom), geom) """ - signed_distance(point, geom)::Real + signed_distance(point, geom, ::Type{T} = Float64)::T Calculates the signed distance from the geometry `geom` to the given point. Points within `geom` have a negative signed distance, and points outside of @@ -95,116 +92,114 @@ Points within `geom` have a negative signed distance, and points outside of within the polygon and is positive otherwise. The value of the distance is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - - The signed distance from a point to a mulitpolygon is negative if the - point is within one of the polygons that make up the multipolygon and is - positive otherwise. The value of the distance is the minimum distance from - the point to an edge of the multipolygon. This includes edges created by - holes of the polygons as well. + - The signed distance from a point to a multigeometry or a geometry + collection is the minimum signed distance between the point and any of the + sub-geometries. + +Result will be of type T, where T is an optional argument with a default value +of Float64. """ -signed_distance(point, geom) = signed_distance( - GI.trait(point), point, - GI.trait(geom), geom, -) +signed_distance(point, geom, ::Type{T} = Float64) where T<:AbstractFloat = + _signed_distance(T, GI.trait(point), point, GI.trait(geom), geom) -# # Distance # Swap argument order to point as first argument -distance(gtrait::GI.AbstractTrait, geom, ptrait::GI.PointTrait, point) = - distance(ptrait, point, gtrait, geom) +_distance( + ::Type{T}, + gtrait::GI.AbstractTrait, geom, + ptrait::GI.PointTrait, point, +) where T = _distance(T, ptrait, point, gtrait, geom) + +_signed_distance( + ::Type{T}, + gtrait::GI.AbstractTrait, geom, + ptrait::GI.PointTrait, point, +) where T = _signed_distance(T, ptrait, point, gtrait, geom) -# Point-Point -distance(::GI.PointTrait, point, ::GI.PointTrait, geom) = - euclid_distance(point, geom) +# Point-Point, Point-Line, Point-LineString, Point-LinearRing +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.PointTrait, geom) where T = + _euclid_distance(T, point, geom) -# Point-MultiPoint -function distance(::GI.PointTrait, point, ::GI.MultiPointTrait, geom) - T = typeof(GI.x(point)) - min_dist = typemax(T) - for p in GI.getpoint(geom) - dist = euclid_distance(point, p) - min_dist = dist < min_dist ? dist : min_dist - end - return min_dist -end +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LineTrait, geom) where T = + _distance_line(T, point, GI.getpoint(geom, 1), GI.getpoint(geom, 2)) -# Point-Line -distance(::GI.PointTrait, point, ::GI.LineTrait, geom) = - _distance_line(point, GI.getpoint(geom, 1), GI.getpoint(geom, 2)) +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LineStringTrait, geom) where T = + _distance_curve(T, point, geom, close_curve = false) -# Point-LineString -distance(::GI.PointTrait, point, ::GI.LineStringTrait, geom) = - _distance_curve(point, geom, close_curve = false) +_distance(::Type{T}, ::GI.PointTrait, point, ::GI.LinearRingTrait, geom) where T = + _distance_curve(T, point, geom, close_curve = true) -# Point-LinearRing -distance(::GI.PointTrait, point, ::GI.LinearRingTrait, geom) = - _distance_curve(point, geom, close_curve = true) +_signed_distance(::Type{T}, ptrait::GI.PointTrait, point, gtrait::GI.AbstractTrait, geom) where T = + _distance(T, ptrait, point, gtrait, geom) # Point-Polygon -function distance(::GI.PointTrait, point, ::GI.PolygonTrait, geom) - T = typeof(GI.x(point)) +function _distance(::Type{T}, ::GI.PointTrait, point, ::GI.PolygonTrait, geom) where T GI.within(point, geom) && return zero(T) - return _distance_polygon(point, geom) + return _distance_polygon(T, point, geom) end -# Point-MultiPolygon -function distance(::GI.PointTrait, point, ::GI.MultiPolygonTrait, geom) - min_dist = distance(point, GI.getpolygon(geom, 1)) - for i in 2:GI.npolygon(geom) - min_dist == 0 && return min_dist # point inside of last polygon checked - dist = distance(point, GI.getpolygon(geom, i)) - min_dist = dist < min_dist ? dist : min_dist - end - return min_dist +function _signed_distance(::Type{T}, ::GI.PointTrait, point, ::GI.PolygonTrait, geom) where T + min_dist = _distance_polygon(T, point, geom) + # negative if point is inside polygon + return GI.within(point, geom) ? -min_dist : min_dist end -# # Signed Distance - -# Swap argument order to point as first argument -signed_distance(gtrait::GI.AbstractTrait, geom, ptrait::GI.PointTrait, point) = - signed_distance(ptrait, point, gtrait, geom) -# Point-Point, Point-Line, Point-LineString, Point-LinearRing -signed_distance(ptrait::GI.PointTrait, point, gtrait::GI.AbstractTrait, geom) = - distance(ptrait, point, gtrait, geom) - -# Point-Polygon -function signed_distance(::GI.PointTrait, point, ::GI.PolygonTrait, geom) - min_dist = _distance_polygon(point, geom) - # should be negative if point is inside polygon - return GI.within(point, geom) ? -min_dist : min_dist +# Point-MultiGeometries / Point-GeometryCollections +function _distance( + ::Type{T}, + ::GI.PointTrait, + point, + ::Union{ + GI.MultiPointTrait, GI.MultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, + geoms, +) where T + min_dist = typemax(T) + for g in GI.getgeom(geoms) + dist = distance(point, g, T) + min_dist = dist < min_dist ? dist : min_dist + end + return min_dist end -# Point-Multipolygon -function signed_distance(::GI.PointTrait, point, ::GI.MultiPolygonTrait, geom) - min_dist = signed_distance(point, GI.getpolygon(geom, 1)) - for i in 2:GI.npolygon(geom) - dist = signed_distance(point, GI.getpolygon(geom, i)) +function _signed_distance( + ::Type{T}, + ::GI.PointTrait, + point, + ::Union{ + GI.MultiPointTrait, GI.MultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, + geoms, +) where T + min_dist = typemax(T) + for g in GI.getgeom(geoms) + dist = signed_distance(point, g, T) min_dist = dist < min_dist ? dist : min_dist end return min_dist end - -""" - euclid_distance(p1::Point, p2::Point)::Real - -Returns the Euclidean distance between two points. -""" -Base.@propagate_inbounds euclid_distance(p1, p2) = _euclid_distance( - GeoInterface.x(p1), GeoInterface.y(p1), - GeoInterface.x(p2), GeoInterface.y(p2), -) +# Returns the Euclidean distance between two points. +Base.@propagate_inbounds _euclid_distance(::Type{T}, p1, p2) where T = + _euclid_distance( + T, + GeoInterface.x(p1), GeoInterface.y(p1), + GeoInterface.x(p2), GeoInterface.y(p2), + ) # Returns the Euclidean distance between two points given their x and y values. -Base.@propagate_inbounds _euclid_distance(x1, y1, x2, y2) = - sqrt((x2 - x1)^2 + (y2 - y1)^2) +Base.@propagate_inbounds _euclid_distance(::Type{T}, x1, y1, x2, y2) where T = + T(sqrt((x2 - x1)^2 + (y2 - y1)^2)) #= Returns the minimum distance from point p0 to the line defined by endpoints p1 and p2. =# -function _distance_line(p0, p1, p2) +function _distance_line(::Type{T}, p0, p1, p2) where T x0, y0 = GeoInterface.x(p0), GeoInterface.y(p0) x1, y1 = GeoInterface.x(p1), GeoInterface.y(p1) x2, y2 = GeoInterface.x(p2), GeoInterface.y(p2) @@ -221,16 +216,16 @@ function _distance_line(p0, p1, p2) c1 = sum(w .* v) if c1 <= 0 # p0 is closest to first endpoint - return _euclid_distance(x0, y0, xfirst, yfirst) + return _euclid_distance(T, x0, y0, xfirst, yfirst) end c2 = sum(v .* v) if c2 <= c1 # p0 is closest to last endpoint - return _euclid_distance(x0, y0, xlast, ylast) + return _euclid_distance(T, x0, y0, xlast, ylast) end b2 = c1 / c2 # projection fraction - return _euclid_distance(x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) + return _euclid_distance(T, x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2])) end @@ -239,19 +234,18 @@ Returns the minimum distance from the given point to the given curve. If close_curve is true, make sure to include the edge from the first to last point of the curve, even if it isn't explicitly repeated. =# -function _distance_curve(point, curve; close_curve = false) - # See if linear ring has explicitly repeated last point in coordinates +function _distance_curve(::Type{T}, point, curve; close_curve = false) where T + # see if linear ring has explicitly repeated last point in coordinates np = GI.npoint(curve) first_last_equal = equals(GI.getpoint(curve, 1), GI.getpoint(curve, np)) close_curve &= first_last_equal np -= first_last_equal ? 1 : 0 - # Find minimum distance - T = typeof(GI.x(point)) + # find minimum distance min_dist = typemax(T) p1 = GI.getpoint(curve, close_curve ? np : 1) for i in (close_curve ? 1 : 2):np p2 = GI.getpoint(curve, i) - dist = _distance_line(point, p1, p2) + dist = _distance_line(T, point, p1, p2) min_dist = dist < min_dist ? dist : min_dist p1 = p2 end @@ -263,10 +257,10 @@ Returns the minimum distance from the given point to an edge of the given polygon, including from edges created by holes. Assumes polygon isn't filled and treats the exterior and each hole as a linear ring. =# -function _distance_polygon(point, poly) - min_dist = _distance_curve(point, GI.getexterior(poly); close_curve = true) +function _distance_polygon(::Type{T}, point, poly) where T + min_dist = _distance_curve(T, point, GI.getexterior(poly); close_curve = true) @inbounds for hole in GI.gethole(poly) - dist = _distance_curve(point, hole; close_curve = true) + dist = _distance_curve(T, point, hole; close_curve = true) min_dist = dist < min_dist ? dist : min_dist end return min_dist diff --git a/test/methods/area.jl b/test/methods/area.jl index a406ebe2e..99227c7cd 100644 --- a/test/methods/area.jl +++ b/test/methods/area.jl @@ -35,7 +35,11 @@ empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY") # Points, lines, and rings have zero area @test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0 -# @test GO.area(empty_pt) == LG.area(empty_pt) == 0 +@test GO.area(empty_pt) == LG.area(empty_pt) == 0 +@test GO.area(pt) isa Float64 +@test GO.signed_area(pt, Float32) isa Float32 +@test GO.signed_area(pt) isa Float64 +@test GO.area(pt, Float32) isa Float32 @test GO.area(mpt) == GO.signed_area(mpt) == LG.area(mpt) == 0 @test GO.area(empty_mpt) == LG.area(empty_mpt) == 0 @test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0 @@ -49,6 +53,9 @@ empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY") # CCW polygons have positive signed area @test GO.area(p1) == GO.signed_area(p1) == LG.area(p1) @test GO.signed_area(p1) > 0 +# Float32 calculations +@test GO.area(p1) isa Float64 +@test GO.area(p1, Float32) isa Float32 # CW polygons have negative signed area a2 = LG.area(p2) @test GO.area(p2) == a2 @@ -65,10 +72,13 @@ a4 = LG.area(p4) # Multipolygon calculations work @test GO.area(mp1) == a2 + a4 +@test GO.area(mp1, Float32) isa Float32 # Empty multipolygon @test GO.area(empty_mp) == LG.area(empty_mp) == 0 + # Geometry collection summed area @test GO.area(c) == LG.area(c) +@test GO.area(c, Float32) isa Float32 # Empty collection @test GO.area(empty_c) == LG.area(empty_c) == 0 diff --git a/test/methods/distance.jl b/test/methods/distance.jl index e4c6499e1..7d66b13f8 100644 --- a/test/methods/distance.jl +++ b/test/methods/distance.jl @@ -10,6 +10,8 @@ pt9 = LG.Point([3.5, 3.1]) pt10 = LG.Point([10.0, 10.0]) pt11 = LG.Point([2.5, 7.0]) +mpt1 = LG.MultiPoint([pt1, pt2, pt3]) + l1 = LG.LineString([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0]]) r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [0.0, 0.0]]) @@ -23,12 +25,17 @@ p2 = LG.Polygon(r5) mp1 = LG.MultiPolygon([p1, p2]) +c1 = LG.GeometryCollection([pt1, r1, p1]) + # Point and Point # Distance from point to same point @test GO.distance(pt1, pt1) == LG.distance(pt1, pt1) # Distance from point to different point @test GO.distance(pt1, pt2) ≈ GO.distance(pt2, pt1) ≈ LG.distance(pt1, pt2) +# Return types +@test GO.distance(pt1, pt1) isa Float64 +@test GO.distance(pt1, pt1, Float32) isa Float32 # Point and Line @@ -40,6 +47,9 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.distance(pt3, l1) ≈ GO.distance(l1, pt3) ≈ LG.distance(pt3, l1) # Point closer to one segment than another @test GO.distance(pt4, l1) ≈ GO.distance(l1, pt4) ≈ LG.distance(pt4, l1) +# Return types +@test GO.distance(pt1, l1) isa Float64 +@test GO.distance(pt1, l1, Float32) isa Float32 # Point and Ring @@ -75,6 +85,14 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.distance(pt8, p1) ≈ LG.distance(pt8, p1) @test GO.signed_distance(pt8, p1) ≈ LG.distance(pt8, p1) @test GO.distance(pt9, p1) ≈ LG.distance(pt9, p1) +# Return types +@test GO.distance(pt1, p1) isa Float64 +@test GO.distance(pt1, p1, Float32) isa Float32 + +# Point and MultiPoint +@test GO.distance(pt4, mpt1) == LG.distance(pt4, mpt1) +@test GO.distance(pt4, mpt1) isa Float64 +@test GO.distance(pt4, mpt1, Float32) isa Float32 # Point and MultiPolygon # Point outside of either polygon @@ -88,3 +106,6 @@ mp1 = LG.MultiPolygon([p1, p2]) @test GO.signed_distance(pt11, mp1) ≈ -(min(LG.distance(pt11, r2), LG.distance(pt11, r3), LG.distance(pt11, r4), LG.distance(pt11, r5))) +# Point and Geometry Collection +@test GO.distance(pt1, c1) == LG.distance(pt1, c1) +