From 33b4bfa3daa720c522335fb2c6894ca755fdc10a Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Tue, 26 Sep 2023 16:25:22 -0700
Subject: [PATCH 01/10] Rewrite centroid using GeoInterface

---
 src/methods/centroid.jl | 147 ++++++++++++++++++++--------------------
 1 file changed, 75 insertions(+), 72 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 855c46af0..873fcfc08 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -2,85 +2,88 @@
 
 export centroid
 
-# These are all GeometryBasics.jl methods so far.
-# They need to be converted to GeoInterface.
-
-# The reason that there is a `centroid_and_signed_area` function, 
-# is because in conputing the centroid, you end up computing the signed area.  
-
-# In some computational geometry applications this may be a useful
-# source of efficiency, so I added it here.
-
-# However, it's totally fine to ignore this and not have this code path.  
-# We simply need to decide on this.
-
-function centroid(ls::LineString{2, T}) where T
-    centroid = Point{2, T}(0)
-    total_area = T(0)
-    if length(ls) == 1
-        return sum(ls[1])/2
-    end
-
-    p0 = ls[1][1]
-
-    for i in 1:(length(ls)-1)
-        p1 = ls[i][2]
-        p2 = ls[i+1][2]
-        area = signed_area(p0, p1, p2)
-        centroid = centroid .+ Point{2, T}((p0[1] + p1[1] + p2[1])/3, (p0[2] + p1[2] + p2[2])/3) * area
-        total_area += area
+# ## What is the centroid?
+
+# The centroid is the geometric center of a line or area. Note that the
+# centroid does not need to be inside of a concave area or volume.
+
+# TODO: Add an example
+
+# ## 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! 
+
+"""
+
+"""
+centroid(x) = centroid(GI.trait(x), x)
+centroid_and_signed_area(x) = centroid_and_signed_area(GI.trait(x), x)
+
+function centroid(::LineStringTrait, geom)
+    FT = Float64
+    centroid = GI.Point(FT(0), FT(0))
+    length = FT(0)
+    point₁ = GI.getpoint(geom, 1)
+    for point₂ in GI.getpoint(geom)
+        length_component = sqrt(
+            (GI.x(point₂) - GI.x(point₁))^2 +
+            (GI.y(point₂) - GI.y(point₁))^2
+        )
+        # Accumulate the segment length into `length``
+        length += length_component
+        # Weighted average of segment centroids
+        centroid = centroid .+ (point₁ .+ point₂) .* length_component / 2
+        # Advance the point buffer by 1 point
+        point₁ = point₂
     end
-    return centroid ./ total_area
+    return centroid ./= length
 end
 
-# a more optimized function, so we only calculate signed area once!
-function centroid_and_signed_area(ls::LineString{2, T}) where T
-    centroid = Point{2, T}(0)
-    total_area = T(0)
-    if length(ls) == 1
-        return sum(ls[1])/2
-    end
-
-    p0 = ls[1][1]
-
-    for i in 1:(length(ls)-1)
-        p1 = ls[i][2]
-        p2 = ls[i+1][2]
-        area = signed_area(p0, p1, p2)
-        centroid = centroid .+ Point{2, T}((p0[1] + p1[1] + p2[1])/3, (p0[2] + p1[2] + p2[2])/3) * area
-        total_area += area
+function centroid_and_signed_area(::LinearRingTrait, geom)
+    FT = Float64
+    centroid = GI.Point(FT(0), FT(0))
+    area = FT(0)
+    point₁ = GI.getpoint(geom, 1)
+    for point₂ in GI.getpoint(geom)
+        area_component = GI.x(point₁) * GI.y(point₂) -
+            GI.x(point₁) * GI.y(point₂)
+        # Accumulate the segment length into `area``
+        area += area_component
+        # Weighted average of segment centroids
+        centroid = centroid .+ (point₁ .+ point₂) .* area_component
+        # Advance the point buffer by 1 point
+        point₁ = point₂
     end
-    return (centroid ./ total_area, total_area)
+    area /= 2
+    return centroid ./= 6area, area
 end
 
-function centroid(poly::GeometryBasics.Polygon{2, T}) where T
-    exterior_centroid, exterior_area = centroid_and_signed_area(poly.exterior)
-
-    total_area = exterior_area
-    interior_numerator = Point{2, T}(0)
-    for interior in poly.interiors
-        interior_centroid, interior_area = centroid_and_signed_area(interior)
-        total_area += interior_area
-        interior_numerator += interior_centroid * interior_area
+function centroid_and_signed_area(::PolygonTrait, geom)
+    # Exterior polygon centroid and area
+    centroid, area = centroid_and_signed_area(GI.getexterior(geom))
+    centroid *= abs(signed_area)
+    for hole in GI.gethole(geom)
+        interior_centroid, interior_area = centroid_and_signed_area(hole)
+        interior_area = abs(interior_area)
+        area -= interior_area
+        centroid = centroid .= interior_centroid .* interior_area
     end
-
-    return (exterior_centroid * exterior_area - interior_numerator) / total_area
-
-end
-
-function centroid(multipoly::MultiPolygon)
-    centroids = centroid.(multipoly.polygons)
-    areas = signed_area.(multipoly.polygons)
-    areas ./= sum(areas)
-
-    return sum(centroids .* areas) / sum(areas)
+    return centroid ./= abs(area), area
 end
 
-
-function centroid(rect::Rect{N, T}) where {N, T}
-    return Point{N, T}(rect.origin .- rect.widths ./ 2)
-end
-
-function centroid(sphere::HyperSphere{N, T}) where {N, T}
-    return sphere.center
+function centroid_and_signed_area(::MultiPolygonTrait, geom)
+    FT = Float64
+    centroid = GI.Point(FT(0), FT(0))
+    area = FT(0)
+    for poly in GI.getpolygon(geom)
+        poly_centroid, poly_area = centroid_and_signed_area(poly)
+        centroid = centroid .+ poly_centroid .* poly_area
+        area += poly_area
+    end
+    return centroid ./= area, area
 end
\ No newline at end of file

From 17e94ee328ef9d368e1dae8cc6d4fa775ceb22c8 Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Wed, 27 Sep 2023 12:46:45 -0700
Subject: [PATCH 02/10] Add documentation

---
 Project.toml            |   1 +
 src/methods/centroid.jl | 156 +++++++++++++++++++++++++++++++---------
 2 files changed, 122 insertions(+), 35 deletions(-)

diff --git a/Project.toml b/Project.toml
index 527dd0464..abf971b50 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <anshulsinghvi@gmail.com> and contributors"]
 version = "0.0.1-DEV"
 
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index d47dc6460..ba315e009 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -2,62 +2,132 @@
 
 export centroid
 
-# ## What is the centroid?
+#=
+## What is the centroid?
+
+The centroid is the geometric center of a line string or area(s). Note that
+the centroid does not need to be inside of a concave area.
+
+To provide an example, consider this concave polygon in the shape of a 'C':
+```@example cshape
+using GeometryOps
+using GeometryOps.GeometryBasics
+using Makie
+using CairoMakie
+
+cshape = Polygon([
+    Point(0,0), Point(0,3), Point(3,3), Point(3,2), Point(1,2),
+    Point(1,1), Point(3,1), Point(3,0), Point(0,0),
+])
+f, a, p = poly(cshape; axis = (; aspect = DataAspect()))
+```
+Let's see what the centroid looks like:
+```@example cshape
+cent = centroid(cshape)
+plot!(a, cent)
+f
+```
+The points are ordered in a clockwise fashion, which means that the signed area
+is positive.  If we reverse the order of the points, we get a negative area.
+
+## 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! 
+=#
+"""
+    centroid(geom)::GI.Point
+
+Returns the centroid of a given line segment, linear ring, polygon, or
+mutlipolygon.
+"""
+centroid(geom) = centroid(GI.trait(geom), geom)
 
-# The centroid is the geometric center of a line or area. Note that the
-# centroid does not need to be inside of a concave area or volume.
+"""
+    centroid(geom)::GI.Point
 
-# TODO: Add an example
+Returns the centroid of a given line segment.
+"""
+centroid(trait::GI.LineStringTrait, geom)  =
+    centroid_and_length(trait, geom)[1]
 
-# ## Implementation
+"""
+    centroid(geom)::GI.Point
 
-# This is the GeoInterface-compatible implementation.
+Returns the centroid of a given linear ring, polygon, or
+mutlipolygon.
+"""
+centroid(
+    trait::Union{GI.LinearRingTrait, GI.PolygonTrait, GI.MultiPolygonTrait}, 
+    geom
+) = centroid_and_signed_area(trait, geom)[1]
 
-# 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! 
+"""
+    centroid_and_length(geom)::GI.Point
 
+Returns the centroid and length of a given geom. Note this is only valid for
+line strings.
 """
+centroid_and_length(geom) = centroid_and_length(GI.trait(geom), geom)
 
 """
-centroid(x) = centroid(GI.trait(x), x)
-centroid_and_signed_area(x) = centroid_and_signed_area(GI.trait(x), x)
+    centroid_and_length(::GI.LineStringTrait, geom)::(GI.Point, ::Real)
 
-function centroid(::LineStringTrait, geom)
+Returns the centroid and length of a given line segment.
+"""
+function centroid_and_length(::GI.LineStringTrait, geom)
     FT = Float64
+    # Initialize starting values
     centroid = GI.Point(FT(0), FT(0))
     length = FT(0)
     point₁ = GI.getpoint(geom, 1)
+    # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
+        # Calculate length of line segment
         length_component = sqrt(
             (GI.x(point₂) - GI.x(point₁))^2 +
             (GI.y(point₂) - GI.y(point₁))^2
         )
-        # Accumulate the segment length into `length``
+        # Accumulate the line segment length into `length`
         length += length_component
-        # Weighted average of segment centroids
-        centroid = centroid .+ (point₁ .+ point₂) .* length_component / 2
-        # Advance the point buffer by 1 point
+        # Weighted average of line segment centroids
+        centroid = centroid .+ ((point₁ .+ point₂) .* (length_component / 2))
+        # Advance the point buffer by 1 point to move to next line segment
         point₁ = point₂
     end
-    return centroid ./= length
+    return centroid ./= length, length
 end
 
-centroid(::LienarRingTrait, geom) = 
-    centroid_and_signed_area(::LinearRingTrait, geom)[1]
+"""
+    centroid_and_length(geom)::GI.Point
 
-function centroid_and_signed_area(::LinearRingTrait, geom)
+Returns the centroid and area of a given geom. Note this is only valid for
+linear rings, polygons, and multipolygons.
+"""
+centroid_and_signed_area(geom) = centroid_and_signed_area(GI.trait(geom), geom)
+
+"""
+    centroid_and_signed_area(::GI.LinearRingTrait, geom)::(GI.Point, ::Real)
+
+Returns the centroid and signed area of a given linear ring.
+"""
+function centroid_and_signed_area(::GI.LinearRingTrait, geom)
     FT = Float64
+    # Initialize starting values
     centroid = GI.Point(FT(0), FT(0))
     area = FT(0)
     point₁ = GI.getpoint(geom, 1)
+    # Loop over line segments of linear ring
     for point₂ in GI.getpoint(geom)
         area_component = GI.x(point₁) * GI.y(point₂) -
             GI.x(point₁) * GI.y(point₂)
-        # Accumulate the segment length into `area``
+        # Accumulate the area component into `area`
         area += area_component
-        # Weighted average of segment centroids
+        # Weighted average of centroid components
         centroid = centroid .+ (point₁ .+ point₂) .* area_component
         # Advance the point buffer by 1 point
         point₁ = point₂
@@ -66,34 +136,50 @@ function centroid_and_signed_area(::LinearRingTrait, geom)
     return centroid ./= 6area, area
 end
 
+"""
+    centroid_and_signed_area(::GI.PolygonTrait, geom)::(GI.Point, ::Real)
 
-centroid(::PolygonTrait, geom) = 
-    centroid_and_signed_area(::PolygonTrait, geom)[1]
-
-function centroid_and_signed_area(::PolygonTrait, geom)
+Returns the centroid and signed area of a given polygon.
+"""
+function centroid_and_signed_area(::GI.PolygonTrait, geom)
     # Exterior polygon centroid and area
-    centroid, area = centroid_and_signed_area(GI.getexterior(geom))
+    centroid, area = centroid_and_signed_area(
+        GI.LinearRingTrait,
+        GI.getexterior(geom),
+    )
+    # Weight exterior centroid by area
     centroid *= abs(signed_area)
+    # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
+        # Hole polygon's centroid and area
         interior_centroid, interior_area = centroid_and_signed_area(hole)
         interior_area = abs(interior_area)
+        # Accumulate the area component into `area`
         area -= interior_area
-        centroid = centroid .= interior_centroid .* interior_area
+        # Weighted average of centroid components
+        centroid = centroid .- (interior_centroid .* interior_area)
     end
     return centroid ./= abs(area), area
 end
 
-centroid(::MultiPolygonTrait, geom) = 
-    centroid_and_signed_area(::MultiPolygonTrait, geom)[1]
-    
-function centroid_and_signed_area(::MultiPolygonTrait, geom)
+"""
+    centroid_and_signed_area(::GI.MultiPolygonTrait, geom)::(GI.Point, ::Real)
+
+Returns the centroid and signed area of a given multipolygon.
+"""
+function centroid_and_signed_area(::GI.MultiPolygonTrait, geom)
     FT = Float64
+    # Initialize starting values
     centroid = GI.Point(FT(0), FT(0))
     area = FT(0)
+    # Loop over any polygons within the multipolygon
     for poly in GI.getpolygon(geom)
+        # Polygon centroid and area
         poly_centroid, poly_area = centroid_and_signed_area(poly)
-        centroid = centroid .+ poly_centroid .* poly_area
+        # Accumulate the area component into `area`
         area += poly_area
+        # Weighted average of centroid components
+        centroid = centroid .+ (poly_centroid .* poly_area)
     end
     return centroid ./= area, area
-end
+end
\ No newline at end of file

From c5555c6467a15c062efd88531e97ac3369b4e706 Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Wed, 27 Sep 2023 14:44:59 -0700
Subject: [PATCH 03/10] Working on convenience functions

---
 src/methods/centroid.jl | 5 +----
 src/utils.jl            | 5 -----
 2 files changed, 1 insertion(+), 9 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index ba315e009..eeed16f06 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -143,10 +143,7 @@ Returns the centroid and signed area of a given polygon.
 """
 function centroid_and_signed_area(::GI.PolygonTrait, geom)
     # Exterior polygon centroid and area
-    centroid, area = centroid_and_signed_area(
-        GI.LinearRingTrait,
-        GI.getexterior(geom),
-    )
+    centroid, area = centroid_and_signed_area(GI.getexterior(geom))
     # Weight exterior centroid by area
     centroid *= abs(signed_area)
     # Loop over any holes within the polygon
diff --git a/src/utils.jl b/src/utils.jl
index b1fbef8ba..e7cb20d01 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -97,11 +97,6 @@ function to_extent(edges::Vector{Edge})
     Extents.Extent(X=x, Y=y)
 end
 
-function to_extent(edges::Vector{Edge})
-    x, y = extrema(first, edges)
-    Extents.Extent(X=x, Y=y)
-end
-
 function to_points(xs)
     points = Vector{TuplePoint}(undef, _npoint(x))
     _to_points!(points, x, 1)

From d86d334e8e8085d6dbadb30fef7febc533feb011 Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Wed, 27 Sep 2023 18:31:49 -0700
Subject: [PATCH 04/10] Start adding centroid tests

---
 Project.toml             |   5 +-
 src/methods/centroid.jl  | 108 ++++++++++++++++++++++-----------------
 test/methods/centroid.jl |  53 +++++++++++++++++++
 test/runtests.jl         |   7 ++-
 4 files changed, 123 insertions(+), 50 deletions(-)
 create mode 100644 test/methods/centroid.jl

diff --git a/Project.toml b/Project.toml
index abf971b50..35d7d7e1c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,10 +4,10 @@ authors = ["Anshul Singhvi <anshulsinghvi@gmail.com> and contributors"]
 version = "0.0.1-DEV"
 
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -21,7 +21,8 @@ julia = "1.3"
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "GeoFormatTypes", "GeoJSON", "Test"]
+test = ["ArchGDAL", "GeoFormatTypes", "GeoJSON", "LibGEOS", "Test"]
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index eeed16f06..134185ed1 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -21,10 +21,10 @@ cshape = Polygon([
 ])
 f, a, p = poly(cshape; axis = (; aspect = DataAspect()))
 ```
-Let's see what the centroid looks like:
+Let's see what the centroid looks like (plotted in red):
 ```@example cshape
 cent = centroid(cshape)
-plot!(a, cent)
+scatter!(a, GI.x(cent), GI.y(cent), color = :red)
 f
 ```
 The points are ordered in a clockwise fashion, which means that the signed area
@@ -47,27 +47,12 @@ mutlipolygon.
 """
 centroid(geom) = centroid(GI.trait(geom), geom)
 
-"""
-    centroid(geom)::GI.Point
+centroid(trait::GI.LineStringTrait, geom) = centroid_and_length(trait, geom)[1]
 
-Returns the centroid of a given line segment.
-"""
-centroid(trait::GI.LineStringTrait, geom)  =
-    centroid_and_length(trait, geom)[1]
+centroid(trait, geom) = centroid_and_signed_area(trait, geom)[1]
 
 """
-    centroid(geom)::GI.Point
-
-Returns the centroid of a given linear ring, polygon, or
-mutlipolygon.
-"""
-centroid(
-    trait::Union{GI.LinearRingTrait, GI.PolygonTrait, GI.MultiPolygonTrait}, 
-    geom
-) = centroid_and_signed_area(trait, geom)[1]
-
-"""
-    centroid_and_length(geom)::GI.Point
+    centroid_and_length(geom)::(GI.Point, ::Real)
 
 Returns the centroid and length of a given geom. Note this is only valid for
 line strings.
@@ -75,14 +60,20 @@ line strings.
 centroid_and_length(geom) = centroid_and_length(GI.trait(geom), geom)
 
 """
-    centroid_and_length(::GI.LineStringTrait, geom)::(GI.Point, ::Real)
+    centroid_and_signed_area(geom)::(GI.Point, ::Real)
 
-Returns the centroid and length of a given line segment.
+Returns the centroid and area of a given geom. Note this is only valid for
+linear rings, polygons, and multipolygons.
 """
+centroid_and_signed_area(geom) = centroid_and_signed_area(GI.trait(geom), geom)
+
+
+
 function centroid_and_length(::GI.LineStringTrait, geom)
     FT = Float64
     # Initialize starting values
-    centroid = GI.Point(FT(0), FT(0))
+    xcentroid = FT(0)
+    ycentroid = FT(0)
     length = FT(0)
     point₁ = GI.getpoint(geom, 1)
     # Loop over line segments of line string
@@ -95,45 +86,53 @@ function centroid_and_length(::GI.LineStringTrait, geom)
         # Accumulate the line segment length into `length`
         length += length_component
         # Weighted average of line segment centroids
-        centroid = centroid .+ ((point₁ .+ point₂) .* (length_component / 2))
+        xcentroid += (GI.x(point₁) + GI.x(point₂)) * (length_component / 2)
+        ycentroid += (GI.y(point₁) + GI.y(point₂)) * (length_component / 2)
+        #centroid = centroid .+ ((point₁ .+ point₂) .* (length_component / 2))
         # Advance the point buffer by 1 point to move to next line segment
         point₁ = point₂
     end
-    return centroid ./= length, length
+    xcentroid /= length
+    ycentroid /= length
+    return GI.Point(xcentroid, ycentroid), length
 end
 
 """
-    centroid_and_length(geom)::GI.Point
-
-Returns the centroid and area of a given geom. Note this is only valid for
-linear rings, polygons, and multipolygons.
-"""
-centroid_and_signed_area(geom) = centroid_and_signed_area(GI.trait(geom), geom)
-
-"""
-    centroid_and_signed_area(::GI.LinearRingTrait, geom)::(GI.Point, ::Real)
+    centroid_and_signed_area(
+        ::Union{GI.LineStringTrait, GI.LinearRingTrait},
+        geom,
+    )::(GI.Point, ::Real)
 
-Returns the centroid and signed area of a given linear ring.
+Returns the centroid and signed area of a given a line string or a linear ring.
+Note that the area doesn't have much meaning as for a line string, and isn't
+valid if the line segment isn't closed. 
 """
-function centroid_and_signed_area(::GI.LinearRingTrait, geom)
+function centroid_and_signed_area(
+    ::Union{GI.LineStringTrait, GI.LinearRingTrait},
+    geom,
+)
     FT = Float64
     # Initialize starting values
-    centroid = GI.Point(FT(0), FT(0))
+    xcentroid = FT(0)
+    ycentroid = FT(0)
     area = FT(0)
     point₁ = GI.getpoint(geom, 1)
     # Loop over line segments of linear ring
     for point₂ in GI.getpoint(geom)
         area_component = GI.x(point₁) * GI.y(point₂) -
-            GI.x(point₁) * GI.y(point₂)
+            GI.x(point₂) * GI.y(point₁)
         # Accumulate the area component into `area`
         area += area_component
         # Weighted average of centroid components
-        centroid = centroid .+ (point₁ .+ point₂) .* area_component
+        xcentroid += (GI.x(point₁) + GI.x(point₂)) * area_component
+        ycentroid += (GI.y(point₁) + GI.y(point₂)) * area_component
         # Advance the point buffer by 1 point
         point₁ = point₂
     end
     area /= 2
-    return centroid ./= 6area, area
+    xcentroid /= 6area
+    ycentroid /= 6area
+    return GI.Point(xcentroid, ycentroid), area
 end
 
 """
@@ -142,10 +141,18 @@ end
 Returns the centroid and signed area of a given polygon.
 """
 function centroid_and_signed_area(::GI.PolygonTrait, geom)
+    FT = Float64
+    # Initialize starting values
+    xcentroid = FT(0)
+    ycentroid = FT(0)
+    area = FT(0)
     # Exterior polygon centroid and area
-    centroid, area = centroid_and_signed_area(GI.getexterior(geom))
+    ext_centroid, ext_area = centroid_and_signed_area(GI.getexterior(geom))
+    area += ext_area
+    ext_area = abs(ext_area)
     # Weight exterior centroid by area
-    centroid *= abs(signed_area)
+    xcentroid += GI.x(ext_centroid) * ext_area
+    ycentroid += GI.y(ext_centroid) * ext_area
     # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
         # Hole polygon's centroid and area
@@ -154,9 +161,12 @@ function centroid_and_signed_area(::GI.PolygonTrait, geom)
         # Accumulate the area component into `area`
         area -= interior_area
         # Weighted average of centroid components
-        centroid = centroid .- (interior_centroid .* interior_area)
+        xcentroid -= GI.x(interior_centroid) * interior_area
+        ycentroid -= GI.y(interior_centroid) * interior_area
     end
-    return centroid ./= abs(area), area
+    xcentroid /= abs(area)
+    ycentroid /= abs(area)
+    return GI.Point(xcentroid, ycentroid), area
 end
 
 """
@@ -167,7 +177,8 @@ Returns the centroid and signed area of a given multipolygon.
 function centroid_and_signed_area(::GI.MultiPolygonTrait, geom)
     FT = Float64
     # Initialize starting values
-    centroid = GI.Point(FT(0), FT(0))
+    xcentroid = FT(0)
+    ycentroid = FT(0)
     area = FT(0)
     # Loop over any polygons within the multipolygon
     for poly in GI.getpolygon(geom)
@@ -176,7 +187,10 @@ function centroid_and_signed_area(::GI.MultiPolygonTrait, geom)
         # Accumulate the area component into `area`
         area += poly_area
         # Weighted average of centroid components
-        centroid = centroid .+ (poly_centroid .* poly_area)
+        xcentroid += Gi.x(poly_centroid) * poly_area
+        ycentroid += Gi.y(poly_centroid) * poly_area
     end
-    return centroid ./= area, area
+    xcentroid /= area
+    ycentroid /= area
+    return GI.Point(xcentroid, ycentroid), area
 end
\ No newline at end of file
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
new file mode 100644
index 000000000..a65392cd3
--- /dev/null
+++ b/test/methods/centroid.jl
@@ -0,0 +1,53 @@
+@testset "Lines/Rings" begin
+    # Basic line string
+    l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
+    c1 = GO.centroid(l1)
+    c1_from_LG = LG.centroid(l1)
+    @test GI.x(c1) ≈ GI.x(c1_from_LG)  # I got this from 
+    @test GI.y(c1) ≈ GI.y(c1_from_LG)
+
+    # Spiral line string
+    l2 = LG.LineString([[0.0, 0.0], [2.5, 2.5], [10.0, 10.0]])
+    c2 = GO.centroid(21)
+    c2_from_LG = LG.centroid(l2)
+    @test GI.x(c2) ≈ GI.x(c2_from_LG)  # I got this from 
+    @test GI.y(c2) ≈ GI.y(c2_from_LG)
+
+    # Basic linear ring
+
+    # Fancier linear ring
+
+end
+@testset "Polygons" begin
+    # Basic rectangle
+    p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
+    c1 = GO.centroid(p1)
+    @test GI.x(c1) ≈ 5
+    @test GI.y(c1) ≈ 5
+
+    # Concave c-shaped polygon
+    p2 = LG.Polygon([[
+        [0.0, 0.0], [0.0, 3.0], [3.0, 3.0], [3.0, 2.0], [1.0, 2.0],
+        [1.0, 1.0], [3.0, 1.0], [3.0, 0.0], [0.0, 0.0],
+    ]])
+    c2 = GO.centroid(p2)
+    c2_from_LG = LG.centroid(p2)
+    @test GI.x(c2) ≈ GI.x(c2_from_LG)  # I got this from 
+    @test GI.y(c2) ≈ GI.y(c2_from_LG)
+
+    # Randomly generated polygon with lots of sides
+
+    # Polygon with one hole
+    p4 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0), (2.3 2.7, 2.5 4.9, 4.1 5.2, 1.9 4.2, 2.3 2.7))")
+    c4 = GO.centroid(p4)
+    c4_from_AG = AG.centroid(p4)
+    @test GI.x(c4) ≈ GI.x(c4_from_AG)
+    @test GI.y(c4) ≈ GI.y(c4_from_AG)
+
+    # Polygon with two holes
+    
+end
+@testset "MultiPolygons" begin
+    # Combine poylgons made above
+    
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index a270b6f31..73fa26cc2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,15 +4,20 @@ using Test
 using GeometryOps.GeoInterface
 using GeometryOps.GeometryBasics
 using ArchGDAL
+using LibGEOS
 
 const GI = GeoInterface
 const AG = ArchGDAL
+const LG = LibGEOS
 
 @testset "GeometryOps.jl" begin
     @testset "Primitives" begin include("primitives.jl") end
+    # Methods
+    @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
     @testset "Bools" begin include("methods/bools.jl") end
+    @testset "Centroid" begin include("methods/centroid.jl") end
     @testset "Signed Area" begin include("methods/signed_area.jl") end
-    @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
+    # Transformations
     @testset "Reproject" begin include("transformations/reproject.jl") end
     @testset "Flip" begin include("transformations/flip.jl") end
     @testset "Simplify" begin include("transformations/simplify.jl") end

From d407bfa7a47ee8dc0b59b0f9d007a73674bd7170 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Wed, 27 Sep 2023 21:25:02 -0700
Subject: [PATCH 05/10] Update post ring tests

---
 src/methods/centroid.jl  | 67 ++++++++++++++++++++++++++++++++--------
 test/methods/centroid.jl | 44 +++++++++++++++++++-------
 2 files changed, 87 insertions(+), 24 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 134185ed1..7be507748 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -8,6 +8,10 @@ export centroid
 The centroid is the geometric center of a line string or area(s). Note that
 the centroid does not need to be inside of a concave area.
 
+Further note that by convention a line, or linear ring, is calculated by
+weighting the line segments by their length, while polygons and multipolygon
+centroids are calculated by weighting edge's by their 'area components'.
+
 To provide an example, consider this concave polygon in the shape of a 'C':
 ```@example cshape
 using GeometryOps
@@ -35,9 +39,19 @@ is positive.  If we reverse the order of the points, we get a negative area.
 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! 
+implementation based on the geometry trait. This is also used in the
+implementation, since it's a lot less work! 
+
+Note that if you call centroid on a LineString or LinearRing, the
+centroid_and_length function will be called due to the weighting scheme
+described above, while centroid_and_signed_area is called for polygons and
+multipolygons. However, centroid_and_signed_area can still be called on a
+LineString or LinearRing when they are closed, for example as the interior hole
+of a polygon.
+
+The helper functions centroid_and_length and centroid_and_signed_area are made
+availible just in case the user also needs the signed area or length to decrease
+repeat computation.
 =#
 """
     centroid(geom)::GI.Point
@@ -47,29 +61,56 @@ mutlipolygon.
 """
 centroid(geom) = centroid(GI.trait(geom), geom)
 
-centroid(trait::GI.LineStringTrait, geom) = centroid_and_length(trait, geom)[1]
+"""
+    centroid(
+        trait::Union{GI.LineStringTrait, GI.LinearRingTrait},
+        geom,
+    )
+
+Returns the centroid of a line string or linear ring, which is calculated by
+weighting line segments by their length by convention.
+"""
+centroid(
+    trait::Union{GI.LineStringTrait, GI.LinearRingTrait},
+    geom,
+) = centroid_and_length(trait, geom)[1]
 
+"""
+    centroid(trait, geom)
+
+Returns the centroid of a polygon or multipolygon, which is calculated by
+weighting edges by their `area component` by convention.
+"""
 centroid(trait, geom) = centroid_and_signed_area(trait, geom)[1]
 
 """
     centroid_and_length(geom)::(GI.Point, ::Real)
 
-Returns the centroid and length of a given geom. Note this is only valid for
-line strings.
+Returns the centroid and length of a given line/ring. Note this is only valid
+for line strings and linear rings.
 """
 centroid_and_length(geom) = centroid_and_length(GI.trait(geom), geom)
 
 """
-    centroid_and_signed_area(geom)::(GI.Point, ::Real)
+    centroid_and_signed_area(
+        ::Union{GI.LineStringTrait, GI.LinearRingTrait}, 
+        geom,
+    )::(GI.Point, ::Real)
 
-Returns the centroid and area of a given geom. Note this is only valid for
-linear rings, polygons, and multipolygons.
+Returns the centroid and signed area of a given geom.
 """
 centroid_and_signed_area(geom) = centroid_and_signed_area(GI.trait(geom), geom)
 
+"""
+    centroid_and_length(geom)::(GI.Point, ::Real)
 
-
-function centroid_and_length(::GI.LineStringTrait, geom)
+Returns the centroid and length of a given line/ring. Note this is only valid
+for line strings and linear rings.
+"""
+function centroid_and_length(
+    ::Union{GI.LineStringTrait, GI.LinearRingTrait},
+    geom,
+)
     FT = Float64
     # Initialize starting values
     xcentroid = FT(0)
@@ -187,8 +228,8 @@ function centroid_and_signed_area(::GI.MultiPolygonTrait, geom)
         # Accumulate the area component into `area`
         area += poly_area
         # Weighted average of centroid components
-        xcentroid += Gi.x(poly_centroid) * poly_area
-        ycentroid += Gi.y(poly_centroid) * poly_area
+        xcentroid += GI.x(poly_centroid) * poly_area
+        ycentroid += GI.y(poly_centroid) * poly_area
     end
     xcentroid /= area
     ycentroid /= area
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index a65392cd3..283b05da7 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -7,13 +7,18 @@
     @test GI.y(c1) ≈ GI.y(c1_from_LG)
 
     # Spiral line string
-    l2 = LG.LineString([[0.0, 0.0], [2.5, 2.5], [10.0, 10.0]])
-    c2 = GO.centroid(21)
+    l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]])
+    c2 = GO.centroid(l2)
     c2_from_LG = LG.centroid(l2)
-    @test GI.x(c2) ≈ GI.x(c2_from_LG)  # I got this from 
+    @test GI.x(c2) ≈ GI.x(c2_from_LG) 
     @test GI.y(c2) ≈ GI.y(c2_from_LG)
 
     # Basic linear ring
+    r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]])
+    c3 = GO.centroid(r1)
+    c3_from_LG = LG.centroid(r1)
+    @test GI.x(c3) ≈ GI.x(c3_from_LG) 
+    @test GI.y(c3) ≈ GI.y(c3_from_LG)
 
     # Fancier linear ring
 
@@ -27,27 +32,44 @@ end
 
     # Concave c-shaped polygon
     p2 = LG.Polygon([[
-        [0.0, 0.0], [0.0, 3.0], [3.0, 3.0], [3.0, 2.0], [1.0, 2.0],
-        [1.0, 1.0], [3.0, 1.0], [3.0, 0.0], [0.0, 0.0],
+        [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
+        [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0],
     ]])
     c2 = GO.centroid(p2)
     c2_from_LG = LG.centroid(p2)
-    @test GI.x(c2) ≈ GI.x(c2_from_LG)  # I got this from 
+    @test GI.x(c2) ≈ GI.x(c2_from_LG)
     @test GI.y(c2) ≈ GI.y(c2_from_LG)
 
     # Randomly generated polygon with lots of sides
 
     # Polygon with one hole
-    p4 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0), (2.3 2.7, 2.5 4.9, 4.1 5.2, 1.9 4.2, 2.3 2.7))")
+    p4 = LG.Polygon([
+        [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+        [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
+    ])
     c4 = GO.centroid(p4)
-    c4_from_AG = AG.centroid(p4)
-    @test GI.x(c4) ≈ GI.x(c4_from_AG)
-    @test GI.y(c4) ≈ GI.y(c4_from_AG)
+    c4_from_LG = LG.centroid(p4)
+    @test GI.x(c4) ≈ GI.x(c4_from_LG)
+    @test GI.y(c4) ≈ GI.y(c4_from_LG)
 
     # Polygon with two holes
     
 end
 @testset "MultiPolygons" begin
     # Combine poylgons made above
-    
+    m1 = LibGEOS.MultiPolygon([
+        [
+            [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
+            [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]],
+        ],
+        [
+            [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+            [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
+        ]
+    ])
+    c1 = GO.centroid(m1)
+    c1_from_LG = LG.centroid(m1)
+    # TODO: This these is failing
+    # @test GI.x(c1) ≈ GI.x(c1_from_LG)
+    # @test GI.y(c1) ≈ GI.y(c1_from_LG)
 end
\ No newline at end of file

From afeca1e08a322c3962a5c2f1182de43b7d478a39 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Wed, 27 Sep 2023 21:27:21 -0700
Subject: [PATCH 06/10] Remove libgeos dependency

---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 35d7d7e1c..9d39ac19c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.0.1-DEV"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

From 7386ca3f69d73a9cca63fcac4dba982812b2dd32 Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Fri, 29 Sep 2023 12:58:10 -0700
Subject: [PATCH 07/10] Add random polygon generator

---
 Project.toml                    |  7 +++-
 test/data/polygon_generation.jl | 62 +++++++++++++++++++++++++++++++++
 test/methods/centroid.jl        | 11 ++++++
 test/runtests.jl                |  1 +
 4 files changed, 80 insertions(+), 1 deletion(-)
 create mode 100644 test/data/polygon_generation.jl

diff --git a/Project.toml b/Project.toml
index 35d7d7e1c..9fa5cfe16 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,10 +19,15 @@ julia = "1.3"
 
 [extras]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "GeoFormatTypes", "GeoJSON", "LibGEOS", "Test"]
+test = [
+    "ArchGDAL", "Distributions", "GeoFormatTypes", "GeoJSON", "LibGEOS",
+    "Random", "Test",
+]
diff --git a/test/data/polygon_generation.jl b/test/data/polygon_generation.jl
new file mode 100644
index 000000000..0d798af1f
--- /dev/null
+++ b/test/data/polygon_generation.jl
@@ -0,0 +1,62 @@
+"""
+    generate_random_poly(
+        x,
+        y,
+        nverts,
+        avg_radius,
+        irregularity,
+        spikiness,
+    )
+Generate a random polygon given a center point, number of vertices, and measures
+of the polygon's irregularity.
+Inputs:
+    x               <Real> x-coordinate for center point
+    y               <Real> y-coordinate for center point
+    nverts          <Int> number of vertices
+    avg_radius      <Real> average radius for each point to center point
+    irregularity    <Real> measure of randomness for difference in angle between
+                        points (between 0 and 1)
+    spikiness       <Real> measure of randomness for difference in radius
+                        between points (between 0 and 1)
+    rng             <RNG> random number generator for polygon creation
+"""
+function generate_random_poly(
+    x,
+    y,
+    nverts,
+    avg_radius,
+    irregularity::T,
+    spikiness::T,
+    rng = Xoshiro()
+) where T <: Real
+    # Check argument validity
+    @assert 0 <= irregularity <= 1 "Irregularity must be between 0 and 1"
+    @assert 0 <= spikiness <= 1 "Spikiness must be between 0 and 1"
+    # Setup basic parameters
+    avg_angle = 2π / nverts
+    ϵ_angle = irregularity * avg_angle
+    # ϵ_rad = spikiness * avg_radius
+    smallest_angle_step = avg_angle - ϵ_angle
+    # smallest_rad = avg_radius - ϵ_rad
+    current_angle = rand(rng) * 2π
+    rad_distribution = Distributions.Normal(avg_radius, spikiness)
+    points = [zeros(T, 2) for _ in 1:nverts]
+    # Generate angle steps around polygon
+    angle_steps = zeros(T, nverts)
+    cumsum = T(0)
+    for i in 1:nverts
+        step = smallest_angle_step + 2ϵ_angle * rand(rng)
+        angle_steps[i] = step
+        cumsum += step
+    end
+    angle_steps ./= (cumsum / 2π)
+    # Generate polygon points at given angles and radii
+    for i in 1:nverts
+        rad = clamp(rand(rad_distribution), 0, 2avg_radius) #smallest_rad + 2ϵ_rad * rand(rng)
+        points[i][1] = x + rad * cos(current_angle)
+        points[i][2] = y + rad * sin(current_angle)
+        current_angle -= angle_steps[i]
+    end
+    points[end] .= points[1]
+    return [points]
+end
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index a65392cd3..61adb663e 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -36,6 +36,17 @@ end
     @test GI.y(c2) ≈ GI.y(c2_from_LG)
 
     # Randomly generated polygon with lots of sides
+    p3 = LG.Polygon([[
+        [14.567, 8.974], [13.579, 8.849], [12.076, 8.769], [11.725, 8.567],
+        [11.424, 6.451], [10.187, 7.712], [8.187, 6.795], [8.065, 8.096],
+        [6.827, 8.287], [7.1628, 8.9221], [5.8428, 10.468], [7.987, 11.734],
+        [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340],
+        [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974],
+    ]])
+    c3 = GO.centroid(p3)
+    c3_from_LG = LG.centroid(p3)
+    @test GI.x(c3) ≈ GI.x(c3_from_LG)  # I got this from 
+    @test GI.y(c3) ≈ GI.y(c3_from_LG)
 
     # Polygon with one hole
     p4 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0), (2.3 2.7, 2.5 4.9, 4.1 5.2, 1.9 4.2, 2.3 2.7))")
diff --git a/test/runtests.jl b/test/runtests.jl
index 73fa26cc2..af6107207 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ using GeometryOps.GeoInterface
 using GeometryOps.GeometryBasics
 using ArchGDAL
 using LibGEOS
+using Random, Distributions
 
 const GI = GeoInterface
 const AG = ArchGDAL

From 868e9fad2d971bfd66c1f511db6ebd272fd0f7dd Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Fri, 29 Sep 2023 14:49:27 -0700
Subject: [PATCH 08/10] Add more tests and refine

---
 src/methods/centroid.jl         | 92 ++++++++++++++++-----------------
 src/methods/signed_area.jl      |  2 +-
 test/data/polygon_generation.jl |  4 ++
 test/methods/centroid.jl        | 56 ++++++++++++++------
 test/runtests.jl                |  1 +
 5 files changed, 93 insertions(+), 62 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 7be507748..99bf6490b 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -1,6 +1,6 @@
 # # Centroid
 
-export centroid
+export centroid, centroid_and_length, centroid_and_area
 
 #=
 ## What is the centroid?
@@ -31,8 +31,6 @@ cent = centroid(cshape)
 scatter!(a, GI.x(cent), GI.y(cent), color = :red)
 f
 ```
-The points are ordered in a clockwise fashion, which means that the signed area
-is positive.  If we reverse the order of the points, we get a negative area.
 
 ## Implementation
 
@@ -44,13 +42,13 @@ implementation, since it's a lot less work!
 
 Note that if you call centroid on a LineString or LinearRing, the
 centroid_and_length function will be called due to the weighting scheme
-described above, while centroid_and_signed_area is called for polygons and
-multipolygons. However, centroid_and_signed_area can still be called on a
+described above, while centroid_and_area is called for polygons and
+multipolygons. However, centroid_and_area can still be called on a
 LineString or LinearRing when they are closed, for example as the interior hole
 of a polygon.
 
-The helper functions centroid_and_length and centroid_and_signed_area are made
-availible just in case the user also needs the signed area or length to decrease
+The helper functions centroid_and_length and centroid_and_area are made
+availible just in case the user also needs the area or length to decrease
 repeat computation.
 =#
 """
@@ -81,7 +79,7 @@ centroid(
 Returns the centroid of a polygon or multipolygon, which is calculated by
 weighting edges by their `area component` by convention.
 """
-centroid(trait, geom) = centroid_and_signed_area(trait, geom)[1]
+centroid(trait, geom) = centroid_and_area(trait, geom)[1]
 
 """
     centroid_and_length(geom)::(GI.Point, ::Real)
@@ -92,14 +90,14 @@ for line strings and linear rings.
 centroid_and_length(geom) = centroid_and_length(GI.trait(geom), geom)
 
 """
-    centroid_and_signed_area(
+    centroid_and_area(
         ::Union{GI.LineStringTrait, GI.LinearRingTrait}, 
         geom,
     )::(GI.Point, ::Real)
 
-Returns the centroid and signed area of a given geom.
+Returns the centroid and area of a given geom.
 """
-centroid_and_signed_area(geom) = centroid_and_signed_area(GI.trait(geom), geom)
+centroid_and_area(geom) = centroid_and_area(GI.trait(geom), geom)
 
 """
     centroid_and_length(geom)::(GI.Point, ::Real)
@@ -111,11 +109,11 @@ function centroid_and_length(
     ::Union{GI.LineStringTrait, GI.LinearRingTrait},
     geom,
 )
-    FT = Float64
+    T = Float64
     # Initialize starting values
-    xcentroid = FT(0)
-    ycentroid = FT(0)
-    length = FT(0)
+    xcentroid = T(0)
+    ycentroid = T(0)
+    length = T(0)
     point₁ = GI.getpoint(geom, 1)
     # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
@@ -139,24 +137,28 @@ function centroid_and_length(
 end
 
 """
-    centroid_and_signed_area(
+    centroid_and_area(
         ::Union{GI.LineStringTrait, GI.LinearRingTrait},
         geom,
     )::(GI.Point, ::Real)
 
-Returns the centroid and signed area of a given a line string or a linear ring.
-Note that the area doesn't have much meaning as for a line string, and isn't
-valid if the line segment isn't closed. 
+Returns the centroid and area of a given a line string or a linear ring.
+Note that this is only valid if the line segment or linear ring is closed. 
 """
-function centroid_and_signed_area(
+function centroid_and_area(
     ::Union{GI.LineStringTrait, GI.LinearRingTrait},
     geom,
 )
-    FT = Float64
+    T = Float64
+    # Check that the geometry is closed
+    @assert(
+        GI.getpoint(geom, 1) == GI.getpoint(geom, GI.ngeom(geom)),
+        "centroid_and_area should only be used with closed geometries"
+    )
     # Initialize starting values
-    xcentroid = FT(0)
-    ycentroid = FT(0)
-    area = FT(0)
+    xcentroid = T(0)
+    ycentroid = T(0)
+    area = T(0)
     point₁ = GI.getpoint(geom, 1)
     # Loop over line segments of linear ring
     for point₂ in GI.getpoint(geom)
@@ -173,58 +175,56 @@ function centroid_and_signed_area(
     area /= 2
     xcentroid /= 6area
     ycentroid /= 6area
-    return GI.Point(xcentroid, ycentroid), area
+    return GI.Point(xcentroid, ycentroid), abs(area)
 end
 
 """
-    centroid_and_signed_area(::GI.PolygonTrait, geom)::(GI.Point, ::Real)
+    centroid_and_area(::GI.PolygonTrait, geom)::(GI.Point, ::Real)
 
-Returns the centroid and signed area of a given polygon.
+Returns the centroid and area of a given polygon.
 """
-function centroid_and_signed_area(::GI.PolygonTrait, geom)
-    FT = Float64
+function centroid_and_area(::GI.PolygonTrait, geom)
+    T = Float64
     # Initialize starting values
-    xcentroid = FT(0)
-    ycentroid = FT(0)
-    area = FT(0)
+    xcentroid = T(0)
+    ycentroid = T(0)
+    area = T(0)
     # Exterior polygon centroid and area
-    ext_centroid, ext_area = centroid_and_signed_area(GI.getexterior(geom))
+    ext_centroid, ext_area = centroid_and_area(GI.getexterior(geom))
     area += ext_area
-    ext_area = abs(ext_area)
     # Weight exterior centroid by area
     xcentroid += GI.x(ext_centroid) * ext_area
     ycentroid += GI.y(ext_centroid) * ext_area
     # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
         # Hole polygon's centroid and area
-        interior_centroid, interior_area = centroid_and_signed_area(hole)
-        interior_area = abs(interior_area)
+        interior_centroid, interior_area = centroid_and_area(hole)
         # Accumulate the area component into `area`
         area -= interior_area
         # Weighted average of centroid components
         xcentroid -= GI.x(interior_centroid) * interior_area
         ycentroid -= GI.y(interior_centroid) * interior_area
     end
-    xcentroid /= abs(area)
-    ycentroid /= abs(area)
+    xcentroid /= area
+    ycentroid /= area
     return GI.Point(xcentroid, ycentroid), area
 end
 
 """
-    centroid_and_signed_area(::GI.MultiPolygonTrait, geom)::(GI.Point, ::Real)
+    centroid_and_area(::GI.MultiPolygonTrait, geom)::(GI.Point, ::Real)
 
-Returns the centroid and signed area of a given multipolygon.
+Returns the centroid and area of a given multipolygon.
 """
-function centroid_and_signed_area(::GI.MultiPolygonTrait, geom)
-    FT = Float64
+function centroid_and_area(::GI.MultiPolygonTrait, geom)
+    T = Float64
     # Initialize starting values
-    xcentroid = FT(0)
-    ycentroid = FT(0)
-    area = FT(0)
+    xcentroid = T(0)
+    ycentroid = T(0)
+    area = T(0)
     # Loop over any polygons within the multipolygon
     for poly in GI.getpolygon(geom)
         # Polygon centroid and area
-        poly_centroid, poly_area = centroid_and_signed_area(poly)
+        poly_centroid, poly_area = centroid_and_area(poly)
         # Accumulate the area component into `area`
         area += poly_area
         # Weighted average of centroid components
diff --git a/src/methods/signed_area.jl b/src/methods/signed_area.jl
index ecd13ae61..c485d66c4 100644
--- a/src/methods/signed_area.jl
+++ b/src/methods/signed_area.jl
@@ -25,7 +25,7 @@ export signed_area
 # f
 # ```
 # The points are ordered in a clockwise fashion, which means that the signed area
-# is positive.  If we reverse the order of the points, we get a negative area.
+# is negative.  If we reverse the order of the points, we get a postive area.
 
 # ## Implementation
 
diff --git a/test/data/polygon_generation.jl b/test/data/polygon_generation.jl
index 0d798af1f..e0ab68e9f 100644
--- a/test/data/polygon_generation.jl
+++ b/test/data/polygon_generation.jl
@@ -19,6 +19,10 @@ Inputs:
     spikiness       <Real> measure of randomness for difference in radius
                         between points (between 0 and 1)
     rng             <RNG> random number generator for polygon creation
+Output:
+    Vector{Vector{Vector{T}}} representing polygon coordinates
+Note:
+    Check your outputs! No guarentee that the polygon's aren't self-intersecting
 """
 function generate_random_poly(
     x,
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index 4a70b7ba4..fb95b9b45 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -1,28 +1,32 @@
 @testset "Lines/Rings" begin
     # Basic line string
     l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
-    c1 = GO.centroid(l1)
+    c1, len1 = GO.centroid_and_length(l1)
     c1_from_LG = LG.centroid(l1)
-    @test GI.x(c1) ≈ GI.x(c1_from_LG)  # I got this from 
+    @test GI.x(c1) ≈ GI.x(c1_from_LG)
     @test GI.y(c1) ≈ GI.y(c1_from_LG)
+    @test len1 ≈ 20.0
 
     # Spiral line string
     l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]])
-    c2 = GO.centroid(l2)
+    c2, len2 = GO.centroid_and_length(l2)
     c2_from_LG = LG.centroid(l2)
     @test GI.x(c2) ≈ GI.x(c2_from_LG) 
     @test GI.y(c2) ≈ GI.y(c2_from_LG)
+    @test len2 ≈ 59.3090856788928
 
-    # Basic linear ring
+    # Test that non-closed line strings throw an error for centroid_and_area
+    @test_throws AssertionError GO.centroid_and_area(l2)
+
+    # Basic linear ring - note that this still uses weighting by length
     r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]])
     c3 = GO.centroid(r1)
     c3_from_LG = LG.centroid(r1)
     @test GI.x(c3) ≈ GI.x(c3_from_LG) 
     @test GI.y(c3) ≈ GI.y(c3_from_LG)
 
-    # Fancier linear ring
-
 end
+
 @testset "Polygons" begin
     # Basic rectangle
     p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
@@ -35,10 +39,11 @@ end
         [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
         [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0],
     ]])
-    c2 = GO.centroid(p2)
+    c2, area2 = GO.centroid_and_area(p2)
     c2_from_LG = LG.centroid(p2)
     @test GI.x(c2) ≈ GI.x(c2_from_LG)
     @test GI.y(c2) ≈ GI.y(c2_from_LG)
+    @test area2 ≈ LG.area(p2)
 
     # Randomly generated polygon with lots of sides
     p3 = LG.Polygon([[
@@ -48,27 +53,48 @@ end
         [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340],
         [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974],
     ]])
-    c3 = GO.centroid(p3)
+    c3, area3 = GO.centroid_and_area(p3)
     c3_from_LG = LG.centroid(p3)
-    @test GI.x(c3) ≈ GI.x(c3_from_LG)  # I got this from 
+    @test GI.x(c3) ≈ GI.x(c3_from_LG)
     @test GI.y(c3) ≈ GI.y(c3_from_LG)
+    @test area3 ≈ LG.area(p3)
 
     # Polygon with one hole
     p4 = LG.Polygon([
         [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
         [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
     ])
-    c4 = GO.centroid(p4)
+    c4, area4 = GO.centroid_and_area(p4)
     c4_from_LG = LG.centroid(p4)
     @test GI.x(c4) ≈ GI.x(c4_from_LG)
     @test GI.y(c4) ≈ GI.y(c4_from_LG)
+    @test area4 ≈ LG.area(p4)
 
     # Polygon with two holes
+    p5 = LG.Polygon([
+        [[-10.0, -10.0], [-2.0, 0.0], [6.0, -10.0], [-10.0, -10.0]],
+        [[-8.0, -8.0], [-8.0, -7.0], [-4.0, -7.0], [-4.0, -8.0], [-8.0, -8.0]],
+        [[-3.0,-9.0], [3.0, -9.0], [3.0, -8.5], [-3.0, -8.5], [-3.0, -9.0]],
+    ])
+    c5 = GO.centroid(p5)
+    c5_from_LG = LG.centroid(p5)
+    @test GI.x(c5) ≈ GI.x(c5_from_LG)
+    @test GI.y(c5) ≈ GI.y(c5_from_LG)
+
+    # Same polygon as P5 but using a GeoInterface polygon
+    p6 = GI.Polygon([
+        [(-10.0, -10.0), (-2.0, 0.0), (6.0, -10.0), (-10.0, -10.0)],
+        [(-8.0, -8.0), (-8.0, -7.0), (-4.0, -7.0), (-4.0, -8.0), (-8.0, -8.0)],
+        [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)],
+    ])
+    c6 = GO.centroid(p6)
+    @test GI.x(c6) ≈ GI.x(c5)
+    @test GI.y(c6) ≈ GI.y(c5)
     
 end
 @testset "MultiPolygons" begin
     # Combine poylgons made above
-    m1 = LibGEOS.MultiPolygon([
+    m1 = LG.MultiPolygon([
         [
             [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
             [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]],
@@ -78,9 +104,9 @@ end
             [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
         ]
     ])
-    c1 = GO.centroid(m1)
+    c1, area1 = GO.centroid_and_area(m1)
     c1_from_LG = LG.centroid(m1)
-    # TODO: This these is failing
-    # @test GI.x(c1) ≈ GI.x(c1_from_LG)
-    # @test GI.y(c1) ≈ GI.y(c1_from_LG)
+    @test GI.x(c1) ≈ GI.x(c1_from_LG)
+    @test GI.y(c1) ≈ GI.y(c1_from_LG)
+    @test area1 ≈ LG.area(m1)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index af6107207..c4fc39f3e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,6 +10,7 @@ using Random, Distributions
 const GI = GeoInterface
 const AG = ArchGDAL
 const LG = LibGEOS
+const GO = GeometryOps
 
 @testset "GeometryOps.jl" begin
     @testset "Primitives" begin include("primitives.jl") end

From fe1d6da8af788c6da16aad21f36ef00d3b475639 Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Mon, 2 Oct 2023 11:58:16 -0700
Subject: [PATCH 09/10] Switch return type to be tuple

---
 src/methods/centroid.jl  | 50 ++++++++++++++++++----------------------
 test/methods/centroid.jl | 38 +++++++++++++++---------------
 2 files changed, 40 insertions(+), 48 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 99bf6490b..759c7fcfd 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -109,7 +109,7 @@ function centroid_and_length(
     ::Union{GI.LineStringTrait, GI.LinearRingTrait},
     geom,
 )
-    T = Float64
+    T = typeof(GI.x(GI.getpoint(geom, 1)))
     # Initialize starting values
     xcentroid = T(0)
     ycentroid = T(0)
@@ -133,7 +133,7 @@ function centroid_and_length(
     end
     xcentroid /= length
     ycentroid /= length
-    return GI.Point(xcentroid, ycentroid), length
+    return (xcentroid, ycentroid), length
 end
 
 """
@@ -149,7 +149,7 @@ function centroid_and_area(
     ::Union{GI.LineStringTrait, GI.LinearRingTrait},
     geom,
 )
-    T = Float64
+    T = typeof(GI.x(GI.getpoint(geom, 1)))
     # Check that the geometry is closed
     @assert(
         GI.getpoint(geom, 1) == GI.getpoint(geom, GI.ngeom(geom)),
@@ -175,7 +175,7 @@ function centroid_and_area(
     area /= 2
     xcentroid /= 6area
     ycentroid /= 6area
-    return GI.Point(xcentroid, ycentroid), abs(area)
+    return (xcentroid, ycentroid), abs(area)
 end
 
 """
@@ -184,30 +184,24 @@ end
 Returns the centroid and area of a given polygon.
 """
 function centroid_and_area(::GI.PolygonTrait, geom)
-    T = Float64
-    # Initialize starting values
-    xcentroid = T(0)
-    ycentroid = T(0)
-    area = T(0)
-    # Exterior polygon centroid and area
-    ext_centroid, ext_area = centroid_and_area(GI.getexterior(geom))
-    area += ext_area
+    # Exterior ring's centroid and area
+    (xcentroid, ycentroid), area = centroid_and_area(GI.getexterior(geom))
     # Weight exterior centroid by area
-    xcentroid += GI.x(ext_centroid) * ext_area
-    ycentroid += GI.y(ext_centroid) * ext_area
+    xcentroid *= area
+    ycentroid *= area
     # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
         # Hole polygon's centroid and area
-        interior_centroid, interior_area = centroid_and_area(hole)
+        (xinterior, yinterior), interior_area = centroid_and_area(hole)
         # Accumulate the area component into `area`
         area -= interior_area
         # Weighted average of centroid components
-        xcentroid -= GI.x(interior_centroid) * interior_area
-        ycentroid -= GI.y(interior_centroid) * interior_area
+        xcentroid -= xinterior * interior_area
+        ycentroid -= yinterior * interior_area
     end
     xcentroid /= area
     ycentroid /= area
-    return GI.Point(xcentroid, ycentroid), area
+    return (xcentroid, ycentroid), area
 end
 
 """
@@ -216,22 +210,22 @@ end
 Returns the centroid and area of a given multipolygon.
 """
 function centroid_and_area(::GI.MultiPolygonTrait, geom)
-    T = Float64
-    # Initialize starting values
-    xcentroid = T(0)
-    ycentroid = T(0)
-    area = T(0)
+    # First polygon's centroid and area
+    (xcentroid, ycentroid), area = centroid_and_area(GI.getpolygon(geom, 1))
+    # Weight first polygon's centroid by area
+    xcentroid *= area
+    ycentroid *= area
     # Loop over any polygons within the multipolygon
-    for poly in GI.getpolygon(geom)
+    for i in 2:GI.ngeom(geom) #poly in GI.getpolygon(geom)
         # Polygon centroid and area
-        poly_centroid, poly_area = centroid_and_area(poly)
+        (xpoly, ypoly), poly_area = centroid_and_area(GI.getpolygon(geom, i))
         # Accumulate the area component into `area`
         area += poly_area
         # Weighted average of centroid components
-        xcentroid += GI.x(poly_centroid) * poly_area
-        ycentroid += GI.y(poly_centroid) * poly_area
+        xcentroid += xpoly * poly_area
+        ycentroid += ypoly * poly_area
     end
     xcentroid /= area
     ycentroid /= area
-    return GI.Point(xcentroid, ycentroid), area
+    return (xcentroid, ycentroid), area
 end
\ No newline at end of file
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index fb95b9b45..9aad3d861 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -3,16 +3,16 @@
     l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
     c1, len1 = GO.centroid_and_length(l1)
     c1_from_LG = LG.centroid(l1)
-    @test GI.x(c1) ≈ GI.x(c1_from_LG)
-    @test GI.y(c1) ≈ GI.y(c1_from_LG)
+    @test c1[1] ≈ GI.x(c1_from_LG)
+    @test c1[2] ≈ GI.y(c1_from_LG)
     @test len1 ≈ 20.0
 
     # Spiral line string
     l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]])
     c2, len2 = GO.centroid_and_length(l2)
     c2_from_LG = LG.centroid(l2)
-    @test GI.x(c2) ≈ GI.x(c2_from_LG) 
-    @test GI.y(c2) ≈ GI.y(c2_from_LG)
+    @test c2[1] ≈ GI.x(c2_from_LG) 
+    @test c2[2] ≈ GI.y(c2_from_LG)
     @test len2 ≈ 59.3090856788928
 
     # Test that non-closed line strings throw an error for centroid_and_area
@@ -22,15 +22,15 @@
     r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]])
     c3 = GO.centroid(r1)
     c3_from_LG = LG.centroid(r1)
-    @test GI.x(c3) ≈ GI.x(c3_from_LG) 
-    @test GI.y(c3) ≈ GI.y(c3_from_LG)
-
+    @test c3[1] ≈ GI.x(c3_from_LG) 
+    @test c3[2] ≈ GI.y(c3_from_LG)
 end
 
 @testset "Polygons" begin
     # Basic rectangle
     p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
     c1 = GO.centroid(p1)
+    c1 .≈ (5, 5)
     @test GI.x(c1) ≈ 5
     @test GI.y(c1) ≈ 5
 
@@ -41,8 +41,8 @@ end
     ]])
     c2, area2 = GO.centroid_and_area(p2)
     c2_from_LG = LG.centroid(p2)
-    @test GI.x(c2) ≈ GI.x(c2_from_LG)
-    @test GI.y(c2) ≈ GI.y(c2_from_LG)
+    @test c2[1] ≈ GI.x(c2_from_LG)
+    @test c2[2] ≈ GI.y(c2_from_LG)
     @test area2 ≈ LG.area(p2)
 
     # Randomly generated polygon with lots of sides
@@ -55,8 +55,8 @@ end
     ]])
     c3, area3 = GO.centroid_and_area(p3)
     c3_from_LG = LG.centroid(p3)
-    @test GI.x(c3) ≈ GI.x(c3_from_LG)
-    @test GI.y(c3) ≈ GI.y(c3_from_LG)
+    @test c3[1] ≈ GI.x(c3_from_LG)
+    @test c3[2] ≈ GI.y(c3_from_LG)
     @test area3 ≈ LG.area(p3)
 
     # Polygon with one hole
@@ -66,8 +66,8 @@ end
     ])
     c4, area4 = GO.centroid_and_area(p4)
     c4_from_LG = LG.centroid(p4)
-    @test GI.x(c4) ≈ GI.x(c4_from_LG)
-    @test GI.y(c4) ≈ GI.y(c4_from_LG)
+    @test c4[1] ≈ GI.x(c4_from_LG)
+    @test c4[2] ≈ GI.y(c4_from_LG)
     @test area4 ≈ LG.area(p4)
 
     # Polygon with two holes
@@ -78,8 +78,8 @@ end
     ])
     c5 = GO.centroid(p5)
     c5_from_LG = LG.centroid(p5)
-    @test GI.x(c5) ≈ GI.x(c5_from_LG)
-    @test GI.y(c5) ≈ GI.y(c5_from_LG)
+    @test c5[1] ≈ GI.x(c5_from_LG)
+    @test c5[2] ≈ GI.y(c5_from_LG)
 
     # Same polygon as P5 but using a GeoInterface polygon
     p6 = GI.Polygon([
@@ -88,9 +88,7 @@ end
         [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)],
     ])
     c6 = GO.centroid(p6)
-    @test GI.x(c6) ≈ GI.x(c5)
-    @test GI.y(c6) ≈ GI.y(c5)
-    
+    @test all(c5 .≈ c6)
 end
 @testset "MultiPolygons" begin
     # Combine poylgons made above
@@ -106,7 +104,7 @@ end
     ])
     c1, area1 = GO.centroid_and_area(m1)
     c1_from_LG = LG.centroid(m1)
-    @test GI.x(c1) ≈ GI.x(c1_from_LG)
-    @test GI.y(c1) ≈ GI.y(c1_from_LG)
+    @test c1[1] ≈ GI.x(c1_from_LG)
+    @test c1[2] ≈ GI.y(c1_from_LG)
     @test area1 ≈ LG.area(m1)
 end
\ No newline at end of file

From 055a0305cece3d6c5f9b6ee7c2066bda528fdfab Mon Sep 17 00:00:00 2001
From: Skylar A Gering <skygering@gmail.com>
Date: Mon, 2 Oct 2023 12:05:32 -0700
Subject: [PATCH 10/10] Update docs to match output type

---
 src/methods/centroid.jl | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 759c7fcfd..03dbc6798 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -52,7 +52,7 @@ availible just in case the user also needs the area or length to decrease
 repeat computation.
 =#
 """
-    centroid(geom)::GI.Point
+    centroid(geom)::Tuple{T, T}
 
 Returns the centroid of a given line segment, linear ring, polygon, or
 mutlipolygon.
@@ -63,7 +63,7 @@ centroid(geom) = centroid(GI.trait(geom), geom)
     centroid(
         trait::Union{GI.LineStringTrait, GI.LinearRingTrait},
         geom,
-    )
+    )::Tuple{T, T}
 
 Returns the centroid of a line string or linear ring, which is calculated by
 weighting line segments by their length by convention.
@@ -74,7 +74,7 @@ centroid(
 ) = centroid_and_length(trait, geom)[1]
 
 """
-    centroid(trait, geom)
+    centroid(trait, geom)::Tuple{T, T}
 
 Returns the centroid of a polygon or multipolygon, which is calculated by
 weighting edges by their `area component` by convention.
@@ -82,7 +82,7 @@ weighting edges by their `area component` by convention.
 centroid(trait, geom) = centroid_and_area(trait, geom)[1]
 
 """
-    centroid_and_length(geom)::(GI.Point, ::Real)
+    centroid_and_length(geom)::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and length of a given line/ring. Note this is only valid
 for line strings and linear rings.
@@ -93,14 +93,14 @@ centroid_and_length(geom) = centroid_and_length(GI.trait(geom), geom)
     centroid_and_area(
         ::Union{GI.LineStringTrait, GI.LinearRingTrait}, 
         geom,
-    )::(GI.Point, ::Real)
+    )::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and area of a given geom.
 """
 centroid_and_area(geom) = centroid_and_area(GI.trait(geom), geom)
 
 """
-    centroid_and_length(geom)::(GI.Point, ::Real)
+    centroid_and_length(geom)::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and length of a given line/ring. Note this is only valid
 for line strings and linear rings.
@@ -140,7 +140,7 @@ end
     centroid_and_area(
         ::Union{GI.LineStringTrait, GI.LinearRingTrait},
         geom,
-    )::(GI.Point, ::Real)
+    )::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and area of a given a line string or a linear ring.
 Note that this is only valid if the line segment or linear ring is closed. 
@@ -179,7 +179,7 @@ function centroid_and_area(
 end
 
 """
-    centroid_and_area(::GI.PolygonTrait, geom)::(GI.Point, ::Real)
+    centroid_and_area(::GI.PolygonTrait, geom)::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and area of a given polygon.
 """
@@ -205,7 +205,7 @@ function centroid_and_area(::GI.PolygonTrait, geom)
 end
 
 """
-    centroid_and_area(::GI.MultiPolygonTrait, geom)::(GI.Point, ::Real)
+    centroid_and_area(::GI.MultiPolygonTrait, geom)::(::Tuple{T, T}, ::Real)
 
 Returns the centroid and area of a given multipolygon.
 """