Add more tests and refine

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

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
 

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,

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

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