diff --git a/.gitignore b/.gitignore index 6977e6f71..0e58cf59b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ /Manifest.toml /docs/build/ -.vscode/ \ No newline at end of file +.vscode/ diff --git a/Project.toml b/Project.toml index f6787b264..35b23ebe6 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,6 @@ Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] -GeometryBasics = "0.4.7" Proj = "1" julia = "1.3" diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c86ac5fef..95a33ed8e 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -23,16 +23,20 @@ include("methods/bools.jl") include("methods/distance.jl") include("methods/area.jl") include("methods/centroid.jl") -include("methods/intersects.jl") -include("methods/contains.jl") -include("methods/crosses.jl") -include("methods/disjoint.jl") -include("methods/overlaps.jl") -include("methods/within.jl") +include("methods/geom_relations/contains.jl") +include("methods/geom_relations/coveredby.jl") +include("methods/geom_relations/covers.jl") +include("methods/geom_relations/crosses.jl") +include("methods/geom_relations/disjoint.jl") +include("methods/geom_relations/geom_geom_processors.jl") +include("methods/geom_relations/intersects.jl") +include("methods/geom_relations/overlaps.jl") +include("methods/geom_relations/within.jl") include("methods/polygonize.jl") include("methods/barycentric.jl") include("methods/equals.jl") - +include("methods/orientation.jl") +include("methods/geom_relations/touches.jl") include("transformations/extent.jl") include("transformations/flip.jl") include("transformations/simplify.jl") diff --git a/src/methods/bools.jl b/src/methods/bools.jl index 30b8716e1..1208bd555 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -1,8 +1,6 @@ # # Boolean conditions export isclockwise, isconcave -export point_on_line, point_in_polygon, point_in_ring -export line_on_line, line_in_polygon, polygon_in_polygon # These are all adapted from Turf.jl. @@ -89,297 +87,3 @@ function isconcave(poly)::Bool return false end - -# """ -# isparallel(line1::LineString, line2::LineString)::Bool - -# Return `true` if each segment of `line1` is parallel to the correspondent segment of `line2` - -# ## Examples -# ```julia -# import GeoInterface as GI, GeometryOps as GO -# julia> line1 = GI.LineString([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)]) -# GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)], nothing, nothing) - -# julia> line2 = GI.LineString([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)]) -# GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)], nothing, nothing) - -# julia> -# GO.isparallel(line1, line2) -# true -# ``` -# """ -# function isparallel(line1, line2)::Bool -# seg1 = linesegment(line1) -# seg2 = linesegment(line2) - -# for i in eachindex(seg1) -# coors2 = nothing -# coors1 = seg1[i] -# coors2 = seg2[i] -# _isparallel(coors1, coors2) == false && return false -# end -# return true -# end - -# @inline function _isparallel(p1, p2) -# slope1 = bearing_to_azimuth(rhumb_bearing(GI.x(p1), GI.x(p2))) -# slope2 = bearing_to_azimuth(rhumb_bearing(GI.y(p1), GI.y(p2))) - -# return slope1 === slope2 -# end - - -""" - point_on_line(point::Point, line::LineString; ignore_end_vertices::Bool=false)::Bool - -Return true if a point is on a line. Accept a optional parameter to ignore the -start and end vertices of the linestring. - -## Examples - -```jldoctest -import GeoInterface as GI, GeometryOps as GO - -point = (1, 1) -line = GI.LineString([(0, 0), (3, 3), (4, 4)]) -GO.point_on_line(point, line) - -# output -true -``` -""" -function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool - line_points = tuple_points(line) - n = length(line_points) - - exclude_boundary = :none - for i in 1:n - 1 - if ignore_end_vertices - if i === 1 - exclude_boundary = :start - elseif i === n - 2 - exclude_boundary = :end - elseif (i === 1 && i + 1 === n - 1) - exclude_boundary = :both - end - end - if point_on_segment(point, (line_points[i], line_points[i + 1]); exclude_boundary) - return true - end - end - 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) - x2, y2 = GI.x(stop), GI.y(stop) - - dxc = x - x1 - dyc = y - y1 - dx1 = x2 - x1 - dy1 = y2 - y1 - - # TODO use better predicate for crossing here - cross = dxc * dy1 - dyc * dx1 - cross != 0 && return false - - # Will constprop optimise these away? - if exclude_boundary === :none - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 <= x && x <= x2 : x2 <= x && x <= x1 - end - return dy1 > 0 ? y1 <= y && y <= y2 : y2 <= y && y <= y1 - elseif exclude_boundary === :start - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 < x && x <= x2 : x2 <= x && x < x1 - end - return dy1 > 0 ? y1 < y && y <= y2 : y2 <= y && y < y1 - elseif exclude_boundary === :end - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 <= x && x < x2 : x2 < x && x <= x1 - end - return dy1 > 0 ? y1 <= y && y < y2 : y2 < y && y <= y1 - elseif exclude_boundary === :both - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 < x && x < x2 : x2 < x && x < x1 - end - return dy1 > 0 ? y1 < y && y < y2 : y2 < y && y < y1 - end - return false -end - -""" - point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignore_boundary::Bool=false)::Bool - -Take a Point and a Polygon and determine if the point -resides inside the polygon. The polygon can be convex or concave. The function accounts for holes. - -## Examples - -```jldoctest -import GeoInterface as GI, GeometryOps as GO - -point = (-77.0, 44.0) -poly = GI.Polygon([[(-81, 41), (-81, 47), (-72, 47), (-72, 41), (-81, 41)]]) -GO.point_in_polygon(point, poly) - -# output -true -``` -""" -point_in_polygon(point, polygon; kw...)::Bool = - point_in_polygon(GI.trait(point), point, GI.trait(polygon), polygon; kw...) -function point_in_polygon( - ::PointTrait, point, - ::PolygonTrait, poly; - ignore_boundary::Bool=false, - check_extent::Bool=false, -)::Bool - # Cheaply check that the point is inside the polygon extent - if check_extent - point_in_extent(point, GI.extent(poly)) || return false - end - - # Then check the point is inside the exterior ring - 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 - 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 - end - # Loop over all points in the ring - for i in 1:(n - 1) - # First point on edge - 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 - ) - 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) - ) - if intersects - inside = !inside - end - end - return inside -end - -function point_in_extent(p, extent::Extents.Extent) - (x1, x2), (y1, y1) = extent.X, extent.Y - return x1 <= GI.x(p) && y1 <= GI.y(p) && x2 >= GI.x(p) && y2 >= GI.y(p) -end - -line_on_line(line1, line2) = line_on_line(trait(line1), line1, trait(line2), line2) -function line_on_line(t1::GI.AbstractCurveTrait, line1, t2::AbstractCurveTrait, line2) - for p in GI.getpoint(line1) - # FIXME: all points being on the line doesn't - # actually mean the whole line is on the line... - point_on_line(p, line2) || return false - end - return true -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 -) - Extents.intersects(GI.extent(poly), GI.extent(line)) || return false - - inside = false - for i in 1:GI.npoint(line) - 1 - p = GI.getpoint(line, i) - p2 = GI.getpoint(line, i + 1) - point_in_polygon(p, poly) || return false - if !inside - inside = point_in_polygon(p, poly; ignore_boundary=true) - end - # FIXME This seems like a hack, we should check for intersections rather than midpoint?? - if !inside - mid = ((GI.x(p) + GI.x(p2)) / 2, (GI.y(p) + GI.y(p2)) / 2) - inside = point_in_polygon(mid, poly; ignore_boundary=true) - end - end - return inside -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 - - # Check the line of poly1 does not intersect the line of poly2 - #intersects(poly1, poly2) && return false - - # poly1 must be in poly2 - return true - end diff --git a/src/methods/contains.jl b/src/methods/contains.jl deleted file mode 100644 index a58e2a48b..000000000 --- a/src/methods/contains.jl +++ /dev/null @@ -1,25 +0,0 @@ -# # Containment - -export contains - -""" - contains(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool - -Return true if the second geometry is completely contained by the first geometry. -The interiors of both geometries must intersect and, the interior and boundary of the secondary (geometry b) -must not intersect the exterior of the primary (geometry a). -`contains` returns the exact opposite result of `within`. - -## Examples - -```jldoctest -import GeometryOps as GO, GeoInterface as GI -line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) -point = (1, 2) - -GO.contains(line, point) -# output -true -``` -""" -contains(g1, g2)::Bool = within(g2, g1) diff --git a/src/methods/crosses.jl b/src/methods/crosses.jl deleted file mode 100644 index f8a580db0..000000000 --- a/src/methods/crosses.jl +++ /dev/null @@ -1,92 +0,0 @@ -# # Crossing checks - -""" - crosses(geom1, geom2)::Bool - -Return `true` if the intersection results in a geometry whose dimension is one less than -the maximum dimension of the two source geometries and the intersection set is interior to -both source geometries. - -TODO: broken - -## Examples -```julia -import GeoInterface as GI, GeometryOps as GO - -line1 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) -line2 = GI.LineString([(-2, 2), (4, 2)]) - -GO.crosses(line1, line2) -# output -true -``` -""" -crosses(g1, g2)::Bool = crosses(trait(g1), g1, trait(g2), g2)::Bool -crosses(t1::FeatureTrait, g1, t2, g2)::Bool = crosses(GI.geometry(g1), g2) -crosses(t1, g1, t2::FeatureTrait, g2)::Bool = crosses(g1, geometry(g2)) -crosses(::MultiPointTrait, g1, ::LineStringTrait, g2)::Bool = multipoint_crosses_line(g1, g2) -crosses(::MultiPointTrait, g1, ::PolygonTrait, g2)::Bool = multipoint_crosses_poly(g1, g2) -crosses(::LineStringTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_crosses_lines(g2, g1) -crosses(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = line_crosses_poly(g1, g2) -crosses(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_line(g1, g2) -crosses(::PolygonTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_crosses_poly(g2, g1) -crosses(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_poly(g2, g1) - -function multipoint_crosses_line(geom1, geom2) - int_point = false - ext_point = false - i = 1 - np2 = GI.npoint(geom2) - - while i < GI.npoint(geom1) && !int_point && !ext_point - for j in 1:GI.npoint(geom2) - 1 - exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both - if point_on_segment(GI.getpoint(geom1, i), (GI.getpoint(geom2, j), GI.getpoint(geom2, j + 1)); exclude_boundary) - int_point = true - else - ext_point = true - end - end - i += 1 - end - - return int_point && ext_point -end - -function line_crosses_line(line1, line2) - np2 = GI.npoint(line2) - 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 - pa = GI.getpoint(line1, i) - pb = GI.getpoint(line1, i + 1) - p = GI.getpoint(line2, j) - point_on_segment(p, (pa, pb); exclude_boundary) && return true - end - end - end - return false -end - -function line_crosses_poly(line, poly) - for l in flatten(AbstractCurveTrait, poly) - intersects(line, l) && return true - end - return false -end - -function multipoint_crosses_poly(mp, poly) - int_point = false - ext_point = false - - for p in GI.getpoint(mp) - if point_in_polygon(p, poly) - int_point = true - else - ext_point = true - end - int_point && ext_point && return true - end - return false -end diff --git a/src/methods/disjoint.jl b/src/methods/disjoint.jl deleted file mode 100644 index 02b7d4e46..000000000 --- a/src/methods/disjoint.jl +++ /dev/null @@ -1,42 +0,0 @@ -# # Disjointness checks - -""" - disjoint(geom1, geom2)::Bool - -Return `true` if the intersection of the two geometries is an empty set. - -# Examples - -```jldoctest -import GeometryOps as GO, GeoInterface as GI - -poly = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]]) -point = (1, 1) -GO.disjoint(poly, point) - -# output -true -``` -""" -disjoint(g1, g2)::Bool = disjoint(trait(g1), g1, trait(g2), g2) -disjoint(::FeatureTrait, g1, ::Any, g2)::Bool = disjoint(GI.geometry(g1), g2) -disjoint(::Any, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2)) -disjoint(::PointTrait, g1, ::PointTrait, g2)::Bool = !point_equals_point(g1, g2) -disjoint(::PointTrait, g1, ::LineStringTrait, g2)::Bool = !point_on_line(g1, g2) -disjoint(::PointTrait, g1, ::PolygonTrait, g2)::Bool = !point_in_polygon(g1, g2) -disjoint(::LineStringTrait, g1, ::PointTrait, g2)::Bool = !point_on_line(g2, g1) -disjoint(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = !line_on_line(g1, g2) -disjoint(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = !line_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::PointTrait, g2)::Bool = !point_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = !line_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool = polygon_disjoint(g2, g1) - -function polygon_disjoint(poly1, poly2) - for point in GI.getpoint(poly1) - point_in_polygon(point, poly2) && return false - end - for point in GI.getpoint(poly2) - point_in_polygon(point, poly1) && return false - end - return !intersects(poly1, poly2) -end diff --git a/src/methods/equals.jl b/src/methods/equals.jl index f4b710b81..c797432d5 100644 --- a/src/methods/equals.jl +++ b/src/methods/equals.jl @@ -9,7 +9,7 @@ The equals function checks if two geometries are equal. They are equal if they share the same set of points and edges to define the same shape. To provide an example, consider these two lines: -```@example cshape +```@example equals using GeometryOps using GeometryOps.GeometryBasics using Makie @@ -24,7 +24,7 @@ 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 +```@example equals equals(l1, l2) # returns false ``` diff --git a/src/methods/geom_relations/contains.jl b/src/methods/geom_relations/contains.jl new file mode 100644 index 000000000..2bd6a391a --- /dev/null +++ b/src/methods/geom_relations/contains.jl @@ -0,0 +1,64 @@ +# # Contains + +export contains + +#= +## What is contains? + +The contains function checks if a given geometry completly contains another +geometry, or in other words, that the second geometry is completly within the +first. This requires that the two interiors intersect and that the interior and +boundary of the second geometry is not in the exterior of the first geometry. + +To provide an example, consider these two lines: +```@example contains +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(0.25, 0.0), (0.75, 0.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 all of the points and edges of l2 are within l1, so l1 contains +l2. However, l2 does not contain l1. +```@example contains +contains(l1, l2) # returns true +contains(l2, l1) # returns false +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +Given that contains is the exact opposite of within, we simply pass the two +inputs variables, swapped in order, to within. +=# + +""" + contains(g1::AbstractGeometry, g2::AbstractGeometry)::Bool + +Return true if the second geometry is completely contained by the first +geometry. The interiors of both geometries must intersect and the interior and +boundary of the secondary (g2) must not intersect the exterior of the first +(g1). + +`contains` returns the exact opposite result of `within`. + +## Examples + +```jldoctest +import GeometryOps as GO, GeoInterface as GI +line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +point = GI.Point((1, 2)) + +GO.contains(line, point) +# output +true +``` +""" +contains(g1, g2) = GeometryOps.within(g2, g1) diff --git a/src/methods/geom_relations/coveredby.jl b/src/methods/geom_relations/coveredby.jl new file mode 100644 index 000000000..45223f122 --- /dev/null +++ b/src/methods/geom_relations/coveredby.jl @@ -0,0 +1,266 @@ +# # CoveredBy + +export coveredby + +#= +## What is coveredby? + +The coveredby function checks if one geometry is covered by another geometry. +This is an extension of within that does not require the interiors of the two +geometries to intersect, but still does require that the interior and boundary +of the first geometry isn't outside of the second geometry. + +To provide an example, consider this point and line: +```@example coveredby +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +p1 = Point(0.0, 0.0) +p2 = Point(1.0, 1.0) +l1 = Line(p1, p2) + +f, a, p = lines([p1, p2]) +scatter!(p1, color = :red) +``` +As we can see, `p1` is on the endpoint of l1. This means it is not `within`, but +it does meet the definition of `coveredby`. +```@example coveredby +coveredby(p1, l1) # true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`coveredby` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 (either through + overlap or crossing for lines) + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are not allowed to be in the exterior of g2 + - no points of g1 are required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# +const COVEREDBY_ALLOWS = (in_allow = true, on_allow = true, out_allow = false) +const COVEREDBY_CURVE_ALLOWS = (over_allow = true, cross_allow = true, on_allow = true, out_allow = false) +const COVEREDBY_CURVE_REQUIRES = (in_require = false, on_require = false, out_require = false) +const COVEREDBY_POLYGON_REQUIRES = (in_require = true, on_require = false, out_require = false,) + +""" + coveredby(g1, g2)::Bool + +Return `true` if the first geometry is completely covered by the second +geometry. The interior and boundary of the primary geometry (g1) must not +intersect the exterior of the secondary geometry (g2). + +Furthermore, `coveredby` returns the exact opposite result of `covers`. They are +equivalent with the order of the arguments swapped. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI +p1 = GI.Point(0.0, 0.0) +p2 = GI.Point(1.0, 1.0) +l1 = GI.Line([p1, p2]) + +GO.coveredby(p1, l1) +# output +true +``` +""" +coveredby(g1, g2) = _coveredby(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_coveredby(::GI.FeatureTrait, g1, ::Any, g2) = coveredby(GI.geometry(g1), g2) +_coveredby(::Any, g1, t2::GI.FeatureTrait, g2) = coveredby(g1, GI.geometry(g2)) + + +# # Points coveredby geometries + +# Point is coveredby another point if those points are equal +_coveredby( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = equals(g1, g2) + +# Point is coveredby a line/linestring if it is on a line vertex or an edge +_coveredby( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + COVEREDBY_ALLOWS..., + closed_curve = false, +) + +# Point is coveredby a linearring if it is on a vertex or an edge of ring +_coveredby( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + COVEREDBY_ALLOWS..., + closed_curve = true, +) + +# Point is coveredby a polygon if it is inside polygon, including edges/vertices +_coveredby( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., +) + +# Points cannot cover any geometry other than points +_coveredby( + ::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::GI.PointTrait, g2, +) = false + + +# # Lines coveredby geometries + +#= Linestring is coveredby a line if all interior and boundary points of the +first line are on the interior/boundary points of the second line. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is coveredby a ring if all interior and boundary points of the +line are on the edges of the ring. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is coveredby a polygon if all interior and boundary points of the +line are in the polygon interior or on its edges, inlcuding hole edges. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, +) + +# # Rings covered by geometries + +#= Linearring is covered by a line if all vertices and edges of the ring are on +the edges and vertices of the line. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, + closed_curve = false, +) + +#= Linearring is covered by another linear ring if all vertices and edges of the +first ring are on the edges/vertices of the second ring. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is coveredby a polygon if all vertices and edges of the ring are +in the polygon interior or on the polygon edges, inlcuding hole edges. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, +) + + +# # Polygons covered by geometries + +#= Polygon is covered by another polygon if if the interior and edges of the +first polygon are in the second polygon interior or on polygon edges, including +hole edges.=# +_coveredby( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_POLYGON_REQUIRES..., +) + +# Polygons cannot covered by any curves +_coveredby( + ::GI.PolygonTrait, g1, + ::GI.AbstractCurveTrait, g2, +) = false + + +# # Geometries coveredby multi-geometry/geometry collections + +#= Geometry is covered by a multi-geometry or a collection if one of the elements +of the collection cover the geometry. =# +function _coveredby( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + coveredby(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections coveredby geometries + +#= Multi-geometry or a geometry collection is covered by a geometry if all +elements of the collection are covered by the geometry. =# +function _coveredby( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !coveredby(sub_g1, g2) && return false + end + return true +end + diff --git a/src/methods/geom_relations/covers.jl b/src/methods/geom_relations/covers.jl new file mode 100644 index 000000000..714d14643 --- /dev/null +++ b/src/methods/geom_relations/covers.jl @@ -0,0 +1,63 @@ +# # Covers + +export covers + +#= +## What is covers? + +The covers function checks if a given geometry completly covers another +geometry. For this to be true, the "contained" geometry's interior and +boundaries must be covered by the "covering" geometry's interior and boundaries. +The interiors do not need to overlap. + +To provide an example, consider these two lines: +```@example covers +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +p1 = Point(0.0, 0.0) +p2 = Point(1.0, 1.0) +l1 = Line(p1, p2) + +f, a, p = lines([p1, p2]) +scatter!(p1, color = :red) +``` + +```@example covers +covers(l1, p1) # returns true +covers(p1, l1) # returns false +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +Given that covers is the exact opposite of coveredby, we simply pass the two +inputs variables, swapped in order, to coveredby. +=# + +""" + covers(g1::AbstractGeometry, g2::AbstractGeometry)::Bool + +Return true if the first geometry is completely covers the second geometry, +The exterior and boundary of the second geometry must not be outside of the +interior and boundary of the first geometry. However, the interiors need not +intersect. + +`covers` returns the exact opposite result of `coveredby`. + +## Examples + +```jldoctest +import GeometryOps as GO, GeoInterface as GI +l1 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +l2 = GI.LineString([(1, 1), (1, 2)]) + +GO.covers(l1, l2) +# output +true +``` +""" +covers(g1, g2)::Bool = GeometryOps.coveredby(g2, g1) diff --git a/src/methods/geom_relations/crosses.jl b/src/methods/geom_relations/crosses.jl new file mode 100644 index 000000000..40468d6d0 --- /dev/null +++ b/src/methods/geom_relations/crosses.jl @@ -0,0 +1,215 @@ +# # Crosses + +export crosses + +#= +## What is crosses? + +The crosses function checks if one geometry is crosses another geometry. +A geometry can only cross another geometry if they are either two lines, or if +the two geometries have different dimensionalities. If checking two lines, they +must meet in one point. If checking two geometries of different dimensions, the +interiors must meet in at least one point and at least one of the geometries +must have a point outside of the other geometry. + +Note that points can't cross any geometries, despite different dimension, due to +their inability to be both interior and exterior to any other shape. + +To provide an example, consider these two lines: +```@example crosses +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = Line([Point(0.0, 0.0), Point(1.0, 0.0)]) +l2 = Line([Point(0.5, 1.0), Point(0.5, -1.0)]) + +f, a, p = lines(l1) +lines!(l2) +``` +We can see that these two lines cross at their midpoints. +```@example crosses +crosses(l1, l2) # true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`crosses` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 (only through + crossing and NOT overlap for lines) + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - at least one point of g1 are required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - at least one point of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const CROSSES_CURVE_ALLOWS = (over_allow = false, cross_allow = true, on_allow = true, out_allow = true) +const CROSSES_POLYGON_ALLOWS = (in_allow = true, on_allow = true, out_allow = true) +const CROSSES_REQUIRES = (in_require = true, on_require = false, out_require = true) + +""" + crosses(geom1, geom2)::Bool + +Return `true` if the first geometry crosses the second geometry. If they are +both lines, they must meet in one point. Otherwise, they must be of different +dimensions, the interiors must intersect, and the interior of the first geometry +must intersect the exterior of the secondary geometry. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI + +l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)]) +l2 = GI.Line([(0.5, 1.0), (0.5, -1.0)]) + +GO.crosses(l1, l2) +# output +true +``` +""" +crosses(g1, g2) = _crosses(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_crosses(::GI.FeatureTrait, g1, ::Any, g2) = crosses(GI.geometry(g1), g2) +_crosses(::Any, g1, t2::GI.FeatureTrait, g2) = crosses(g1, GI.geometry(g2)) + + +# # Non-specified geometries + +# Points and geometries with the same dimensions D where D ≂̸ 1 default to false +_crosses( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g2, +) = false + + +# # Lines cross geometries + +#= Linestring crosses another linestring if the intersection of the two lines +is exlusivly points (only cross intersections) =# +_crosses( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + CROSSES_CURVE_ALLOWS..., + CROSSES_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring crosses a linearring if the intersection of the line and ring is +exlusivly points (only cross intersections) =# +_crosses( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + CROSSES_CURVE_ALLOWS..., + CROSSES_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring crosses a polygon if at least some of the line interior is in the +polygon interior and some of the line interior is exterior to the polygon. =# +_crosses( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + CROSSES_POLYGON_ALLOWS..., + CROSSES_REQUIRES..., + closed_line = false, +) + + +# # Rings cross geometries + +#= Linearring crosses a linestring if the intersection of the line and ring is +exlusivly points (only cross intersections) =# +_crosses( + trait1::GI.LinearRingTrait, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _crosses(trait2, g2, trait1, g1) + +#= Linearring crosses another ring if the intersection of the two rings is +exlusivly points (only cross intersections) =# +_crosses( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + CROSSES_CURVE_ALLOWS..., + CROSSES_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring crosses a polygon if at least some of the ring interior is in the +polygon interior and some of the ring interior is exterior to the polygon. =# +_crosses( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + CROSSES_POLYGON_ALLOWS..., + CROSSES_REQUIRES..., + closed_line = true, +) + + +# # Polygons cross geometries + +#= Polygon crosses a curve if at least some of the curve interior is in the +polygon interior and some of the curve interior is exterior to the polygon.=# +_crosses( + trait1::GI.PolygonTrait, g1, + trait2::GI.AbstractCurveTrait, g2 +) = _crosses(trait2, g2, trait1, g1) + + +# # Geometries cross multi-geometry/geometry collections + +#= Geometry crosses a multi-geometry or a collection if the geometry crosses +one of the elements of the collection. =# +function _crosses( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + crosses(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections cross geometries + +#= Multi-geometry or a geometry collection crosses a geometry one elements of +the collection crosses the geometry. =# +function _crosses( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + crosses(sub_g1, g2) && return true + end + return false +end \ No newline at end of file diff --git a/src/methods/geom_relations/disjoint.jl b/src/methods/geom_relations/disjoint.jl new file mode 100644 index 000000000..f3e051c05 --- /dev/null +++ b/src/methods/geom_relations/disjoint.jl @@ -0,0 +1,247 @@ +# # Disjoint + +export disjoint +#= +## What is disjoint? + +The disjoint function checks if one geometry is outside of another geometry, +without sharing any boundaries or interiors. + +To provide an example, consider these two lines: +```@example disjoint +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(2.0, 0.0), (2.75, 0.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 none of the edges or vertices of l1 interact with l2 so they are +disjoint. +```@example disjoint +disjoint(l1, l2) # returns true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`disjoint` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are not allowed to be in the interior of g2 + - points of g1 are not allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - no points required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const DISJOINT_ALLOWS = (in_allow = false, on_allow = false, out_allow = true) +const DISJOINT_CURVE_ALLOWS = (over_allow = false, cross_allow = false, on_allow = false, out_allow = true) +const DISJOINT_REQUIRES = (in_require = false, on_require = false, out_require = false) +""" + disjoint(geom1, geom2)::Bool + +Return `true` if the first geometry is disjoint from the second geometry. + +Return `true` if the first geometry is disjoint from the second geometry. The +interiors and boundaries of both geometries must not intersect. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI + +line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +point = (2, 2) +GO.disjoint(point, line) + +# output +true +``` +""" +disjoint(g1, g2)::Bool = _disjoint(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +disjoint(::FeatureTrait, g1, ::Any, g2)::Bool = disjoint(GI.geometry(g1), g2) +disjoint(::Any, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2)) + + +# # Point disjoint geometries + +# Point is disjoint from another point if the points are not equal. +_disjoint( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = !equals(g1, g2) + +# Point is disjoint from a linestring if it is not on the line's edges/vertices. +_disjoint( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + DISJOINT_ALLOWS..., + closed_curve = false, +) + +# Point is disjoint from a linearring if it is not on the ring's edges/vertices. +_disjoint( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + DISJOINT_ALLOWS..., + closed_curve = true, +) + +#= Point is disjoint from a polygon if it is not on any edges, vertices, or +within the polygon's interior. =# +_disjoint( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., +) + +#= Geometry is disjoint from a point if the point is not in the interior or on +the boundary of the geometry. =# +_disjoint( + trait1::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + trait2::GI.PointTrait, g2, +) = _disjoint(trait2, g2, trait1, g1) + + +# # Lines disjoint geometries + +#= Linestring is disjoint from another line if they do not share any interior +edge/vertex points or boundary points. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is disjoint from a linearring if they do not share any interior +edge/vertex points or boundary points. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is disjoint from a polygon if the interior and boundary points of +the line are not in the polygon's interior or on the polygon's boundary. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, +) + +#= Geometry is disjoint from a linestring if the line's interior and boundary +points don't intersect with the geometrie's interior and boundary points. =# +_disjoint( + trait1::Union{GI.LinearRingTrait, GI.PolygonTrait}, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _disjoint(trait2, g2, trait1, g1) + + +# # Rings disjoint geometries + +#= Linearrings is disjoint from another linearring if they do not share any +interior edge/vertex points or boundary points.=# +_disjoint( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is disjoint from a polygon if the interior and boundary points of +the ring are not in the polygon's interior or on the polygon's boundary. =# +_disjoint( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = true, +) + +# # Polygon disjoint geometries + +#= Polygon is disjoint from another polygon if they do not share any edges or +vertices and if their interiors do not intersect, excluding any holes. =# +_disjoint( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., +) + + +# # Geometries disjoint multi-geometry/geometry collections + +#= Geometry is disjoint from a multi-geometry or a collection if all of the +elements of the collection are disjoint from the geometry. =# +function _disjoint( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + !disjoint(g1, sub_g2) && return false + end + return true +end + +# # Multi-geometry/geometry collections coveredby geometries + +#= Multi-geometry or a geometry collection is covered by a geometry if all +elements of the collection are covered by the geometry. =# +function _disjoint( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !disjoint(sub_g1, g2) && return false + end + return true +end diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl new file mode 100644 index 000000000..3068f023f --- /dev/null +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -0,0 +1,733 @@ +#= Code is based off of DE-9IM Standards (https://en.wikipedia.org/wiki/DE-9IM) +and attempts a standardized solution for most of the functions. +=# + +@enum PointOrientation point_in=1 point_on=2 point_out=3 + +@enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4 + + +#= +Determines if a point meets the given checks with respect to a curve. + +If in_allow is true, the point can be on the curve interior. +If on_allow is true, the point can be on the curve boundary. +If out_allow is true, the point can be disjoint from the curve. + +If the point is in an "allowed" location, return true. Else, return false. + +If closed_curve is true, curve is treated as a closed curve where the first and +last point are connected by a segment. +=# +function _point_curve_process( + point, curve; + in_allow, on_allow, out_allow, + closed_curve = false, +) + # Determine if curve is closed + n = GI.npoint(curve) + first_last_equal = equals(GI.getpoint(curve, 1), GI.getpoint(curve, n)) + closed_curve |= first_last_equal + n -= first_last_equal ? 1 : 0 + # Loop through all curve segments + p_start = GI.getpoint(curve, closed_curve ? n : 1) + @inbounds for i in (closed_curve ? 1 : 2):n + p_end = GI.getpoint(curve, i) + seg_val = _point_segment_orientation(point, p_start, p_end) + seg_val == point_in && return in_allow + if seg_val == point_on + if !closed_curve # if point is on curve endpoints, it is "on" + i == 2 && equals(point, p_start) && return on_allow + i == n && equals(point, p_end) && return on_allow + end + return in_allow + end + p_start = p_end + end + return out_allow +end + +#= +Determines if a point meets the given checks with respect to a polygon. + +If in_allow is true, the point can be within the polygon interior +If on_allow is true, the point can be on the polygon boundary. +If out_allow is true, the point can be disjoint from the polygon. + +If the point is in an "allowed" location, return true. Else, return false. +=# +function _point_polygon_process( + point, polygon; + in_allow, on_allow, out_allow, +) + # Check interaction of geom with polygon's exterior boundary + ext_val = _point_filled_curve_orientation(point, GI.getexterior(polygon)) + # If a point is outside, it isn't interacting with any holes + ext_val == point_out && return out_allow + # if a point is on an external boundary, it isn't interacting with any holes + ext_val == point_on && return on_allow + + # If geom is within the polygon, need to check interactions with holes + for hole in GI.gethole(polygon) + hole_val = _point_filled_curve_orientation(point, hole) + # If a point in in a hole, it is outside of the polygon + hole_val == point_in && return out_allow + # If a point in on a hole edge, it is on the edge of the polygon + hole_val == point_on && return on_allow + end + + # Point is within external boundary and on in/on any holes + return in_allow +end + +#= +Determines if a line meets the given checks with respect to a curve. + +If over_allow is true, segments of the line and curve can be co-linear. +If cross_allow is true, segments of the line and curve can cross. +If on_allow is true, endpoints of either the line or curve can intersect a + segment of the other geometry. +If cross_allow is true, segments of the line and curve can be disjoint. + +If in_require is true, the interiors of the line and curve must meet in at least + one point. +If on_require is true, the bounday of one of the two geometries can meet the + interior or boundary of the other geometry in at least one point. +If out_require is true, there must be at least one point of the given line that + is exterior of the curve. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. Same with closed_curve. +=# +function _line_curve_process( + line, curve; + over_allow, cross_allow, on_allow, out_allow, + in_require, on_require, out_require, + closed_line = false, + closed_curve = false, +) + # Set up requirments + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Determine curve endpoints + nl = GI.npoint(line) + nc = GI.npoint(curve) + first_last_equal_line = equals(GI.getpoint(line, 1), GI.getpoint(line, nl)) + first_last_equal_curve = equals(GI.getpoint(curve, 1), GI.getpoint(curve, nc)) + nl -= first_last_equal_line ? 1 : 0 + nc -= first_last_equal_curve ? 1 : 0 + closed_line |= first_last_equal_line + closed_curve |= first_last_equal_curve + # Loop over each line segment + l_start = GI.getpoint(line, closed_line ? nl : 1) + i = closed_line ? 1 : 2 + while i ≤ nl + l_end = GI.getpoint(line, i) + c_start = GI.getpoint(curve, closed_curve ? nc : 1) + # Loop over each curve segment + for j in (closed_curve ? 1 : 2):nc + c_end = GI.getpoint(curve, j) + # Check if line and curve segments meet + seg_val = _segment_segment_orientation( + (l_start, l_end), + (c_start, c_end), + ) + # If segments are co-linear + if seg_val == line_over + !over_allow && return false + # at least one point in, meets requirments + in_req_met = true + point_val = _point_segment_orientation( + l_start, + c_start, c_end, + ) + # If entire segment isn't covered, consider remaining section + if point_val != point_out + if _point_segment_orientation( + l_end, + c_start, c_end, + ) != point_out + l_start = l_end + i += 1 + break + elseif _point_segment_orientation( + c_start, + l_start, l_end, + ) != point_out && !equals(l_start, c_start) + l_start = c_start + break + elseif _point_segment_orientation( + c_end, + l_start, l_end, + ) != point_out && !equals(l_start, c_end) + l_start = c_end + break + end + end + else + if seg_val == line_cross + !cross_allow && return false + in_req_met = true + elseif seg_val == line_hinge # could cross or overlap + # Determine location of intersection point on each segment + _, fracs = _intersection_point( + (_tuple_point(l_start), _tuple_point(l_end)), + (_tuple_point(c_start), _tuple_point(c_end)) + ) + (α, β) = + if !isnothing(fracs) + fracs + else # line and curve segments are parallel + if equals(l_start, c_start) + (0, 0) + elseif equals(l_start, c_end) + (0, 1) + elseif equals(l_end, c_start) + (1, 0) + else # equals(l_end, c_end) + (1, 1) + end + end + + if ( # don't consider edges of curves as they can't cross + (β == 0 && !closed_curve && j == 2) || + (β == 1 && !closed_curve && j == nc) || + (α == 0 && !closed_line && i == 2) || + (α == 1 && !closed_line && i == nl) + ) + !on_allow && return false + on_req_met = true + else + in_req_met = true + # If needed, determine if hinge actually crosses + if (!cross_allow || !over_allow) && α != 0 && β != 0 + l, c = if β == 1 + if α == 1 + ( + (l_end, GI.getpoint(line, i + 1)), + (c_end, GI.getpoint(curve, j + 1)) + ) + else + ( + (l_start, l_end), + (c_end, GI.getpoint(curve, j + 1)) + ) + end + else # β ≠ 1 and α == 1 + ( + (l_end, GI.getpoint(line, i + 1)), + (c_start, c_end) + ) + end + if _segment_segment_orientation(l, c) == line_hinge + !cross_allow && return false + else + !over_allow && return false + end + end + end + end + # no overlap for a give segment + if j == nc + !out_allow && return false + out_req_met = true + end + end + c_start = c_end + if j == nc + i += 1 + l_start = l_end + end + end + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a line meets the given checks with respect to a polygon. + +If in_allow is true, segments of the line can be in the polygon interior. +If on_allow is true, segments of the line can be on the polygon's boundary. +If out_allow is true, segments of the line can be outside of the polygon. + +If in_require is true, the interiors of the line and polygon must meet in at + least one point. +If on_require is true, the line must have at least one point on the polygon'same + boundary. +If out_require is true, the line must have at least one point outside of the + polygon. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. +=# +function _line_polygon_process( + line, polygon; + in_allow, on_allow, out_allow, + in_require, on_require, out_require, + closed_line = false, +) + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Check interaction of line with polygon's exterior boundary + in_curve, on_curve, out_curve = _line_filled_curve_interactions( + line, GI.getexterior(polygon); + closed_line = closed_line, + ) + if on_curve + !on_allow && return false + on_req_met = true + end + if out_curve + !out_allow && return false + out_req_met = true + end + # If no points within the polygon, the line is disjoint and we are done + !in_curve && return in_req_met && on_req_met && out_req_met + + # Loop over polygon holes + for hole in GI.gethole(polygon) + in_hole, on_hole, out_hole =_line_filled_curve_interactions( + line, hole; + closed_line = closed_line, + ) + if in_hole # line in hole is equivalent to being out of polygon + !out_allow && return false + out_req_met = true + end + if on_hole # hole bounday is polygon boundary + !on_allow && return false + on_req_met = true + end + if !out_hole # entire line is in/on hole, can't be in/on other holes + in_curve = false + break + end + end + if in_curve # entirely of curve isn't within a hole + !in_allow && return false + in_req_met = true + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a polygon meets the given checks with respect to a polygon. + +If in_allow is true, the polygon's interiors must intersect. +If on_allow is true, the one of the polygon's boundaries must either interact + with the other polygon's boundary or interior. +If out_allow is true, the first polygon must have interior regions outside of + the second polygon. + +If in_require is true, the polygon interiors must meet in at least one point. +If on_require is true, one of the polygon's must have at least one boundary + point in or on the other polygon. +If out_require is true, the first polygon must have at least one interior point + outside of the second polygon. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. +=# +function _polygon_polygon_process( + poly1, poly2; + in_allow, on_allow, out_allow, + in_require, on_require, out_require, +) + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Check if exterior of poly1 is within poly2 + ext1 = GI.getexterior(poly1) + ext2 = GI.getexterior(poly2) + # Check if exterior of poly1 is in polygon 2 + e1_in_p2, e1_on_p2, e1_out_p2 = _line_polygon_interactions( + ext1, poly2; + closed_line = true, + ) + if e1_on_p2 + !on_allow && return false + on_req_met = true + end + if e1_out_p2 + !out_allow && return false + out_req_met = true + end + + if !e1_in_p2 + # if exterior ring isn't in poly2, check if it surrounds poly2 + _, _, e2_out_e1 = _line_filled_curve_interactions( + ext2, ext1; + closed_line = true, + ) # if they really are disjoint, we are done + e2_out_e1 && return in_req_met && on_req_met && out_req_met + end + # If interiors interact, check if poly2 interacts with any of poly1's holes + for h1 in GI.gethole(poly1) + h1_in_p2, h1_on_p2, h1_out_p2 = _line_polygon_interactions( + h1, poly2; + closed_line = true, + ) + if h1_on_p2 + !on_allow && return false + on_req_met = true + end + if h1_out_p2 + !out_allow && return false + out_req_met = true + end + if !h1_in_p2 + # If hole isn't in poly2, see if poly2 is in hole + _, _, e2_out_h1 = _line_filled_curve_interactions( + ext2, h1; + closed_line = true, + ) + # hole encompasses all of poly2 + !e2_out_h1 && return in_req_met && on_req_met && out_req_met + break + end + end + #= + poly2 isn't outside of poly1 and isn't in a hole, poly1 interior must + interact with poly2 interior + =# + !in_allow && return false + in_req_met = true + + # If any of poly2 holes are within poly1, part of poly1 is exterior to poly2 + for h2 in GI.gethole(poly2) + h2_in_p1, h2_on_p1, _ = _line_polygon_interactions( + h2, poly1; + closed_line = true, + ) + if h2_on_p1 + !on_allow && return false + on_req_met = true + end + if h2_in_p1 + !out_allow && return false + out_req_met = true + end + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a point is in, on, or out of a segment. If the point is `on` the +segment it is on one of the segments endpoints. If it is `in`, it is on any +other point of the segment. If the point is not on any part of the segment, it +is `out` of the segment. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. +=# +function _point_segment_orientation( + point, start, stop; + in::T = point_in, on::T = point_on, out::T = point_out, +) where {T} + # 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) + Δx_seg = x2 - x1 + Δy_seg = y2 - y1 + Δx_pt = x - x1 + Δy_pt = y - y1 + if (Δx_pt == 0 && Δy_pt == 0) || (Δx_pt == Δx_seg && Δy_pt == Δy_seg) + # If point is equal to the segment start or end points + return on + else + #= + Determine if the point is on the segment -> see if vector from segment + start to point is parallel to segment and if point is between the + segment endpoints + =# + on_line = _isparallel(Δx_seg, Δy_seg, Δx_pt, Δy_pt) + !on_line && return out + between_endpoints = + (x2 > x1 ? x1 <= x <= x2 : x2 <= x <= x1) && + (y2 > y1 ? y1 <= y <= y2 : y2 <= y <= y1) + !between_endpoints && return out + end + return in +end + +#= +Determine if point is in, on, or out of a closed curve, which includes the space +enclosed by the closed curve. + +`In` means the point is within the closed curve (excluding edges and vertices). +`On` means the point is on an edge or a vertex of the closed curve. +`Out` means the point is outside of the closed curve. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait, that is assumed to be closed, regardless of +repeated last point. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. + +Note that this uses the Algorithm by Hao and Sun (2018): +https://doi.org/10.3390/sym10100477 +Paper seperates orientation of point and edge into 26 cases. For each case, it +is either a case where the point is on the edge (returns on), where a ray from +the point (x, y) to infinity along the line y = y cut through the edge (k += 1), +or the ray does not pass through the edge (do nothing and continue). If the ray +passes through an odd number of edges, it is within the curve, else outside of +of the curve if it didn't return 'on'. +See paper for more information on cases denoted in comments. +=# +function _point_filled_curve_orientation( + point, curve; + in::T = point_in, on::T = point_on, out::T = point_out, +) where {T} + x, y = GI.x(point), GI.y(point) + n = GI.npoint(curve) + n -= equals(GI.getpoint(curve, 1), GI.getpoint(curve, n)) ? 1 : 0 + k = 0 # counter for ray crossings + p_start = GI.getpoint(curve, n) + @inbounds for i in 1:n + p_end = GI.getpoint(curve, i) + v1 = GI.y(p_start) - y + v2 = GI.y(p_end) - y + if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26 + u1 = GI.x(p_start) - x + u2 = GI.x(p_end) - x + c1 = u1 * v2 # first element of cross product summation + c2 = u2 * v1 # second element of cross product summation + f = c1 - c2 + if v2 > 0 && v1 ≤ 0 # Case 3, 9, 16, 21, 13, or 24 + (c1 ≈ c2) && return on # Case 16 or 21 + f > 0 && (k += 1) # Case 3 or 9 + elseif v1 > 0 && v2 ≤ 0 # Case 4, 10, 19, 20, 12, or 25 + (c1 ≈ c2) && return on # Case 19 or 20 + f < 0 && (k += 1) # Case 4 or 10 + elseif v2 == 0 && v1 < 0 # Case 7, 14, or 17 + (c1 ≈ c2) && return on # Case 17 + elseif v1 == 0 && v2 < 0 # Case 8, 15, or 18 + (c1 ≈ c2) && return on # Case 18 + elseif v1 == 0 && v2 == 0 # Case 1, 2, 5, 6, 22, or 23 + u2 ≤ 0 && u1 ≥ 0 && return on # Case 1 + u1 ≤ 0 && u2 ≥ 0 && return on # Case 2 + end + end + p_start = p_end + end + return iseven(k) ? out : in +end + +#= +Determines the type of interaction between two line segments. If the segments +`cross`, this means that they have a single intersection point that isn't on +either of their enpoints. If they form a `hinge`, they meet at one of the +segments endpoints. If they are `over`, then they are co-linear for at least +some of the length of the segments. Finally, if they are `out`, then the +segments are disjoint. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. +=# +function _segment_segment_orientation( + (a_point, b_point), (c_point, d_point); + cross::T = line_cross, hinge::T = line_hinge, + over::T = line_over, out::T = line_out, +) where T + (ax, ay) = _tuple_point(a_point) + (bx, by) = _tuple_point(b_point) + (cx, cy) = _tuple_point(c_point) + (dx, dy) = _tuple_point(d_point) + meet_type = ExactPredicates.meet((ax, ay), (bx, by), (cx, cy), (dx, dy)) + # Lines meet at one point within open segments + meet_type == 1 && return cross + # Lines don't meet at any points + meet_type == -1 && return out + # Lines meet at one or more points within closed segments + if _isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy))) + min_x, max_x = cx < dx ? (cx, dx) : (dx, cx) + min_y, max_y = cy < dy ? (cy, dy) : (dy, cy) + if ( + ((ax ≤ min_x && bx ≤ min_x) || (ax ≥ max_x && bx ≥ max_x)) && + ((ay ≤ min_y && by ≤ min_y) || (ay ≥ max_y && by ≥ max_y)) + ) + # a_point and b_point are on the same side of segment, don't overlap + return hinge + else + return over + end + end + # if lines aren't parallel then they must hinge + return hinge +end + +#= +Determines if a line meets the given checks with respect to a filled-in curve. +By filled-in curve, I am referring to the exterior ring of a poylgon, for +example. + +If over_allow is true, segments of the line and curve can be co-linear. +If cross_allow is true, segments of the line and curve can cross. +If on_allow is true, endpoints of either the line or curve can intersect a + segment of the other geometry. +If cross_allow is true, segments of the line and curve can be disjoint. + +If in_require is true, the interiors of the line and curve must meet in at least + one point. +If on_require is true, the bounday of one of the two geometries can meet the + interior or boundary of the other geometry in at least one point. +If out_require is true, there must be at least one point of the given line that + is exterior of the curve. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. Same with closed_curve. +=# +function _line_filled_curve_interactions( + line, curve; + closed_line = false, +) + in_curve = false + on_curve = false + out_curve = false + + # Determine number of points in curve and line + nl = GI.npoint(line) + nc = GI.npoint(curve) + first_last_equal_line = equals(GI.getpoint(line, 1), GI.getpoint(line, nl)) + first_last_equal_curve = equals(GI.getpoint(curve, 1), GI.getpoint(curve, nc)) + nl -= first_last_equal_line ? 1 : 0 + nc -= first_last_equal_curve ? 1 : 0 + closed_line |= first_last_equal_line + + # See if first point is in an acceptable orientation + l_start = GI.getpoint(line, closed_line ? nl : 1) + point_val = _point_filled_curve_orientation(l_start, curve) + if point_val == point_in + in_curve = true + elseif point_val == point_on + on_curve = true + else # point_val == point_out + out_curve = true + end + + # Check for any intersections between line and curve + for i in (closed_line ? 1 : 2):nl + l_end = GI.getpoint(line, i) + c_start = GI.getpoint(curve, nc) + # If already interacted with all regions of curve, can stop + in_curve && on_curve && out_curve && break + # Check next segment of line against curve + for j in 1:nc + c_end = GI.getpoint(curve, j) + # Check if two line and curve segments meet + seg_val = _segment_segment_orientation( + (l_start, l_end), + (c_start, c_end), + ) + if seg_val != line_out + # If line and curve meet, then at least one point is on boundary + on_curve = true + if seg_val == line_cross + # When crossing boundary, line is both in and out of curve + in_curve = true + out_curve = true + else + if seg_val == line_over + sp = _point_segment_orientation(l_start, c_start, c_end) + lp = _point_segment_orientation(l_end, c_start, c_end) + if sp != point_in || lp != point_in + #= + Line crosses over segment endpoint, creating a hinge + with another segment. + =# + seg_val = line_hinge + end + end + if seg_val == line_hinge + #= + Can't determine all types of interactions (in, out) with + hinge as it could pass through multiple other segments + so calculate if segment endpoints and intersections are + in/out of filled curve + =# + ipoints = intersection_points( + GI.Line([l_start, l_end]), + curve + ) + npoints = length(ipoints) # since hinge, at least one + sort!(ipoints, by = p -> _euclid_distance(Float64, p, l_start)) + p_start = _tuple_point(l_start) + for i in 1:(npoints + 1) + p_end = i ≤ npoints ? + ipoints[i] : + _tuple_point(l_end) + mid_val = _point_filled_curve_orientation( + (p_start .+ p_end) ./ 2, + curve, + ) + if mid_val == point_in + in_curve = true + elseif mid_val == point_out + out_curve = true + end + end + # already checked segment against whole filled curve + l_start = l_end + break + end + end + end + c_start = c_end + end + l_start = l_end + end + return in_curve, on_curve, out_curve +end + +function _line_polygon_interactions( + line, polygon; + closed_line = false, +) + in_ext, on_ext, out_ext = _line_filled_curve_interactions( + line, GI.getexterior(polygon); + closed_line = closed_line, + ) + !in_ext && return (in_ext, on_ext, out_ext) + # Loop over polygon holes + for hole in GI.gethole(polygon) + in_hole, on_hole, out_hole =_line_filled_curve_interactions( + line, hole; + closed_line = closed_line, + ) + if in_hole + out_ext = true + end + if on_hole + on_ext = true + end + if !out_hole # entire line is in/on hole, can't be in/on other holes + in_ext = false + return (in_ext, on_ext, out_ext) + end + end + return in_ext, on_ext, out_ext +end + +function _point_in_extent(p, extent::Extents.Extent) + (x1, x2), (y1, y2) = extent.X, extent.Y + return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 +end \ No newline at end of file diff --git a/src/methods/intersects.jl b/src/methods/geom_relations/intersects.jl similarity index 75% rename from src/methods/intersects.jl rename to src/methods/geom_relations/intersects.jl index 2efaf1b78..834745a33 100644 --- a/src/methods/intersects.jl +++ b/src/methods/geom_relations/intersects.jl @@ -5,7 +5,10 @@ export intersects, intersection, intersection_points #= ## What is `intersects` vs `intersection` vs `intersection_points`? -The `intersects` methods check whether two geometries intersect with each other. +The intersects function checks if a given geometry intersects with another +geometry, or in other words, the either the interiors or boundaries of the two +geometries intersect. + 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. @@ -44,18 +47,16 @@ f 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. +Given that intersects is the exact opposite of disjoint, we simply pass the two +inputs variables, swapped in order, to disjoint. =# """ intersects(geom1, geom2)::Bool -Check if two geometries intersect, returning true if so and false otherwise. +Return true if the interiors or boundaries of the two geometries interact. + +`intersects` returns the exact opposite result of `disjoint`. ## Example @@ -70,77 +71,8 @@ GO.intersects(line1, line2) true ``` """ -intersects(geom1, geom2) = intersects( - GI.trait(geom1), - geom1, - GI.trait(geom2), - geom2 -) - -""" - intersects(::GI.LineTrait, a, ::GI.LineTrait, b)::Bool - -Returns true if two line segments intersect and false otherwise. -""" -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 == 0 || meet_type == 1 -end - -""" - 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, convert the geometry to a list of edges -and cross compare the edges for intersections. -""" -function intersects( - 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 +intersects(geom1, geom2) = !disjoint(geom1, geom2) -""" - _line_intersects( - edges_a::Vector{Edge}, - edges_b::Vector{Edge} - )::Bool - -Returns true if there is at least one intersection between edges within the -two lists of edges. -""" -function _line_intersects( - edges_a::Vector{Edge}, - 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 - _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} @@ -229,8 +161,8 @@ 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, + trait_a::GI.AbstractGeometryTrait, geom_a, + trait_b::GI.AbstractGeometryTrait, geom_b, ) @assert( false, @@ -269,7 +201,7 @@ 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) +function intersection_points(::GI.AbstractGeometryTrait, a, ::GI.AbstractGeometryTrait, 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 diff --git a/src/methods/geom_relations/old_code.jl b/src/methods/geom_relations/old_code.jl new file mode 100644 index 000000000..e69de29bb diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl new file mode 100644 index 000000000..9133abf63 --- /dev/null +++ b/src/methods/geom_relations/overlaps.jl @@ -0,0 +1,241 @@ +# # Overlaps + +export overlaps + +#= +## What is overlaps? + +The overlaps function checks if two geometries overlap. Two geometries overlap +if they have the same dimension, and if they overlap then their interiors +interact, but they both also need interior points exterior to the other +geometry. + +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 (crosses don't count for interior points interacting). + +To provide an example, consider these two lines: +```@example overlaps +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 overlaps +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. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`overlaps` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - at least one point of g1 is required to be in the interior of g2 + - at least one point of g2 is required to be in the interior of g1 + - no points of g1 is required to be on the boundary of g2 + - at least one point of g1 is required to be in the exterior of g2 + - at least one point of g2 is required to be in the exterior of g1 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const OVERLAPS_CURVE_ALLOWS = (over_allow = true, cross_allow = false, on_allow = true, out_allow = true) +const OVERLAPS_POLYGON_ALLOWS = (in_allow = true, on_allow = true, out_allow = true) +const OVERLAPS_REQUIRES = (in_require = true, on_require = false, out_require = true) + +""" + overlaps(geom1, geom2)::Bool + +Compare two Geometries of the same dimension and return true if their interiors +interact, but they both also have interior points exterior to the other +geometry. Lines crossing doesn't count for interiors interacting as overlaps +of curves must be of dimension one. + +## 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([[(1,1), (1,6), (6,6), (6,1), (1,1)]]) + +GO.overlaps(poly1, poly2) +# output +true +``` +""" +overlaps(g1, g2)::Bool = _overlaps(GI.trait(g1), g1, GI.trait(g2), g2) + + +# # Convert features to geometries +_overlaps(::GI.FeatureTrait, g1, ::Any, g2) = overlaps(GI.geometry(g1), g2) +_overlaps(::Any, g1, t2::GI.FeatureTrait, g2) = overlaps(g1, GI.geometry(g2)) + + +# # Non-specified geometries + +# Geometries of different dimensions and points cannot overlap and return false +_overlaps( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g2, +) = false + + +# # Lines cross curves + +#= Linestring overlaps with another linestring when they share co-linear +segments (interiors interacting), but both have interior points exterior to the +other line. =# +_overlaps( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = false, + closed_curve = false, +) && _line_curve_process( + g2, g1; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = false, + closed_curve = false, + ) + +#= Linestring overlaps with a linearring when they share co-linear segments +(interiors interacting), but both have interior points exterior to the other. =# +_overlaps( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = false, + closed_curve = true, +) && _line_curve_process( + g2, g1; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = true, + closed_curve = false, + ) + + +# # Rings cross curves + +#= Linearring overlaps with a linestring when they share co-linear segments +(interiors interacting), but both have interior points exterior to the other. =# +_overlaps( + trait1::GI.LinearRingTrait, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _overlaps(trait2, g2, trait1, g1) + +#= Linearring overlaps with another linearring when they share co-linear +segments (interiors interacting), but both have interior points exterior to the +other line. =# +_overlaps( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = true, + closed_curve = true, +) && _line_curve_process( + g2, g1; + OVERLAPS_CURVE_ALLOWS..., + OVERLAPS_REQUIRES..., + closed_line = true, + closed_curve = true, + ) + + +# # Polygons cross polygons + +#= Polygon overlaps with another polygon when their interiors intersect, but +both have interior points exterior to the other polygon. =# +_overlaps( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + OVERLAPS_POLYGON_ALLOWS..., + OVERLAPS_REQUIRES..., +) && _polygon_polygon_process( + g2, g1; + OVERLAPS_POLYGON_ALLOWS..., + OVERLAPS_REQUIRES..., +) + +# # Geometries disjoint multi-geometry/geometry collections + +# Multipoints overlap with other multipoints if only some sub-points are shared +function _overlaps( + ::GI.MultiPointTrait, g1, + ::GI.MultiPointTrait, g2, +) + 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(g1) + match_point = false + for p2 in GI.getpoint(g2) + 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 + +#= Geometry overlaps a multi-geometry or a collection if the geometry overlaps +at least one of the elements of the collection. =# +function _overlaps( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + overlaps(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections cross geometries + +#= Multi-geometry or a geometry collection overlaps a geometry if at least one +elements of the collection overlaps the geometry. =# +function _overlaps( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + overlaps(sub_g1, g2) && return true + end + return false +end diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl new file mode 100644 index 000000000..24bf4b2ea --- /dev/null +++ b/src/methods/geom_relations/touches.jl @@ -0,0 +1,250 @@ +# # Touches + +export touches + +#= +## What is touches? + +The touches function checks if one geometry touches another geometry. In other +words, the interiors of the two geometries don't interact, but one of the +geometries must have a boundary point that interacts with either the other +geometies interior or boundary. + + +To provide an example, consider these two lines: +```@example touches +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = Line([Point(0.0, 0.0), Point(1.0, 0.0)]) +l2 = Line([Point(1.0, 0.0), Point(1.0, -1.0)]) + +f, a, p = lines(l1) +lines!(l2) +``` +We can see that these two lines touch only at their endpoints. +```@example touches +touches(l1, l2) # true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`touches` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are not allowed to be in the interior of g2 + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - no points of g1 are required to be in the interior of g2 + - at least one point of g1 is required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const TOUCHES_POINT_ALLOWED = (in_allow = false, on_allow = true, out_allow = false) +const TOUCHES_CURVE_ALLOWED = (over_allow = false, cross_allow = false, on_allow = true, out_allow = true) +const TOUCHES_POLYGON_ALLOWS = (in_allow = false, on_allow = true, out_allow = true) +const TOUCHES_REQUIRES = (in_require = false, on_require = true, out_require = false) + +""" + touches(geom1, geom2)::Bool + +Return `true` if the first geometry touches the second geometry. In other words, +the two interiors cannot interact, but one of the geometries must have a +boundary point that interacts with either the other geometies interior or +boundary. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI + +l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)]) +l2 = GI.Line([(0.5, 1.0), (0.5, -1.0)]) + +GO.crosses(l1, l2) +# output +true +``` +""" +touches(g1, g2)::Bool = _touches(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_touches(::GI.FeatureTrait, g1, ::Any, g2) = touches(GI.geometry(g1), g2) +_touches(::Any, g1, t2::GI.FeatureTrait, g2) = touches(g1, GI.geometry(g2)) + + +# # Point touches geometries + +# Point cannot touch another point as if they are equal, interiors interact +_touches( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = false + +# Point touches a linestring if it equal to the first of last point of the line +function _touches( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) + n = GI.npoint(g2) + p1 = GI.getpoint(g2, 1) + pn = GI.getpoint(g2, n) + equals(p1, pn) && return false + return equals(g1, p1) || equals(g1, pn) +end + +# Point cannot 'touch' a linearring given that the ring has no boundary points +_touches( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = false + +# Point touches a polygon if it is on the boundary of that polygon +_touches( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + TOUCHES_POINT_ALLOWED..., +) + +#= Geometry touches a point if the point is on the geometry boundary. =# +_touches( + trait1::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + trait2::GI.PointTrait, g2, +) = _touches(trait2, g2, trait1, g1) + + +# # Lines touching geometries + +#= Linestring touches another line if at least one bounday point interacts with +the bounday of interior of the other line, but the interiors don't interact. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + TOUCHES_CURVE_ALLOWED..., + TOUCHES_REQUIRES..., + closed_line = false, + closed_curve = false, +) + + +#= Linestring touches a linearring if at least one of the boundary points of the +line interacts with the linear ring, but their interiors can't interact. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + TOUCHES_CURVE_ALLOWED..., + TOUCHES_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring touches a polygon if at least one of the boundary points of the +line interacts with the boundary of the polygon. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., + closed_line = false, +) + + +# # Rings touch geometries + +#= Linearring touches a linestring if at least one of the boundary points of the +line interacts with the linear ring, but their interiors can't interact. =# +_touches( + trait1::GI.LinearRingTrait, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _touches(trait2, g2, trait1, g1) + +#= Linearring cannot touch another linear ring since they are both exclusively +made up of interior points and no bounday points =# +_touches( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = false + +#= Linearring touches a polygon if at least one of the points of the ring +interact with the polygon bounday and non are in the polygon interior. =# +_touches( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., + closed_line = true, +) + + +# # Polygons touch geometries + +#= Polygon touches a curve if at least one of the curve bounday points interacts +with the polygon's bounday and no curve points interact with the interior.=# +_touches( + trait1::GI.PolygonTrait, g1, + trait2::GI.AbstractCurveTrait, g2 +) = _touches(trait2, g2, trait1, g1) + + +#= Polygon touches another polygon if they share at least one boundary point and +no interior points. =# +_touches( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., +) + +# # Geometries touch multi-geometry/geometry collections + +#= Geometry touch a multi-geometry or a collection if the geometry touches at +least one of the elements of the collection. =# +function _touches( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + !touches(g1, sub_g2) && return false + end + return true +end + +# # Multi-geometry/geometry collections cross geometries + +#= Multi-geometry or a geometry collection touches a geometry if at least one +elements of the collection touches the geometry. =# +function _touches( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !touches(sub_g1, g2) && return false + end + return true +end \ No newline at end of file diff --git a/src/methods/geom_relations/within.jl b/src/methods/geom_relations/within.jl new file mode 100644 index 000000000..9b28c0ee7 --- /dev/null +++ b/src/methods/geom_relations/within.jl @@ -0,0 +1,269 @@ +# # Within + +export within + +#= +## What is within? + +The within function checks if one geometry is inside another geometry. This +requires that the two interiors intersect and that the interior and +boundary of the first geometry is not in the exterior of the second geometry. + +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), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(0.25, 0.0), (0.75, 0.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 all of the points and edges of l2 are within l1, so l2 is +within l1, but l1 is not within l2 +```@example cshape +within(l1, l2) # returns false +within(l2, l1) # returns true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`within` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 (either through + overlap or crossing for lines) + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are not allowed to be in the exterior of g2 + - at least one point of g1 is required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const WITHIN_POINT_ALLOWS = (in_allow = true, on_allow = false, out_allow = false) +const WITHIN_CURVE_ALLOWS = (over_allow = true, cross_allow = true, on_allow = true, out_allow = false) +const WITHIN_POLYGON_ALLOWS = (in_allow = true, on_allow = true, out_allow = false) +const WITHIN_REQUIRES = (in_require = true, on_require = false, out_require = false) + +""" + within(geom1, geom2)::Bool + +Return `true` if the first geometry is completely within the second geometry. +The interiors of both geometries must intersect and the interior and boundary of +the primary geometry (geom1) must not intersect the exterior of the secondary +geometry (geom2). + +Furthermore, `within` returns the exact opposite result of `contains`. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI + +line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +point = (1, 2) +GO.within(point, line) + +# output +true +``` +""" +within(g1, g2) = _within(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_within(::GI.FeatureTrait, g1, ::Any, g2) = within(GI.geometry(g1), g2) +_within(::Any, g1, t2::GI.FeatureTrait, g2) = within(g1, GI.geometry(g2)) + + +# # Points within geometries + +# Point is within another point if those points are equal. +_within( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = equals(g1, g2) + +#= Point is within a linestring if it is on a vertex or an edge of that line, +excluding the start and end vertex if the line is not closed. =# +_within( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + WITHIN_POINT_ALLOWS..., + closed_curve = false, +) + +# Point is within a linearring if it is on a vertex or an edge of that ring. +_within( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + WITHIN_POINT_ALLOWS..., + closed_curve = true, +) + +#= Point is within a polygon if it is inside of that polygon, excluding edges, +vertices, and holes. =# +_within( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + WITHIN_POINT_ALLOWS..., +) + +# No geometries other than points can be within points +_within( + ::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::GI.PointTrait, g2, +) = false + + +# # Lines within geometries + +#= Linestring is within another linestring if their interiors intersect and no +points of the first line are in the exterior of the second line. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is within a linear ring if their interiors intersect and no points +of the line are in the exterior of the ring. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is within a polygon if their interiors intersect and no points of +the line are in the exterior of the polygon, although they can be on an edge. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, +) + + +# # Rings covered by geometries + +#= Linearring is within a linestring if their interiors intersect and no points +of the ring are in the exterior of the line. =# +_within( + ::GI.LinearRingTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, + closed_curve = false, +) + +#= Linearring is within another linearring if their interiors intersect and no +points of the first ring are in the exterior of the second ring. =# +_within( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is within a polygon if their interiors intersect and no points of +the ring are in the exterior of the polygon, although they can be on an edge. =# +_within( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, +) + + +# # Polygons within geometries + +#= Polygon is within another polygon if the interior of the first polygon +intersects with the interior of the second and no points of the first polygon +are outside of the second polygon. =# +_within( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., +) + +# Polygons cannot be within any curves +_within( + ::GI.PolygonTrait, g1, + ::GI.AbstractCurveTrait, g2, +) = false + + +# # Geometries within multi-geometry/geometry collections + +#= Geometry is within a multi-geometry or a collection if the geometry is within +at least one of the collection elements. =# +function _within( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + within(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections within geometries + +#= Multi-geometry or a geometry collection is within a geometry if all +elements of the collection are within the geometry. =# +function _within( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !within(sub_g1, g2) && return false + end + return true +end diff --git a/src/methods/orientation.jl b/src/methods/orientation.jl new file mode 100644 index 000000000..8a97dfb85 --- /dev/null +++ b/src/methods/orientation.jl @@ -0,0 +1,43 @@ +# """ +# isparallel(line1::LineString, line2::LineString)::Bool + +# Return `true` if each segment of `line1` is parallel to the correspondent segment of `line2` + +# ## Examples +# ```julia +# import GeoInterface as GI, GeometryOps as GO +# julia> line1 = GI.LineString([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)]) +# GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)], nothing, nothing) + +# julia> line2 = GI.LineString([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)]) +# GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)], nothing, nothing) + +# julia> +# GO.isparallel(line1, line2) +# true +# ``` +# """ +# function isparallel(line1, line2)::Bool +# seg1 = linesegment(line1) +# seg2 = linesegment(line2) + +# for i in eachindex(seg1) +# coors2 = nothing +# coors1 = seg1[i] +# coors2 = seg2[i] +# _isparallel(coors1, coors2) == false && return false +# end +# return true +# end + +# @inline function _isparallel(p1, p2) +# slope1 = bearing_to_azimuth(rhumb_bearing(GI.x(p1), GI.x(p2))) +# slope2 = bearing_to_azimuth(rhumb_bearing(GI.y(p1), GI.y(p2))) + +# return slope1 === slope2 +# end + +_isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy))) = + _isparallel(bx - ax, by - ay, dx - cx, dy - cy) + +_isparallel(Δx1, Δy1, Δx2, Δy2) = (Δx1 * Δy2 == Δy1 * Δx2) \ No newline at end of file diff --git a/src/methods/overlaps.jl b/src/methods/overlaps.jl deleted file mode 100644 index f99b75c9d..000000000 --- a/src/methods/overlaps.jl +++ /dev/null @@ -1,233 +0,0 @@ -# # Overlaps - -export overlaps - -#= -## What is overlaps? - -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. This means one geometry cannot be within or contain the other and -they cannot be equal - -## 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([[(1,1), (1,6), (6,6), (6,1), (1,1)]]) - -GO.overlaps(poly1, poly2) -# output -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 - -""" - 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 - -""" - 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 - -""" - 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/src/methods/within.jl b/src/methods/within.jl deleted file mode 100644 index 16366f944..000000000 --- a/src/methods/within.jl +++ /dev/null @@ -1,43 +0,0 @@ -# # Containment/withinness - -export within - - -""" - within(geom1, geom)::Bool - -Return `true` if the first geometry is completely within the second geometry. -The interiors of both geometries must intersect and, the interior and boundary of the primary (geometry a) -must not intersect the exterior of the secondary (geometry b). -`within` returns the exact opposite result of `contains`. - -## Examples -```jldoctest setup=:(using GeometryOps, GeometryBasics) -import GeometryOps as GO, GeoInterface as GI - -line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) -point = (1, 2) -GO.within(point, line) - -# output -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(::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.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 cb1ff945c..ef86f6c58 100644 --- a/test/methods/bools.jl +++ b/test/methods/bools.jl @@ -1,117 +1,30 @@ -using Test, GeometryOps -import GeoInterface as GI -import GeometryOps as GO - -@testset "booleans" begin - line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]]) - line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]]) - line3 = GI.LineString([ - (-111.544189453125, 24.186847428521244), - (-110.687255859375, 24.966140159912975), - (-110.4510498046875, 24.467150664739002), - (-109.9951171875, 25.180087808990645) - ]) - line4 = GI.LineString([ - (-111.4617919921875, 24.05148034322011), - (-110.8795166015625, 24.681961205014595), - (-110.841064453125, 24.14174098050432), - (-109.97863769531249, 24.617057340809524) - ]) - - # @test isparallel(line1, line2) == true - # @test isparallel(line3, line4) == false - - poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]]) - poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]]) - - @test isconcave(poly1) == true - @test isconcave(poly2) == false - - l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]]) - l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]]) - - @test isclockwise(l1) == true - @test isclockwise(l2) == false - - l3 = GI.LineString([[0, 0], [3, 3], [4, 4]]) - p1 = GI.Point([1,1]) - - l4 = GI.LineString([[0, 0], [3, 3]]) - p2 = GI.Point([0, 0]) - - p3 = GI.Point([20, 20]) - l5 = GI.LineString([[0, 0], [3, 3], [38.32, 5.96]]) - - @test GO.point_on_line(p2, l4; ignore_end_vertices=true) == false - @test GO.point_on_line(p3, l5; ignore_end_vertices=true) == false - @test GO.point_on_line(p1, l3) == true - - pt = (-77, 44) - poly = GI.Polygon([[[-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]]) - - @test point_in_polygon(pt, poly) == true - - poly3 = GI.Polygon([[(1, 1), (1, 10), (10, 10), (10, 1), (1, 1)]]) - poly4 = GI.Polygon([[(1, 1), (2, 2), (3, 2), (1, 1)]]) - line5 = GI.LineString([(1.0, 1.0), (2.0, 3.0), (2.0, 3.5)]) - - line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)]) - poly5 = GI.Polygon([[(1.0, 1.0), (1.0, 20.0), (1.0, 3.0), (1.0, 4.0), (1.0, 1.0)]]) - line7 = GI.LineString([(1.0, 2.0), (1.0, 3.0), (1.0, 3.5)]) - - @test GO.contains(poly3, poly4) == true - @test GO.contains(poly3, line5) == true - @test GO.contains(line6, (1, 2)) == true - @test GO.contains(poly3, poly5) == false - @test GO.contains(poly3 , line7) == false - - @test GO.within(poly4, poly3) == true - @test GO.within(line5, poly3) == true - @test GO.within(poly5, poly3) == false - @test GO.within((1, 2), line6) == true - @test GO.within(line7, poly3) == false - - poly6 = GI.Polygon([[(-11, -12), (-13, -12), (-13, -13), (-11, -13), (-11, -12)]]) - poly7 = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]]) - poly8 = GI.Polygon([[(-1, 2), (-13, -12), (-13, -13), (-11, -13), (-1, 2)]]) - - @test GO.disjoint(poly7, poly6) == true - @test GO.disjoint(poly7, (1, 1)) == true - @test GO.disjoint(poly7, GI.LineString([(0, 0), (12, 2), (12, 3), (12, 4)])) == true - @test GO.disjoint(poly8, poly7) == false - - 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)) - - line10 = GI.LineString([ - (142.03125, -11.695273), - (138.691406, -16.804541), - (136.40625, -14.604847), - (135.966797, -12.039321), - (131.308594, -11.436955), - (128.232422, -15.36895), - (125.947266, -13.581921), - (121.816406, -18.729502), - (117.421875, -20.632784), - (113.378906, -23.402765), - (114.169922, -26.667096), - ]) - line11 = GI.LineString([ - (117.861328, -15.029686), - (122.124023, -24.886436), - (132.583008, -22.309426), - (132.890625, -7.754537), - ]) - - points = GO.intersection(line10, line11) - @test all(points[1] .≈ (119.832884, -19.58857)) - @test all(points[2] .≈ (132.808697, -11.6309378)) - - @test GO.crosses(GI.LineString([(-2, 2), (4, 2)]), line6) == true - @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true - @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 -end +# @testset "booleans" begin +# line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]]) +# line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]]) +# line3 = GI.LineString([ +# (-111.544189453125, 24.186847428521244), +# (-110.687255859375, 24.966140159912975), +# (-110.4510498046875, 24.467150664739002), +# (-109.9951171875, 25.180087808990645) +# ]) +# line4 = GI.LineString([ +# (-111.4617919921875, 24.05148034322011), +# (-110.8795166015625, 24.681961205014595), +# (-110.841064453125, 24.14174098050432), +# (-109.97863769531249, 24.617057340809524) +# ]) + +# # @test isparallel(line1, line2) == true +# # @test isparallel(line3, line4) == false + +# poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]]) +# poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]]) + +# @test isconcave(poly1) == true +# @test isconcave(poly2) == false + +# l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]]) +# l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]]) + +# @test isclockwise(l1) == true +# @test isclockwise(l2) == false \ No newline at end of file diff --git a/test/methods/equals.jl b/test/methods/equals.jl index ea1eb73fb..f24f54908 100644 --- a/test/methods/equals.jl +++ b/test/methods/equals.jl @@ -1,3 +1,5 @@ +# # Point and Geometry + @testset "Points/MultiPoints" begin p1 = LG.Point([0.0, 0.0]) p2 = LG.Point([0.0, 1.0]) diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl new file mode 100644 index 000000000..3fa793913 --- /dev/null +++ b/test/methods/geom_relations.jl @@ -0,0 +1,178 @@ +# Tests of DE-9IM Methods + +pt1 = LG.Point([0.0, 0.0]) +pt2 = LG.Point([5.0, 5.0]) +pt3 = LG.Point([1.0, 0.0]) +pt4 = LG.Point([0.5, 0.0]) +pt5 = LG.Point([0.5, 0.25]) +pt6 = LG.Point([0.6, 0.4]) +pt7 = LG.Point([0.4, 0.8]) + +l1 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]) +l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0]]) +l3 = LG.LineString([[-1.0, 0.0], [0.0, 0.0], [1.0, 0.0]]) +l4 = LG.LineString([[0.5, 0.0], [1.0, 0.0], [1.0, 0.5]]) +l5 = LG.LineString([[0.0, 0.0], [-1.0, -1.0]]) +l6 = LG.LineString([[2.0, 2.0], [0.0, 1.0]]) +l7 = LG.LineString([[0.5, 1.0], [0.5, -1.0]]) +l8 = LG.LineString([[0.0, 0.0], [0.5, 0.0], [0.5, 0.5], [1.0, -0.5]]) +l9 = LG.LineString([[0.0, 1.0], [0.0, -1.0], [1.0, 1.0]]) +l10 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]]) +l11 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0], [-1.0, 0.0]]) +l12 = LG.LineString([[0.6, 0.5], [0.6, 0.9]]) +l13 = LG.LineString([[2.0, 2.0], [3.0, 3.0]]) +l14 = LG.LineString([[0.6, 0.25], [0.6, 0.35]]) +l15 = LG.LineString([[0.0, 3.0], [4.0, 3.0]]) +l16 = LG.LineString([[0.3, -0.7], [1.0, 0.0], [3.0, 0.6]]) + +r1 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]]) +r2 = LG.LinearRing([[0.5, 0.2], [0.6, 0.4], [0.7, 0.2], [0.5, 0.2]]) +r3 = LG.LinearRing([[0.2, 0.7], [0.4, 0.9], [0.5, 0.6], [0.2, 0.7]]) +r4 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]) +r5 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [2.0, 1.0], [1.0, 2.0], [1.0, 1.0], [0.0, 0.0]]) +r6 = LG.LinearRing([[0.0, 0.0], [-1.0, 0.0], [-1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]) +r7 = LG.LinearRing([[0.5, 0.5], [1.5, 0.5], [1.5, 1.5], [0.5, 1.5], [0.5, 0.5]]) +r8 = LG.LinearRing([[0.1, 0.1], [0.2, 0.1], [0.2, 0.2], [0.1, 0.2], [0.1, 0.1]]) +r9 = LG.LinearRing([[1.0, -0.5], [0.5, -1.0], [1.5, -1.0], [1.0, -0.5]]) +r10 = LG.LinearRing([[0.5, 0.2], [0.6, 0.4], [0.3, 0.3], [0.5, 0.2]]) +r11 = LG.LinearRing([[0.55, 0.21], [0.55, 0.23], [0.65, 0.23], [0.66, 0.21], [0.55, 0.21]]) + +p1 = LG.Polygon(r4, [r2, r3]) +p2 = LG.Polygon(r4) +p3 = LG.Polygon(r2) +p4 = LG.Polygon(r6) +p5 = LG.Polygon(r9) +p6 = LG.Polygon(r11) +p7 = LG.Polygon(r7) +p8 = LG.Polygon([[[0.0, 0.0], [0.0, 3.0], [1.0, 2.0], [2.0, 3.0], [3.0, 2.0], [4.0, 3.0], [4.0, 0.0], [0.0, 0.0]]]) +p9 = LG.Polygon([[[0.0, 0.0], [0.0, 3.0], [1.0, 4.0], [2.0, 3.0], [3.0, 4.0], [4.0, 3.0], [4.0, 0.0], [0.0, 0.0]]]) +p10 = LG.Polygon([ + [[0.1, 0.5], [0.1, 0.99], [0.6, 0.99], [0.6, 0.5], [0.1, 0.5]], + [[0.15, 0.55], [0.15, 0.95], [0.55, 0.95], [0.55, 0.55], [0.15, 0.55]] +]) +p11 = LG.Polygon(r3) + +mpt1 = LG.MultiPoint([pt1, pt2]) +mpt2 = LG.MultiPoint([pt2, pt3]) +mpt3 = LG.MultiPoint([pt4, pt5]) +ml1 = LG.MultiLineString([l5, l6, l7]) +ml2 = LG.MultiLineString([l1]) +mp1 = LG.MultiPolygon([p1]) +mp2 = LG.MultiPolygon([p6, p7]) +gc1 = LG.GeometryCollection([pt1, l5, p6]) + +test_pairs = [ + # Points and geometries + (pt1, pt1, "pt1", "pt1"), # Same point + (pt1, pt2, "pt1", "pt2"), # Different point + (pt1, l1, "pt1", "l1"), # Point on line endpoint + (pt2, l1, "pt2", "l1"), # Point outside line + (pt3, l1, "pt3", "l1"), # Point on line segment + (pt4, l1, "pt4", "l1"), # Point on line vertex between segments + (l1, pt3, "l1", "pt3"), # oint on line segment (order swapped) + (pt1, r1, "pt1", "r1"), # Point on ring "endpoint" + (pt2, r1, "pt2", "r1"), # Point outside ring + (pt3, r1, "pt3", "r1"), # Point on ring segment + (pt4, r1, "pt4", "r1"), # Point on ring vertex between segments + (r1, pt3, "r1", "pt3"), # Point on ring segment (order swapped) + (p1, pt1, "p1", "pt1"), # Point on vertex of polygon + (pt2, p1, "pt2", "p1"), # Point outside of polygon's external ring + (pt4, p1, "pt4", "p1"), # Point on polygon's edge + (pt5, p1, "pt5", "p1"), # Point inside of polygon + (pt6, p1, "pt6", "p1"), # Point on hole edge + (pt7, p1, "pt7", "p1"), # Point inside of polygon hole + (p1, pt5, "p1", "pt5"), # Point inside of polygon (order swapped) + # # Lines and geometries + (l1, l1, "l1", "l1"), # Same line + (l2, l1, "l2", "l1"), # L2 is one segment of l1 + (l3, l1, "l3", "l1"), # L3 shares one segment with l1 and has one segment outside + (l4, l1, "l4", "l1"), # L4 shares half of each of l1 segments + (l5, l1, "l5", "l1"), # L5 shares one endpoint with l1 but not segments + (l6, l1, "l6", "l1"), # Lines are disjoint + (l7, l1, "l7", "l1"), # L7 crosses through one of l1's segments + (l8, l1, "l8", "l1"), # Overlaps one segment and crosses through another segment + (l9, l1, "l9", "l1"), # Two segments touch and two segments crosses + (l16, l1, "l16", "l1"), # L16 bounces off of l1's corner + (l1, r1, "l1", "r1"), # Line inside of ring + (l3, r1, "l3", "r1"), # Line covers one edge of linear ring and has segment outside + (l5, r1, "l5", "r1"), # Line and linear ring are only covered by vertex + (l6, r1, "l6", "r1"), # Line and ring are disjoint + (l7, r1, "l7", "r1"), # Line crosses through two ring edges + (l8, r1, "l8", "r1"), # Line crosses through two ring edges and touches third edge + (l10, r1, "l10", "r1"), # Line is equal to linear ring + (l11, r1, "l11", "r1"), # Line covers linear ring and then has extra segment + (l1, p1, "l1", "p1"), # Line on polygon edge + (l3, p1, "l3", "p1"), # Line on polygon edge and extending beyond polygon edge + (l5, p1, "l5", "p1"), # Line outside polygon connected by a vertex + (l7, p1, "l7", "p1"), # Line through polygon cutting to the outside + (l12, p1, "l12", "p1"), # Line inside polygon + (l13, p1, "l13", "p1"), # Line outside of polygon + (l14, p1, "l14", "p1"), # Line in polygon hole + (l15, p8, "l15", "p8"), # Line outside crown-shaped polygon but touching edges + (l15, p9, "l15", "p9"), # Line within crown-shaped polygon but touching edges + # Ring and geometries + (r1, l1, "r1", "l1"), # Line is within linear ring + (r1, l3, "r1", "l3"), # Line covers one edge of linear ring and has segment outside + (r1, l5, "r1", "l5"), # Line and linear ring are only connected at vertex + (r1, l6, "r1", "l6"), # Line and linear ring are disjoint + (r1, l7, "r1", "l7"), # Line crosses though two ring edges + (r1, l8, "r1", "l8"), # Line crosses through two ring edges and touches third edge + (r1, l10, "r1", "l10"), # Line is equal to linear ring + (r1, l11, "r1", "l11"), # Line covers linear ring and then has extra segment + (r1, r1, "r1", "r1"), # Same rings + (r2, r1, "r2", "r1"), # Disjoint ring with one "inside" of hole created + (r3, r1, "r3", "r1"), # Disjoint ring with one "outside" of hole created + (r4, r1, "r4", "r1"), # Rings share two sides and rest of sides dont touch + (r1, r5, "r1", "r5"), # Ring shares all edges with other ring, plus an extra loop + (r1, r6, "r1", "r6"), # Rings share just one vertex + (r1, r7, "r1", "r7"), # Rings cross one another + (r4, p1, "r4", "p1"), # Ring on boundary of polygon + (r1, p1, "r1", "p1"), # Ring on boundary and cutting through polygon + (r2, p1, "r2", "p1"), # Ring on hole bounday + (r6, p1, "r6", "p1"), # Ring touches polygon at one vertex + (r7, p1, "r7", "p1"), # Ring crosses through polygon + (r8, p1, "r8", "p1"), # Ring inside of polygon + (r9, p1, "r9", "p1"), # Ring outside of polygon + (r10, p1, "r10", "p1"), # Ring inside of polygon and shares hole's edge + (r11, p1, "r11", "p1"), # Ring inside of polygon hole + # Polygon and geometries + (p1, p1, "p1", "p1"), # Same polygons + (p1, p2, "p1", "p2"), # P1 and p2 are the same but p1 has holes + (p2, p1, "p2", "p1"), # P1 and p2 are the same but p1 has holes (order swapped) + (p3, p1, "p3", "p1"), # P3 is equal to one of p1's holes + (p4, p1, "p4", "p1"), # Polygon's share just one vertex + (p5, p1, "p5", "p1"), # Polygon outside of other polygon + (p6, p1, "p6", "p1"), # Polygon inside of other polygon's hole + (p7, p1, "p7", "p1"), # Polygons overlap + (p10, p1, "p10", "p1"), # Polygon's with nested holes + # Multigeometries + (mpt1, mpt1, "mpt1", "mpt1"), # Same set of points for multipoints + (mpt1, mpt2, "mpt1", "mpt2"), # Some point matches, others are different + (mpt1, mpt3, "mpt1", "mpt3"), # No shared points + (ml1, ml2, "ml1", "ml2"), # Lines in ml1 cross and touch ml2 + (mp1, mp2, "mp1", "mp2"), # Polygons in mp1 are inside hole and overlap + (gc1, ml1, "gc1", "ml1"), # Make sure collection works with multi-geom +] + +function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) + for (g1, g2, sg1, sg2) in test_pairs + if swap_points + g1, g2 = g2, g1 + sg1, sg2 = sg2, sg1 + end + go_val = GO_f(g1, g2) + lg_val = LG_f(g1, g2) + @test go_val == lg_val + go_val != lg_val && println("\n" * "↑ TEST INFO: " * sg1 * " " * f_name * " " * sg2 * "\n\n") + end +end + +@testset "Contains" begin test_geom_relation(GO.contains, LG.contains, "contains"; swap_points = true) end +@testset "Covered By" begin test_geom_relation(GO.coveredby, LG.coveredby, "coveredby") end +@testset "Covers" begin test_geom_relation(GO.covers, LG.covers, "covers"; swap_points = true) end +@testset "Crosses" begin test_geom_relation(GO.crosses, LG.crosses, "crosses") end +@testset "Disjoint" begin test_geom_relation(GO.disjoint, LG.disjoint, "disjoint")end +@testset "Intersect" begin test_geom_relation(GO.intersects, LG.intersects, "intersects") end +@testset "Overlaps" begin test_geom_relation(GO.overlaps, LG.overlaps, "overlaps") end +@testset "Touches" begin test_geom_relation(GO.touches, LG.touches, "touches") end +@testset "Within" begin test_geom_relation(GO.within, LG.within, "within") end \ No newline at end of file diff --git a/test/methods/intersects.jl b/test/methods/intersects.jl deleted file mode 100644 index 4251d45a8..000000000 --- a/test/methods/intersects.jl +++ /dev/null @@ -1,145 +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) - @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) - @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) - @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) - @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) - 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) - 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) - 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) - @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) - @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) - 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) - 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) - @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) - @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) - 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) - @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) - @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/methods/overlaps.jl b/test/methods/overlaps.jl deleted file mode 100644 index 0c2dab96d..000000000 --- a/test/methods/overlaps.jl +++ /dev/null @@ -1,104 +0,0 @@ -@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) - @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 - @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) - - 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 779cfa15e..84a88663d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,13 +16,12 @@ const GO = GeometryOps @testset "Primitives" begin include("primitives.jl") end # Methods @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 "Geometry Relations" begin include("methods/geom_relations.jl") end + @testset "Bools" begin include("methods/bools.jl") end @testset "Area" begin include("methods/area.jl") end @testset "Distance" begin include("methods/distance.jl") end - @testset "Overlaps" begin include("methods/overlaps.jl") end # Transformations @testset "Embed Extent" begin include("transformations/extent.jl") end @testset "Reproject" begin include("transformations/reproject.jl") end