diff --git a/Project.toml b/Project.toml index 0fe3b53d8..318c474d8 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" @@ -27,4 +26,7 @@ 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/bools.jl b/src/methods/bools.jl index fd6cffa6f..aac4f8075 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -265,66 +265,63 @@ 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 vs opne rings - if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end) - n -= 1 + + # 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 end + # Loop over all points in the ring - for i in 1:(n - 1) - # First point on edge + for i in 1:l - 1 + j = i + 1 + p_i = GI.getpoint(ring, i) - 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 - ) + 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) + on_boundary && return !ignore_boundary - # Check if ray from point passes through edge - intersects = ( - (yi > y) !== (yj > y) && - (x < (xj - xi) * (y - yi) / (yj - yi) + xi) - ) + + intersects = ((yi > GI.y(pt)) !== (yj > GI.y(pt))) && + (GI.x(pt) < (xj - xi) * (GI.y(pt) - yi) / (yj - yi) + xi) + if intersects inside = !inside end end + return inside end @@ -344,7 +341,6 @@ 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 @@ -369,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 the line of poly1 does not intersect the line of poly2 - #intersects(poly1, 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 - # poly1 must be in poly2 - return true + # Check the line of poly1 does not intersect the line of poly2 + line_intersects(poly1, poly2) && return false + + # poly1 must be in poly2 + return true end diff --git a/src/methods/crosses.jl b/src/methods/crosses.jl index 3aa62d62e..7c215c857 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 line_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) - intersects(line, l) && return true + line_intersects(line, l) && return true end return false end diff --git a/src/methods/disjoint.jl b/src/methods/disjoint.jl index 02b7d4e46..b51e5ab66 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 !intersects(poly1, poly2) + return !line_intersects(poly1, poly2) end diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl index 78f5784f1..e0476bd81 100644 --- a/src/methods/intersects.jl +++ b/src/methods/intersects.jl @@ -1,67 +1,21 @@ # # Intersection checks -export intersects, intersection, intersection_points +export intersects, intersection -#= -## What is `intersects` vs `intersection` vs `intersection_points`? +# This code checks whether geometries intersect with each other. -The `intersects` methods check whether two geometries intersect with each other. -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. +# !!! note +# This does not compute intersections, only checks if they exist. -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. Finally, the `intersection_points` method returns a list -of tuple points. - -To provide an example, consider these two lines: -```@example intersects_intersection -using GeometryOps -using GeometryOps.GeometryBasics -using Makie -using CairoMakie -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]) -``` -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 -``` - -## Implementation - -This is the GeoInterface-compatible implementation. - -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_CLOSED = 0 const MEETS_OPEN = 1 +const MEETS_CLOSED = 0 """ - intersects(geom1, geom2; kw...)::Bool + line_intersects(line_a, line_b) -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. +Check if `line_a` intersects with `line_b`. + +These can be `LineTrait`, `LineStringTrait` or `LinearRingTrait` ## Example @@ -70,85 +24,41 @@ 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.intersects(line1, line2) +GO.line_intersects(line1, line2) # output true ``` """ -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) +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) a1 = _tuple_point(GI.getpoint(a, 1)) - a2 = _tuple_point(GI.getpoint(a, 2)) b1 = _tuple_point(GI.getpoint(b, 1)) + a2 = _tuple_point(GI.getpoint(a, 2)) b2 = _tuple_point(GI.getpoint(b, 2)) - meet_type = ExactPredicates.meet(a1, a2, b1, b2) - return meet_type == MEETS_OPEN || meet_type == meets + return ExactPredicates.meet(a1, a2, b1, b2) == meets end - -""" - 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( - trait_a::GI.AbstractTrait, a, - trait_b::GI.AbstractTrait, b; - kw..., -) +function line_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...) || - within(trait_a, a, trait_b, b) || within(trait_b, b, trait_a, a) + return line_intersects(edges_a, edges_b; kw...) end - -""" - _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, -) +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 - meet_type = ExactPredicates.meet(edge_a..., edge_b...) - (meet_type == MEETS_OPEN || meet_type == meets) && return true + ExactPredicates.meet(edge_a..., edge_b...) == meets && return true end end return false end """ - intersection(geom_a, geom_b)::Union{Tuple{::Real, ::Real}, ::Nothing} + line_intersection(line_a, line_b) + +Find a point that intersects LineStrings with two coordinates each. -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 +Returns `nothing` if no point is found. ## Example @@ -157,193 +67,59 @@ 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.intersection(line1, line2) +GO.line_intersection(line1, line2) # output (125.58375366067547, -14.83572303404496) ``` """ -intersection(geom_a, geom_b) = - intersection(GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) - -""" - intersection( - ::GI.LineTrait, line_a, - ::GI.LineTrait, line_b, - )::Union{ - ::Tuple{::Real, ::Real}, - ::Nothing - } - -Calculates the intersection between two line segments. Return nothing if -there isn't one. -""" -function intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) - # Get start and end points for both lines +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) + 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) + isnothing(x) || push!(result, x) + end + end + return result +end +function line_intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b) a1 = GI.getpoint(line_a, 1) - a2 = GI.getpoint(line_a, 2) b1 = GI.getpoint(line_b, 1) + a2 = GI.getpoint(line_a, 2) b2 = GI.getpoint(line_b, 2) - # Determine the intersection 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 + return _line_intersection((a1, a2), (b1, b2)) 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 = 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 - 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 +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 result + return nothing end - return nothing -end - -""" - _intersection_point( - (a1, a2)::Tuple, - (b1, b2)::Tuple, - ) -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: - (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 - Both are ::Nothing if point doesn't exist! + ã = a / d + b̃ = b / d -Calculation derivation can be found here: - https://stackoverflow.com/questions/563198/ -""" -function _intersection_point((a1, a2)::Tuple, (b1, b2)::Tuple) - # First line runs from p to p + r - 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(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 - Δ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) + if ã >= 0 && ã <= 1 && b̃ >= 0 && b̃ <= 1 + x = x1 + (ã * (x2 - x1)) + y = y1 + (ã * (y2 - y1)) + return (x, y) end - return nothing, nothing + + return nothing end diff --git a/src/methods/overlaps.jl b/src/methods/overlaps.jl index 6d84f393b..b846e43de 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 intersects(g1, g2) + return line_intersects(g1, g2) end function overlaps(t1::MultiPolygonTrait, mp, t2::PolygonTrait, p1)::Bool for p2 in GI.getgeom(mp) - overlaps(p1, p2) && return true + overlaps(p1, thp2) && return true end end function overlaps(::MultiPolygonTrait, g1, ::MultiPolygonTrait, g2)::Bool diff --git a/src/methods/within.jl b/src/methods/within.jl index 16366f944..c930ce62f 100644 --- a/src/methods/within.jl +++ b/src/methods/within.jl @@ -23,21 +23,11 @@ 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, GI.geometry(g2)) -# Points in geometries +within(::Any, g1, t2::GI.FeatureTrait, g2)::Bool = within(g1, geometry(g2)) 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) -# Polygons within geometries +within(::GI.LineStringTrait, g1, ::GI.LineStringTrait, g2)::Bool = line_on_line(g1, g2) 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/src/transformations/extent.jl b/src/transformations/extent.jl index a5e180d76..2e230672d 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) -extent_applicator(::Union{AbstractCurveTrait,MultiPointTrait}, point) = point +function 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 b95b78172..dc6c078eb 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},Tuple{Float64,Float64}}` edges. +vector of `Tuple{Tuple{Float64,Float64},{Float64,Float64}}` edges. """ function to_edges(x) edges = Vector{Edge}(undef, _nedge(x)) diff --git a/test/methods/bools.jl b/test/methods/bools.jl index 791f0598e..b7650cb87 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.intersection(line8, line9)[1] .≈ (125.583754, -14.835723)) + @test all(GO.line_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.intersection(line10, line11) + points = GO.line_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) == true # this was false before... why? + @test GO.overlaps(pl3, pl4) == false mp1 = GI.MultiPoint([ (-36.05712890625, 26.480407161007275), diff --git a/test/methods/intersects.jl b/test/methods/intersects.jl deleted file mode 100644 index f3d35c68f..000000000 --- a/test/methods/intersects.jl +++ /dev/null @@ -1,162 +0,0 @@ -@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 - 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 - 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 - 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 - - # Multi-polygon that intersects with linestring - -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7c96de785..c4fc39f3e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,7 +18,6 @@ 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