Sg/de9im functions

src/GeometryOps.jl

@@ -26,11 +26,15 @@ include("methods/bools.jl")
 include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
-include("methods/geom_relations/intersects.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/touches.jl")
 include("methods/geom_relations/within.jl")
 include("methods/polygonize.jl")
 

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.
 
@@ -128,258 +126,15 @@ end
 
 #     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)  
 
-"""
-    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

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
 ```
 

src/methods/geom_relations/contains.jl

@@ -1,25 +1,64 @@
-# # Containment
+# # 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(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
+    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).
 
-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)
+point = GI.Point((1, 2))
 
 GO.contains(line, point)
 # output
 true
 ```
 """
-contains(g1, g2)::Bool = within(g2, g1)
+contains(g1, g2) = GeometryOps.within(g2, g1)