[WIP] Experimental adjacency graph method + perimeter

docs/Project.toml

@@ -27,6 +27,8 @@ GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
+GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"

docs/src/experiments/adjacency_graph.jl

@@ -0,0 +1,184 @@
+import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+using SortTileRecursiveTree, Tables, DataFrames
+using NaturalEarth # data source
+
+using GeoMakie # enable plotting
+# We have to monkeypatch GeoMakie.to_multipoly for geometry collections:
+function GeoMakie.to_multipoly(::GI.GeometryCollectionTrait, geom)
+    geoms = collect(GI.getgeom(geom))
+    poly_and_multipoly_s = filter(x -> GI.trait(x) isa GI.PolygonTrait || GI.trait(x) isa GI.MultiPolygonTrait, geoms)
+    if isempty(poly_and_multipoly_s)
+        return GeometryBasics.MultiPolygon([GeometryBasics.Polygon(Point{2 + GI.hasz(geom) + GI.hasm(geom), Float64}[])])
+    else
+        final_multipoly = reduce((x, y) -> GO.union(x, y; target = GI.MultiPolygonTrait()), poly_and_multipoly_s)
+        return GeoMakie.to_multipoly(final_multipoly)
+    end
+end
+
+
+function _geom_vector(object)
+    if GI.trait(object) isa GI.FeatureCollectionTrait
+        return GI.geometry.(GI.getfeature(object))
+    elseif Tables.istable(object)
+        return Tables.getcolumn(object, first(GI.geometrycolumns(object)))
+    elseif object isa AbstractVector
+        return object
+    else
+        error("Given object was neither a FeatureCollection, Table, nor AbstractVector.  Could not extract a vector of geometries.  Type was $(typeof(object))")
+    end
+end
+
+function adjacency_matrix(weight_f, predicate_f, source, target; self_intersection = false, use_strtree = false)
+    @info "Config: use_strtree = $use_strtree, self_intersection = $self_intersection"
+
+    source_geoms = _geom_vector(source)
+    target_geoms = _geom_vector(target)
+
+    local target_rtree
+    if use_strtree
+        # Create an STRTree on the target geometries
+        target_rtree = SortTileRecursiveTree.STRtree(target_geoms)
+    end
+
+    # create an adjacency matrix
+    adjacency_matrix = zeros(length(source_geoms), length(target_geoms))
+
+    if use_strtree
+        # loop over the source tiles and target tiles and fill in the adjacency matrix
+        # only call weight_f on those geometries which pass:
+        # (a) the STRTree test and 
+        # (b) the predicate_f test
+        # WARNING: sometimes, STRTree may not provide correct results.
+        # Maybe use LibSpatialIndex instead?
+        for (i, source_geom) in enumerate(source_geoms)
+            targets_in_neighbourhood = SortTileRecursiveTree.query(target_rtree, source_geom)
+            for target_index in targets_in_neighbourhood
+                if predicate_f(source_geom, target_geoms[target_index])
+                    adjacency_matrix[i, target_index] = weight_f(source_geom, target_geoms[target_index])
+                end
+            end
+        end
+    else
+        # loop over the source tiles and target tiles and fill in the adjacency matrix
+        # only call weight_f on those geometries which pass:
+        # (a) the predicate_f test
+        # TODO: potential optimizations:
+        # - if weight_f is symmetric, then skip the computation if that index of the matrix
+        #   is already nonzero, i.e., full
+        for (i, source_geom) in enumerate(source_geoms)
+            for (j, target_geom) in enumerate(target_geoms)
+                if !self_intersection && i == j # self intersection
+                    continue
+                end
+                if predicate_f(source_geom, target_geom)
+                    adjacency_matrix[i, j] = weight_f(source_geom, target_geom)
+                end
+            end
+        end
+    end
+
+    if self_intersection
+        # fill in the identity line
+        for i in 1:length(source_geoms)
+            adjacency_matrix[i, i] = 1.0
+        end
+    end
+
+    return adjacency_matrix
+end
+
+# basic test using NaturalEarth US states
+
+# get the US states as a GeoInterface FeatureCollection
+us_states = DataFrame(naturalearth("admin_1_states_provinces", 110)) # 110m is only US states but 50m is all states, be careful then about filtering.
+# We also have to make all geometries valid so that they don't cause problems for the predicates!
+us_states.geometry = LG.makeValid.(GI.convert.((LG,), us_states.geometry))
+
+@time adjmat = adjacency_matrix(LG.touches, us_states, us_states) do source_geom, target_geom
+    return 1.0
+    # this could be some measure of arclength, for example.
+end
+
+f, a, p = poly(us_states.geometry |> GeoMakie.to_multipoly; color = fill(:black, size(us_states, 1)))
+
+record(f, "adjacency_matrix.mp4", 1:size(us_states, 1); framerate = 1) do i
+    a.title = "$(us_states[i, :name])"
+    colors = fill(:gray, size(us_states, 1))
+    colors[(!iszero).(view(adjmat, :, i))] .= :yellow
+    colors[i] = :red
+    p.color = colors
+end
+
+using GraphMakie
+using Graphs
+
+g = SimpleDiGraph(adjmat)
+
+abbrevs = getindex.(us_states.iso_3166_2, (4:5,))
+
+GraphMakie.graphplot(
+    g; 
+    ilabels = abbrevs, 
+    figure = (; size = (1000, 1000))
+)
+
+GraphMakie.graphplot(
+    g; 
+    layout = GraphMakie.Spring(; iterations = 100, C = 3), 
+    ilabels = abbrevs, 
+    figure = (; size = (1000, 1000))
+)
+
+GraphMakie.graphplot(
+    g; 
+    layout = GraphMakie.Spring(; iterations = 100, initialpos = GO.tuples(LG.centroid.(us_states.geometry))), 
+    ilabels = abbrevs, 
+    figure = (; size = (1000, 1000))
+)
+
+
+# Now try actually weighting by intersection distance
+
+@time adjmat = adjacency_matrix(LG.touches, us_states, us_states) do source_geom, target_geom
+    return GO.perimeter(LG.intersection(source_geom, target_geom))
+    # this could be some measure of arclength, for example.
+end
+
+f, a, p = GraphMakie.graphplot(
+    g; 
+    layout = GraphMakie.Spring(; iterations = 10000, initialpos = GO.tuples(LG.centroid.(us_states.geometry))), 
+    ilabels = abbrevs, 
+    figure = (; size = (1000, 1000))
+)
+
+record(f, "graph_tightening.mp4", exp10.(LinRange(log10(10), log10(10_000), 100)); framerate = 24) do i
+    niters = round(Int, i)
+    a.title = "$niters iterations"
+    p.layout = GraphMakie.Spring(; iterations = niters, initialpos = GO.tuples(LG.centroid.(us_states.geometry)))
+end
+
+
+
+const _ALASKA_EXTENT = GI.Extent(X = (-171.79110717773438, -129.97999572753906), Y = (54.4041748046875, 71.3577651977539))
+const _HAWAII_EXTENT = Extent(X = (-159.80050659179688, -154.80740356445312), Y = (18.916189193725586, 22.23617935180664))
+function albers_usa_projection(lon, lat)
+    if lon in (..)(_ALASKA_EXTENT.X...) && lat in (..)(_ALASKA_EXTENT.Y...)
+        return _alaska_projection(lon, lat)
+    elseif lon in (..)(_HAWAII_EXTENT.X...) && lat in (..)(_HAWAII_EXTENT.Y...)
+        return _hawaii_projection(lon, lat)
+    else
+        return _albers_projection(lon, lat)
+    end
+end
+
+function _alaska_projection(lon, lat)
+    return (lon, lat)
+end
+
+function _hawaii_projection(lon, lat)
+    return (lon, lat)
+end
+
+function _albers_projection(lon, lat)
+    return (lon, lat)
+end

ext/GeometryOpsProjExt/GeometryOpsProjExt.jl

@@ -2,6 +2,7 @@ module GeometryOpsProjExt
 
 using GeometryOps, Proj
 
+include("perimeter.jl")
 include("reproject.jl")
 include("segmentize.jl")
 

ext/GeometryOpsProjExt/perimeter.jl

@@ -0,0 +1,14 @@
+
+function GeometryOps.GeodesicDistance(; equatorial_radius::Real=6378137, flattening::Real=1/298.257223563, geodesic::Proj.geod_geodesic = Proj.geod_geodesic(equatorial_radius, flattening))
+    GeometryOps.GeodesicDistance{Proj.geod_geodesic}(geodesic)
+end
+
+function GeometryOps.point_distance(alg::GeometryOps.GeodesicDistance, p1, p2, ::Type{T}) where T <: Number
+    lon1 = Base.convert(Float64, GI.x(p1))
+    lat1 = Base.convert(Float64, GI.y(p1))
+    lon2 = Base.convert(Float64, GI.x(p2))
+    lat2 = Base.convert(Float64, GI.y(p2))
+
+    dist, _azi1, _azi2 = Proj.geod_inverse(alg.geodesic, lat1, lon1, lat2, lon2)
+    return T(dist)
+end

src/GeometryOps.jl

@@ -48,6 +48,7 @@ include("methods/geom_relations/overlaps.jl")
 include("methods/geom_relations/touches.jl")
 include("methods/geom_relations/within.jl")
 include("methods/orientation.jl")
+include("methods/perimeter.jl")
 include("methods/polygonize.jl")
 
 include("transformations/extent.jl")
@@ -73,7 +74,9 @@ function __init__()
     # Handle all available errors!
     Base.Experimental.register_error_hint(_reproject_error_hinter, MethodError)
     Base.Experimental.register_error_hint(_geodesic_segments_error_hinter, MethodError)
+    Base.Experimental.register_error_hint(_geodesic_distance_error_hinter, MethodError)
     Base.Experimental.register_error_hint(_buffer_error_hinter, MethodError)
+
 end
 
 end

src/methods/perimeter.jl

@@ -0,0 +1,122 @@
+#=
+# Perimeter
+
+Computes the perimeter of a geometry.
+
+Perimeter is not defined for points and multipoints.
+
+For linestrings and linearrings, it is the sum of the lengths of the segments.
+For polygons, it is the sum of the lengths of the exterior and holes of the polygon.
+For multipolygons, it is the sum of the perimeters of the polygons in the multipolygon.
+For geometry collections, it is the sum of the perimeters of the geometries in the collection. 
+The geometry collections cannot have point or multipoint geometries.
+
+TODOs:
+- Add support for point geometries
+- Add support for a "true 3D" perimeter
+
+```@example perimeter
+import GeoInterface as GI, GeometryOps as GO
+
+outer = GI.LinearRing([(0,0),(10,0),(10,10),(0,10),(0,0)])
+hole1 = GI.LinearRing([(1,1),(1,2),(2,2),(2,1),(1,1)])
+hole2 = GI.LinearRing([(5,5),(5,6),(6,6),(6,5),(5,5)])
+
+p = GI.Polygon([outer, hole1, hole2])
+mp = GI.MultiPolygon([
+    p, 
+    GO.transform(x -> x .+ 12, GI.Polygon([outer, hole1]))
+])
+
+(p, mp)
+```
+```@example perimeter
+GO.perimeter(p) # should be 48
+```
+```@example perimeter
+GO.perimeter(mp) # should be 92
+```
+=#
+
+const _PERIMETER_TARGETS = TraitTarget{GI.AbstractCurveTrait}()
+
+"""
+    abstract type DistanceAlgorithm
+
+An abstract supertype for a distance algorithm that computes the distance between two points.
+
+Currently used in [`GO.perimeter`](@ref GeometryOps.perimeter), but should be used in GO.distance and other functions as well...
+
+See also: [`LinearDistance`](@ref), [`GeodesicDistance`](@ref), [`RhumbDistance`](@ref).
+"""
+abstract type DistanceAlgorithm end
+
+"""
+    LinearDistance()
+
+A linear distance algorithm that uses simple 2D Euclidean distance between points.
+"""
+struct LinearDistance <: DistanceAlgorithm end
+
+"""
+    GeodesicDistance()
+
+A geodesic distance algorithm that uses the geodesic distance between points.
+
+Requires the Proj.jl package to be loaded, uses Proj's GeographicLib.
+"""
+struct GeodesicDistance{T} <: DistanceAlgorithm 
+    geodesic::T# ::Proj.geod_geodesic
+end
+
+"""
+    RhumbDistance()
+
+A rhumb distance algorithm that uses the rhumb distance between points.
+"""
+struct RhumbDistance <: DistanceAlgorithm end
+
+"""
+    perimeter([alg::DistanceAlgorithm], geom, [T::Type{<: Number} = Float64]; threaded = false)
+
+Computes the perimeter of a geometry using the specified distance algorithm.  
+
+Allowed distance algorithms are: [`LinearDistance`](@ref) (default), [`GeodesicDistance`](@ref), [`RhumbDistance`](@ref).  The latter two are not yet implemented.
+"""
+perimeter(geom, T::Type{_T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where _T <: Number = perimeter(LinearDistance(), geom, _T; threaded = threaded)
+
+function perimeter(alg::DistanceAlgorithm, geom, T::Type{_T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where _T <: Number
+    _threaded = _booltype(threaded)
+    find_perimeter(geom) = _perimeter(alg, geom, T)
+    return applyreduce(find_perimeter, +, _PERIMETER_TARGETS, geom; threaded = _threaded, init = zero(T))
+end
+
+
+
+_perimeter(alg::DistanceAlgorithm, geom, ::Type{T}) where T <: Number = _perimeter(GI.trait(geom), alg, geom, T)
+
+function _perimeter(::GI.AbstractCurveTrait, alg::DistanceAlgorithm, geom, ::Type{T}) where T <: Number
+    ret = T(0)
+    prev = GI.getpoint(geom, 1)
+    for point in GI.getpoint(geom)
+        ret += point_distance(alg, prev, point, T)
+        prev = point
+    end
+    return ret
+end
+
+point_distance(alg::DistanceAlgorithm, p1, p2) = point_distance(alg, p1, p2, Float64)
+point_distance(::LinearDistance, p1, p2, ::Type{T}) where T <: Number = T(hypot(GI.x(p2) - GI.x(p1), GI.y(p2) - GI.y(p1)))
+point_distance(alg::DistanceAlgorithm, p1, p2, ::Type{T}) where T <: Number = error("Not implemented yet for alg $alg")
+
+# point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = ...
+
+# Add an error hint for Geodesicdistance if Proj is not loaded!
+function _geodesic_distance_error_hinter(io, exc, argtypes, kwargs)
+    if isnothing(Base.get_extension(GeometryOps, :GeometryOpsProjExt)) && exc.f == GeodesicDistance
+        print(io, "\n\nThe `GeodesicDistance` method requires the Proj.jl package to be explicitly loaded.\n")
+        print(io, "You can do this by simply typing ")
+        printstyled(io, "using Proj"; color = :cyan, bold = true)
+        println(io, " in your REPL, \nor otherwise loading Proj.jl via using or import.")
+    end
+end

test/methods/perimeter.jl

@@ -0,0 +1,23 @@
+using Test
+
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+import Proj # make sure the GO extension on Proj is active
+
+# First, test against R
+
+@testset "R/SF readme" begin
+
+    outer = GI.LinearRing([(0,0),(10,0),(10,10),(0,10),(0,0)])
+    hole1 = GI.LinearRing([(1,1),(1,2),(2,2),(2,1),(1,1)])
+    hole2 = GI.LinearRing([(5,5),(5,6),(6,6),(6,5),(5,5)])
+
+    p = GI.Polygon([outer, hole1, hole2])
+    mp = GI.MultiPolygon([
+        p, 
+        GO.transform(x -> x .+ 12, GI.Polygon([outer, hole1]))
+    ])
+
+    @test GO.perimeter(p) == 48
+    @test GO.perimeter(mp) == 92
+    @test GO.perimeter(GI.GeometryCollection([p, mp])) == 48+92
+end

test/runtests.jl

@@ -25,6 +25,7 @@ const GO = GeometryOps
     @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
     @testset "Distance" begin include("methods/distance.jl") end
     @testset "Equals" begin include("methods/equals.jl") end
+    @testset "Perimeter" begin include("methods/perimeter.jl") end
     # Clipping
     @testset "Coverage" begin include("methods/clipping/coverage.jl") end
     @testset "Cut" begin include("methods/clipping/cut.jl") end