Update equals and overlaps

src/GeometryOps.jl

@@ -31,6 +31,7 @@ include("methods/overlaps.jl")
 include("methods/within.jl")
 include("methods/polygonize.jl")
 include("methods/barycentric.jl")
+include("methods/equals.jl")
 
 include("transformations/flip.jl")
 include("transformations/simplify.jl")

src/methods/bools.jl

@@ -11,7 +11,8 @@ export line_on_line, line_in_polygon, polygon_in_polygon
 """
     isclockwise(line::Union{LineString, Vector{Position}})::Bool
 
-Take a ring and return true or false whether or not the ring is clockwise or counter-clockwise.
+Take a ring and return true or false whether or not the ring is clockwise or
+counter-clockwise.
 
 ## Example
 
@@ -26,6 +27,7 @@ true
 ```
 """
 isclockwise(geom)::Bool = isclockwise(GI.trait(geom), geom)
+
 function isclockwise(::AbstractCurveTrait, line)::Bool
     sum = 0.0
     prev = GI.getpoint(line, 1)
@@ -88,30 +90,6 @@ function isconcave(poly)::Bool
     return false
 end
 
-equals(geo1, geo2) = _equals(trait(geo1), geo1, trait(geo2), geo2)
-
-_equals(::T, geo1, ::T, geo2) where T = error("Cant compare $T yet")
-function _equals(::T, p1, ::T, p2) where {T<:PointTrait}
-    GI.ncoord(p1) == GI.ncoord(p2) || return false
-    GI.x(p1) == GI.x(p2) || return false
-    GI.y(p1) == GI.y(p2) || return false
-    if GI.is3d(p1)
-        GI.z(p1) == GI.z(p2) || return false 
-    end
-    return true
-end
-function _equals(::T, l1, ::T, l2) where {T<:AbstractCurveTrait}
-    # Check line lengths match
-    GI.npoint(l1) == GI.npoint(l2) || return false
-
-    # Then check all points are the same
-    for (p1, p2) in zip(GI.getpoint(l1), GI.getpoint(l2))
-        equals(p1, p2) || return false
-    end
-    return true
-end
-_equals(t1, geo1, t2, geo2) = false
-
 # """
 #     isparallel(line1::LineString, line2::LineString)::Bool
 
@@ -193,6 +171,26 @@ function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool
     return false
 end
 
+function point_on_seg(point, start, stop)
+    # Parse out points
+    x, y = GI.x(point), GI.y(point)
+    x1, y1 = GI.x(start), GI.y(start)
+    x2, y2 = GI.x(stop), GI.y(stop)
+    Δxl = x2 - x1
+    Δyl = y2 - y1
+    # Determine if point is on segment
+    cross = (x - x1) * Δyl - (y - y1) * Δxl
+    if cross == 0  # point is on line extending to infinity
+        # is line between endpoints
+        if abs(Δxl) >= abs(Δyl)  # is line between endpoints
+            return Δxl > 0 ? x1 <= x <= x2 : x2 <= x <= x1
+        else
+            return Δyl > 0 ? y1 <= y <= y2 : y2 <= y <= y1
+        end
+    end
+    return false
+end
+
 function point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)::Bool
     x, y = GI.x(point), GI.y(point)
     x1, y1 = GI.x(start), GI.y(start)

src/methods/centroid.jl

@@ -216,7 +216,7 @@ function centroid_and_area(::GI.MultiPolygonTrait, geom)
     xcentroid *= area
     ycentroid *= area
     # Loop over any polygons within the multipolygon
-    for i in 2:GI.ngeom(geom) #poly in GI.getpolygon(geom)
+    for i in 2:GI.ngeom(geom)
         # Polygon centroid and area
         (xpoly, ypoly), poly_area = centroid_and_area(GI.getpolygon(geom, i))
         # Accumulate the area component into `area`

src/methods/crosses.jl

@@ -55,7 +55,7 @@ end
 
 function line_crosses_line(line1, line2)
     np2 = GI.npoint(line2)
-    if intersects(line1, line2; meets=MEETS_CLOSED)
+    if intersects(line1, line2)
         for i in 1:GI.npoint(line1) - 1
             for j in 1:GI.npoint(line2) - 1
                 exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both

src/methods/equals.jl

@@ -0,0 +1,192 @@
+# # Equals
+
+export equals
+
+#=
+## What is equals?
+
+The equals function checks if two geometries are equal. They are equal if they
+share the same set of points and edges.
+
+To provide an example, consider these two lines:
+```@example cshape
+using GeometryOps
+using GeometryOps.GeometryBasics
+using Makie
+using CairoMakie
+
+l1 = GI.LineString([(0.0, 0.0), (0.0, 10.0)])
+l2 = GI.LineString([(0.0, -10.0), (0.0, 3.0)])
+f, a, p = lines(GI.getpoint(l1), color = :blue)
+scatter!(GI.getpoint(l1), color = :blue)
+lines!(GI.getpoint(l2), color = :orange)
+scatter!(GI.getpoint(l2), color = :orange)
+```
+We can see that the two lines do not share a commen set of points and edges in
+the plot, so they are not equal:
+```@example cshape
+equals(l1, l2)  # returns false
+```
+
+## Implementation
+
+This is the GeoInterface-compatible implementation.
+
+First, we implement a wrapper method that dispatches to the correct
+implementation based on the geometry trait. This is also used in the
+implementation, since it's a lot less work! 
+
+Note that while we need the same set of points and edges, they don't need to be
+provided in the same order for polygons. For for example, we need the same set
+points for two multipoints to be equal, but they don't have to be saved in the
+same order. This requires checking every point against every other point in the
+two geometries we are comparing.  
+=#
+
+"""
+    equals(geom1, geom2)::Bool
+
+Compare two Geometries return true if they are the same geometry.
+
+## Examples
+```jldoctest
+import GeometryOps as GO, GeoInterface as GI
+poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
+poly2 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
+
+GO.equals(poly1, poly2)
+# output
+true
+```
+"""
+equals(geom_a, geom_b) = equals(
+    GI.trait(geom_a), geom_a,
+    GI.trait(geom_b), geom_b,
+)
+
+"""
+    equals(::T, geom_a, ::T, geom_b)::Bool
+
+Two geometries of the same type, which don't have a equals function to dispatch
+off of should throw an error.
+"""
+equals(::T, geom_a, ::T, geom_b) where T = error("Cant compare $T yet")
+
+"""
+    equals(trait_a, geom_a, trait_b, geom_b)
+
+Two geometries which are not of the same type cannot be equal so they always
+return false.
+"""
+equals(trait_a, geom_a, trait_b, geom_b) = false
+
+"""
+    equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)::Bool
+
+Two points are the same if they have the same x and y (and z if 3D) coordinates.
+"""
+function equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)
+    GI.ncoord(p1) == GI.ncoord(p2) || return false
+    GI.x(p1) == GI.x(p2) || return false
+    GI.y(p1) == GI.y(p2) || return false
+    if GI.is3d(p1)
+        GI.z(p1) == GI.z(p2) || return false 
+    end
+    return true
+end
+
+"""
+    equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)::Bool
+
+Two multipoints are equal if they share the same set of points.
+"""
+function equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)
+    GI.npoint(mp1) == GI.npoint(mp2) || return false
+    for p1 in GI.getpoint(mp1)
+        has_match = false  # if point has a matching point in other multipoint
+        for p2 in GI.getpoint(mp2)
+            if equals(p1, p2)
+                has_match = true
+                break
+            end
+        end
+        has_match || return false  # if no matching point, can't be equal
+    end
+    return true  # all points had a match
+end
+
+"""
+    equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait} ::Bool
+
+Two curves are equal if they share the same set of points going around the
+curve. 
+"""
+function equals(::T, l1, ::T, l2) where {T<:GI.AbstractCurveTrait}
+    # Check line lengths match
+    n1 = GI.npoint(l1)
+    n2 = GI.npoint(l2)
+    # TODO: do we need to account for repeated last point??
+    n1 == n2 || return false
+
+    # Find first matching point if it exists
+    p1 = GI.getpoint(l1, 1)
+    offset = findfirst(p2 -> equals(p1, p2), GI.getpoint(l2))
+    isnothing(offset) && return false
+    offset -= 1
+
+    # Then check all points are the same wrapping around line
+    for i in 1:n1
+        pi = GI.getpoint(l1, i)
+        j = i + offset
+        j = j <= n1 ? j : (j - n1)
+        pj = GI.getpoint(l2, j)
+        equals(pi, pj) || return false
+    end
+    return true
+end
+
+"""
+    equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool
+
+Two polygons are equal if they share the same exterior edge and holes.
+"""
+function equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)
+    # Check if exterior is equal
+    equals(GI.getexterior(geom_a), GI.getexterior(geom_b)) || return false
+    # Check if number of holes are equal
+    GI.nhole(geom_a) == GI.nhole(geom_b) || return false
+    # Check if holes are equal
+    for ihole in GI.gethole(geom_a)
+        has_match = false
+        for jhole in GI.gethole(geom_b)
+            if equals(ihole, jhole)
+                has_match = true
+                break
+            end
+        end
+        has_match || return false
+    end
+    return true
+end
+
+"""
+    equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool
+
+Two multipolygons are equal if they share the same set of polygons.
+"""
+function equals(::GI.MultiPolygonTrait, geom_a, ::GI.MultiPolygonTrait, geom_b)
+    # Check if same number of polygons
+    GI.npolygon(geom_a) == GI.npolygon(geom_b) || return false
+    # Check if each polygon has a matching polygon
+    for poly_a in GI.getpolygon(geom_a)
+        has_match = false
+        for poly_b in GI.getpolygon(geom_b)
+            if equals(poly_a, poly_b)
+                has_match = true
+                break
+            end
+        end
+        has_match || return false
+    end
+    return true
+end