overload hash and equality (#147)

* overload hash and equality * Update src/geos_types.jl * bump version * wip * fix hash * implement == via equals * implement coordinatewise equality * overload Base.isapprox * more isapprox tests
JuliaGeo · Feb 15, 2023 · ffe4d6e · ffe4d6e
1 parent f56f3f9
commit ffe4d6e
Show file tree

Hide file tree

Showing 4 changed files with 169 additions and 1 deletion.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "LibGEOS"
 uuid = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 license = "MIT"
-version = "0.7.5"
+version = "0.7.6"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"

diff --git a/src/LibGEOS.jl b/src/LibGEOS.jl
@@ -127,6 +127,13 @@ mutable struct GEOSContext
     end
 end
 
+function Base.:(==)(c1::GEOSContext, c2::GEOSContext)
+    (c1.ptr == c2.ptr)
+end
+function Base.hash(c::GEOSContext, h::UInt)
+    hash(c.ptr, h)
+end
+
 "Get a copy of a string from GEOS, freeing the GEOS managed memory."
 function string_copy_free(s::Cstring, context::GEOSContext = get_global_context())::String
     copy = unsafe_string(s)

diff --git a/src/geos_types.jl b/src/geos_types.jl
@@ -309,6 +309,97 @@ const geomtypes = [
     GeometryCollection,
 ]
 
+Base.@kwdef struct IsApprox
+    atol::Float64=0.0
+    rtol::Float64=sqrt(eps(Float64))
+end
+
+function Base.:(==)(geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
+    compare(==, geo1, geo2)
+end
+function Base.isequal(geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
+    compare(isequal, geo1, geo2)
+end
+function Base.isapprox(geo1::AbstractGeometry, geo2::AbstractGeometry; kw...)::Bool
+    compare(IsApprox(;kw...), geo1, geo2)
+end
+function (cmp::IsApprox)(geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
+    compare(cmp, geo1, geo2)
+end
+function compare(≅, geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
+    (typeof(geo1) === typeof(geo2)) || return false
+    ng1 = ngeom(geo1)
+    ng2 = ngeom(geo2)
+    ng1 == ng2 || return false
+    for i in 1:ng1
+        (getgeom(geo1, i) ≅ getgeom(geo2, i)) || return false
+    end
+    return true
+end
+function compare(≅, pt1::Point, pt2::Point)::Bool
+    is3d = GeoInterface.is3d(pt1)
+    is3d === GeoInterface.is3d(pt2) || return false
+    GeoInterface.getcoord(pt1,1) ≅ GeoInterface.getcoord(pt2,1) || return false
+    GeoInterface.getcoord(pt1,2) ≅ GeoInterface.getcoord(pt2,2) || return false
+    if is3d
+        GeoInterface.getcoord(pt1,3) ≅ GeoInterface.getcoord(pt2,3) || return false
+    end
+    return true
+end
+
+function _norm(x,y,z)
+    return sqrt(x^2 + y^2 + z^2)
+end
+
+function compare(cmp::IsApprox, pt1::Point, pt2::Point)::Bool
+    is3d = GeoInterface.is3d(pt1)
+    is3d === GeoInterface.is3d(pt2) || return false
+    x1 = GeoInterface.getcoord(pt1,1)
+    x2 = GeoInterface.getcoord(pt2,1)
+    y1 = GeoInterface.getcoord(pt1,2)
+    y2 = GeoInterface.getcoord(pt2,2)
+    if is3d
+        z1 = GeoInterface.getcoord(pt1,3) 
+        z2 = GeoInterface.getcoord(pt2,3)
+    else
+        z1 = 0.0
+        z2 = 0.0
+    end
+    lhs = _norm(x1 - x2, y1 - y2, z1 - z2)
+    rhs = cmp.atol + cmp.rtol * max(_norm(x1,y1,z2), _norm(x2,y2,z2))
+    return lhs <= rhs
+end
+
+typesalt(::Type{GeometryCollection} ) = 0xd1fd7c6403c36e5b
+typesalt(::Type{PreparedGeometry}   ) = 0xbc1a26fe2f5b7537
+typesalt(::Type{LineString}         ) = 0x712352fe219fca15
+typesalt(::Type{LinearRing}         ) = 0xac7644fd36955ef1
+typesalt(::Type{MultiLineString}    ) = 0x85aff0a53a2f2a32
+typesalt(::Type{MultiPoint}         ) = 0x6213e67dbfd3b570
+typesalt(::Type{MultiPolygon}       ) = 0xff2f957b4cdb5832
+typesalt(::Type{Point}              ) = 0x4b5c101d3843160e
+typesalt(::Type{Polygon}            ) = 0xa5c895d62ef56723
+
+function Base.hash(geo::AbstractGeometry, h::UInt)::UInt
+    salt = typesalt(typeof(geo))
+    h = hash(salt, h)
+    for i in 1:ngeom(geo)
+        g = getgeom(geo,i)
+        h = hash(g, h)
+    end
+    return h
+end
+
+function Base.hash(pt::Point, h::UInt)::UInt
+    h = hash(getcoord(pt,1), h)
+    h = hash(getcoord(pt,2), h)
+    if GeoInterface.is3d(pt)
+        h = hash(getcoord(pt,3), h)
+    end
+    h = hash(getcoord(pt,2), h)
+    return h
+end
+
 # teach ccall how to get the pointer to pass to libgeos
 # this way the Julia compiler will track the lifetimes for us
 Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractGeometry) = x.ptr

diff --git a/test/test_misc.jl b/test/test_misc.jl
@@ -1,5 +1,6 @@
 using Test
 using LibGEOS
+import GeoInterface
 
 @testset "allow_global_context!" begin
     Point(1,2,3)
@@ -45,6 +46,75 @@ end
     @test LibGEOS.intersects(p2, q1, ctx3)
 end
 
+@testset "hash eq" begin
+    ctx1 = LibGEOS.GEOSContext()
+    ctx2 = LibGEOS.GEOSContext()
+    @test ctx1 != ctx2
+    @test ctx1 == ctx1
+    @test hash(ctx1) != hash(ctx2)
+    pt1 = readgeom("POINT(0.12345 2.000 0.1)")
+    pt2 = readgeom("POINT(0.12345 2.000 0.2)")
+    @test pt1 != pt2
+    @test hash(pt1) != hash(pt2)
+    @test !isequal(pt1, pt2)
+    @test !isapprox(pt1, pt2)
+    @test !isapprox(pt1, pt2, atol=0, rtol=0)
+    @test isapprox(pt1, pt2, atol=0.2)
+    @test isapprox(pt1, pt2, rtol=0.1)
+    @test readgeom("LINESTRING (130 240, 650 240)") != readgeom("LINESTRING (130 240, -650 240)")
+    @test !(readgeom("LINESTRING (130 240, 650 240)") ≈ readgeom("LINESTRING (130 240, -650 240)"))
+    @test readgeom("LINESTRING (130 240, 650 240)") ≈ readgeom("LINESTRING (130 240, -650 240)") atol=1300
+    @test readgeom("LINESTRING (130 240, 650 240)") ≈ readgeom("LINESTRING (130 240, 650 240.00000001)")
+
+    pt = readgeom("POINT(0 NaN)")
+    @test isequal(pt, pt)
+    @test pt != pt
+    @test !(isapprox(pt, pt))
+    @test !(isapprox(pt, pt, atol=Inf))
+
+    geo = readgeom("POLYGON((1 1,1 2,2 2,2 NaN,1 1))")
+    @test isequal(geo,geo)
+    @test geo != geo
+    @test !(isapprox(geo, geo))
+    @test !(isapprox(geo, geo, atol=Inf))
+
+    geos = [
+        readgeom("POINT(0.12345 2.000 0.1)"),
+        readgeom("POINT(0.12345 2.000)"),
+        readgeom("POINT(0.12345 2.000 0.2)"),
+        readgeom("POLYGON EMPTY"),
+        readgeom("LINESTRING EMPTY"),
+        readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))"),
+        readgeom("POLYGON((1 1,1 2,2 2,2 1,1 1))"),
+        readgeom("POLYGON((1 1,1 2,2 2,2.1 1,1 1))"),
+        readgeom("LINESTRING (130 240, 650 240)"),
+        readgeom("MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))"),
+        readgeom("GEOMETRYCOLLECTION (LINESTRING (1 2, 2 2), LINESTRING (2 1, 1 1), POLYGON ((0.5 1, 1 2, 1 1, 0.5 1)), POLYGON ((9 2, 9.5 1, 2 1, 2 2, 9 2)))"),
+        readgeom("MULTIPOINT(0 0, 5 0, 10 0)"),
+    ]
+    for g1 in geos
+        for g2 in geos
+            if g1 === g2
+                @test g1 == g2
+                @test isequal(g1,g2)
+                @test hash(g1) == hash(g2)
+                @test isapprox(g1,g2)
+            else
+                @test g1 != g2
+                @test !isequal(g1,g2)
+                @test !isapprox(g1,g2)
+                @test hash(g1) != hash(g2)
+            end
+        end
+    end
+    for g in geos
+        @test g == LibGEOS.clone(g)
+        @test g ≈ LibGEOS.clone(g)
+        @test isequal(g, LibGEOS.clone(g))
+        @test hash(g) == hash(LibGEOS.clone(g))
+    end
+end
+
 @testset "show it like you build it" begin
     for geo in [
         readgeom("POINT(0 0)")