From e3a94b49d412ce36cb141ebd6ec68f1e32c4a27a Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 20:33:22 -0400
Subject: [PATCH 01/14] Add an experimental adjacency graph method

---
 docs/src/experiments/adjacency_graph.jl | 110 ++++++++++++++++++++++++
 1 file changed, 110 insertions(+)
 create mode 100644 docs/src/experiments/adjacency_graph.jl

diff --git a/docs/src/experiments/adjacency_graph.jl b/docs/src/experiments/adjacency_graph.jl
new file mode 100644
index 0000000000..a9a54c2527
--- /dev/null
+++ b/docs/src/experiments/adjacency_graph.jl
@@ -0,0 +1,110 @@
+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

From 855f6256e573e31fc9f1c9669cbb7e84e1a917dc Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 21:04:36 -0400
Subject: [PATCH 02/14] Basic framework for perimeter finding

---
 src/methods/perimeter.jl | 91 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 91 insertions(+)
 create mode 100644 src/methods/perimeter.jl

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
new file mode 100644
index 0000000000..af23dc560c
--- /dev/null
+++ b/src/methods/perimeter.jl
@@ -0,0 +1,91 @@
+#=
+# 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}()
+
+function perimeter(geom, ::Type{T}; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: Number
+    _threaded = booltype(threaded)
+    return applyreduce(_perimeter, +, geom; threaded = _threaded)
+end
+
+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, ::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)
+    curr = prev
+    for point in GI.getpoint(geom)
+        ret += point_distance(alg, prev, point, T)
+        prev = point
+    end
+    return ret
+end
+
+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")
+#=
+function GO.point_distance(::GO.GeodesicDistance, p1, p2, ::Type{T}) where T <: Number
+
+end
+# point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = T(Geod.geod(p1, p2)[])
\ No newline at end of file

From 260c95fe967e519a64c0124ed86b01e93d4186ff Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 21:05:24 -0400
Subject: [PATCH 03/14] Flesh out the commented-out geodesic definition

---
 src/methods/perimeter.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index af23dc560c..e777d9e98a 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -52,6 +52,7 @@ abstract type DistanceAlgorithm end
 A linear distance algorithm that uses simple 2D Euclidean distance between points.
 """
 struct LinearDistance <: DistanceAlgorithm end
+
 """
     GeodesicDistance()
 
@@ -62,6 +63,7 @@ Requires the Proj.jl package to be loaded, uses Proj's GeographicLib.
 struct GeodesicDistance{T} <: DistanceAlgorithm 
     geodesic::T# ::Proj.geod_geodesic
 end
+
 """
     RhumbDistance()
 
@@ -85,7 +87,13 @@ end
 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")
 #=
-function GO.point_distance(::GO.GeodesicDistance, p1, p2, ::Type{T}) where T <: Number
-
+function GO.point_distance(alg::GO.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, _ = Proj.geod_inverse(alg.geodesic, lon1, lat1, lon2, lat2)
+    return T(dist)
 end
 # point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = T(Geod.geod(p1, p2)[])
\ No newline at end of file

From c87a878e9cba3dbf5c511fd9b4dbca64f5b9779c Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 21:09:08 -0400
Subject: [PATCH 04/14] fix

---
 src/methods/perimeter.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index e777d9e98a..88d41f74f3 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -93,7 +93,7 @@ function GO.point_distance(alg::GO.GeodesicDistance, p1, p2, ::Type{T}) where T
     lon2 = Base.convert(Float64, GI.x(p2))
     lat2 = Base.convert(Float64, GI.y(p2))
 
-    dist, _ = Proj.geod_inverse(alg.geodesic, lon1, lat1, lon2, lat2)
+    dist, _azi1, _azi2 = Proj.geod_inverse(alg.geodesic, lon1, lat1, lon2, lat2)
     return T(dist)
 end
 # point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = T(Geod.geod(p1, p2)[])
\ No newline at end of file

From 4ea0353b44f55cbcf24edf4d622e70d47f6e9044 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 21:27:18 -0400
Subject: [PATCH 05/14] more fix

---
 src/methods/perimeter.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index 88d41f74f3..c1749042dd 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -42,7 +42,8 @@ const _PERIMETER_TARGETS = TraitTarget{GI.AbstractCurveTrait}()
 
 function perimeter(geom, ::Type{T}; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: Number
     _threaded = booltype(threaded)
-    return applyreduce(_perimeter, +, geom; threaded = _threaded)
+    find_perimeter(geom) = _perimeter(alg, geom, T)
+    return applyreduce(find_perimeter, +, geom; threaded = _threaded)
 end
 
 abstract type DistanceAlgorithm end
@@ -76,7 +77,6 @@ _perimeter(alg::DistanceAlgorithm, geom, ::Type{T}) where T <: Number = _perimet
 function _perimeter(::GI.AbstractCurveTrait, alg::DistanceAlgorithm, geom, ::Type{T}) where T <: Number
     ret = T(0)
     prev = GI.getpoint(geom, 1)
-    curr = prev
     for point in GI.getpoint(geom)
         ret += point_distance(alg, prev, point, T)
         prev = point
@@ -96,4 +96,6 @@ function GO.point_distance(alg::GO.GeodesicDistance, p1, p2, ::Type{T}) where T
     dist, _azi1, _azi2 = Proj.geod_inverse(alg.geodesic, lon1, lat1, lon2, lat2)
     return T(dist)
 end
-# point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = T(Geod.geod(p1, p2)[])
\ No newline at end of file
+=#
+
+# point_distance(::RhumbDistance, p1, p2, ::Type{T}) where T <: Number = ...

From 11b90994526baea42c1bbbbd94d6fbb3706a5d96 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 22:04:56 -0400
Subject: [PATCH 06/14] Construct a graph and visualize it

---
 docs/Project.toml                       |  4 ++-
 docs/src/experiments/adjacency_graph.jl | 48 +++++++++++++++++++++++++
 2 files changed, 51 insertions(+), 1 deletion(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 268f762492..01d492ef92 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,15 +12,17 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233"
-ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoDatasets = "ddc7317b-88db-5cb5-a849-8449e5df04f9"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoInterfaceMakie = "0edc0954-3250-4c18-859d-ec71c1660c08"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 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"
diff --git a/docs/src/experiments/adjacency_graph.jl b/docs/src/experiments/adjacency_graph.jl
index a9a54c2527..1b60bab835 100644
--- a/docs/src/experiments/adjacency_graph.jl
+++ b/docs/src/experiments/adjacency_graph.jl
@@ -108,3 +108,51 @@ record(f, "adjacency_matrix.mp4", 1:size(us_states, 1); framerate = 1) do i
     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 = 1000), 
+    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)
+end
\ No newline at end of file

From 55f122ab4e1a367e815f1f62b1b3a4ee136b7fb4 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 22:05:03 -0400
Subject: [PATCH 07/14] Include the perimeter file

---
 src/GeometryOps.jl        |  1 +
 src/methods/perimeter.jl  | 10 +++++++---
 test/methods/perimeter.jl |  0
 3 files changed, 8 insertions(+), 3 deletions(-)
 create mode 100644 test/methods/perimeter.jl

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 1aced0e820..cee8d7998a 100644
--- a/src/GeometryOps.jl
+++ b/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")
diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index c1749042dd..01280878ee 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -40,10 +40,11 @@ GO.perimeter(mp) # should be 92
 
 const _PERIMETER_TARGETS = TraitTarget{GI.AbstractCurveTrait}()
 
-function perimeter(geom, ::Type{T}; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: Number
-    _threaded = booltype(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, +, geom; threaded = _threaded)
+    return applyreduce(find_perimeter, +, _PERIMETER_TARGETS, geom; threaded = _threaded, init = zero(T))
 end
 
 abstract type DistanceAlgorithm end
@@ -54,6 +55,9 @@ A linear distance algorithm that uses simple 2D Euclidean distance between point
 """
 struct LinearDistance <: DistanceAlgorithm end
 
+perimeter(geom, _T::Type{T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: Number = perimeter(LinearDistance(), geom, T; threaded = threaded)
+
+
 """
     GeodesicDistance()
 
diff --git a/test/methods/perimeter.jl b/test/methods/perimeter.jl
new file mode 100644
index 0000000000..e69de29bb2

From d0828a963d9fe54f9cf76340d44b20c17b1ecd9b Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 22:11:33 -0400
Subject: [PATCH 08/14] Fix loading

---
 src/methods/perimeter.jl | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index 01280878ee..b97c10d447 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -41,12 +41,6 @@ GO.perimeter(mp) # should be 92
 const _PERIMETER_TARGETS = TraitTarget{GI.AbstractCurveTrait}()
 
 
-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
-
 abstract type DistanceAlgorithm end
 """
     LinearDistance()
@@ -55,9 +49,6 @@ A linear distance algorithm that uses simple 2D Euclidean distance between point
 """
 struct LinearDistance <: DistanceAlgorithm end
 
-perimeter(geom, _T::Type{T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: Number = perimeter(LinearDistance(), geom, T; threaded = threaded)
-
-
 """
     GeodesicDistance()
 
@@ -76,6 +67,16 @@ A rhumb distance algorithm that uses the rhumb distance between points.
 """
 struct RhumbDistance <: DistanceAlgorithm end
 
+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

From 1e9b97f884a9fde28592514e016355fbc9f255d2 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 22:38:27 -0400
Subject: [PATCH 09/14] Skeleton for Albers USA projection

---
 docs/src/experiments/adjacency_graph.jl | 26 +++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/docs/src/experiments/adjacency_graph.jl b/docs/src/experiments/adjacency_graph.jl
index 1b60bab835..bbd79e9e1b 100644
--- a/docs/src/experiments/adjacency_graph.jl
+++ b/docs/src/experiments/adjacency_graph.jl
@@ -155,4 +155,30 @@ record(f, "graph_tightening.mp4", exp10.(LinRange(log10(10), log10(10_000), 100)
     niters = round(Int, i)
     a.title = "$niters iterations"
     p.layout = GraphMakie.Spring(; iterations = niters)
+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
\ No newline at end of file

From e8659cdbab1a9abeeae95436c05b2e30206cd453 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 12 Jun 2024 22:46:45 -0400
Subject: [PATCH 10/14] Add some docs

---
 src/methods/perimeter.jl | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index b97c10d447..96e2d15f11 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -40,8 +40,17 @@ 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()
 
@@ -67,7 +76,14 @@ A rhumb distance algorithm that uses the rhumb distance between points.
 """
 struct RhumbDistance <: DistanceAlgorithm end
 
-perimeter(geom, T::Type{_T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where _T <: Number = perimeter(LinearDistance(), geom, _T; threaded = threaded)
+"""
+    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(alg::DistanceAlgorithm, 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)

From 5675c91b6fff4e396ebae8b134582d6a091debbc Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Thu, 13 Jun 2024 08:35:50 -0400
Subject: [PATCH 11/14] Fully implement the geodesic method

---
 ext/GeometryOpsProjExt/GeometryOpsProjExt.jl |  1 +
 ext/GeometryOpsProjExt/perimeter.jl          | 14 ++++++++++++
 src/GeometryOps.jl                           |  2 ++
 src/methods/perimeter.jl                     | 23 ++++++++++----------
 4 files changed, 28 insertions(+), 12 deletions(-)
 create mode 100644 ext/GeometryOpsProjExt/perimeter.jl

diff --git a/ext/GeometryOpsProjExt/GeometryOpsProjExt.jl b/ext/GeometryOpsProjExt/GeometryOpsProjExt.jl
index aa89b9f967..ec0e0f4727 100644
--- a/ext/GeometryOpsProjExt/GeometryOpsProjExt.jl
+++ b/ext/GeometryOpsProjExt/GeometryOpsProjExt.jl
@@ -2,6 +2,7 @@ module GeometryOpsProjExt
 
 using GeometryOps, Proj
 
+include("perimeter.jl")
 include("reproject.jl")
 include("segmentize.jl")
 
diff --git a/ext/GeometryOpsProjExt/perimeter.jl b/ext/GeometryOpsProjExt/perimeter.jl
new file mode 100644
index 0000000000..903ab3636c
--- /dev/null
+++ b/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, lon1, lat1, lon2, lat2)
+    return T(dist)
+end
\ No newline at end of file
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index cee8d7998a..2b2898d832 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -74,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
diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index 96e2d15f11..074a1b7dc4 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -83,7 +83,7 @@ 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(alg::DistanceAlgorithm, geom, T::Type{_T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where _T <: Number = perimeter(LinearDistance(), geom, _T; threaded = threaded)
+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)
@@ -107,16 +107,15 @@ end
 
 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")
-#=
-function GO.point_distance(alg::GO.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, lon1, lat1, lon2, lat2)
-    return T(dist)
-end
-=#
 
 # 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
\ No newline at end of file

From f8ce98ba3ccfeb6810fccbe7162c79fc1c3f1829 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Thu, 13 Jun 2024 11:10:31 -0400
Subject: [PATCH 12/14] Allow initial positions to be graph centroids

---
 docs/src/experiments/adjacency_graph.jl | 4 ++--
 test/runtests.jl                        | 1 +
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/adjacency_graph.jl b/docs/src/experiments/adjacency_graph.jl
index bbd79e9e1b..3847b84e48 100644
--- a/docs/src/experiments/adjacency_graph.jl
+++ b/docs/src/experiments/adjacency_graph.jl
@@ -146,7 +146,7 @@ end
 
 f, a, p = GraphMakie.graphplot(
     g; 
-    layout = GraphMakie.Spring(; iterations = 1000), 
+    layout = GraphMakie.Spring(; iterations = 10000, initialpos = GO.tuples(LG.centroid.(us_states.geometry))), 
     ilabels = abbrevs, 
     figure = (; size = (1000, 1000))
 )
@@ -154,7 +154,7 @@ f, a, p = GraphMakie.graphplot(
 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)
+    p.layout = GraphMakie.Spring(; iterations = niters, initialpos = GO.tuples(LG.centroid.(us_states.geometry)))
 end
 
 
diff --git a/test/runtests.jl b/test/runtests.jl
index f00beb439a..8ad9f5863f 100644
--- a/test/runtests.jl
+++ b/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

From 7f88e7ffd7911f293e7222674ac0213f3cfdabc9 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Thu, 13 Jun 2024 12:11:38 -0400
Subject: [PATCH 13/14] Fix the order in perimeter

---
 ext/GeometryOpsProjExt/perimeter.jl | 2 +-
 src/methods/perimeter.jl            | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/ext/GeometryOpsProjExt/perimeter.jl b/ext/GeometryOpsProjExt/perimeter.jl
index 903ab3636c..e5c9501543 100644
--- a/ext/GeometryOpsProjExt/perimeter.jl
+++ b/ext/GeometryOpsProjExt/perimeter.jl
@@ -9,6 +9,6 @@ function GeometryOps.point_distance(alg::GeometryOps.GeodesicDistance, p1, p2, :
     lon2 = Base.convert(Float64, GI.x(p2))
     lat2 = Base.convert(Float64, GI.y(p2))
 
-    dist, _azi1, _azi2 = Proj.geod_inverse(alg.geodesic, lon1, lat1, lon2, lat2)
+    dist, _azi1, _azi2 = Proj.geod_inverse(alg.geodesic, lat1, lon1, lat2, lon2)
     return T(dist)
 end
\ No newline at end of file
diff --git a/src/methods/perimeter.jl b/src/methods/perimeter.jl
index 074a1b7dc4..91e3281498 100644
--- a/src/methods/perimeter.jl
+++ b/src/methods/perimeter.jl
@@ -105,6 +105,7 @@ function _perimeter(::GI.AbstractCurveTrait, alg::DistanceAlgorithm, geom, ::Typ
     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")
 

From 55f4a12b85d8eead1bc2764dfd0bc1ea44a90ed6 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Thu, 13 Jun 2024 12:11:59 -0400
Subject: [PATCH 14/14] Add an actual test

---
 test/methods/perimeter.jl | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/test/methods/perimeter.jl b/test/methods/perimeter.jl
index e69de29bb2..dd4f023808 100644
--- a/test/methods/perimeter.jl
+++ b/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