From 868e9fad2d971bfd66c1f511db6ebd272fd0f7dd Mon Sep 17 00:00:00 2001 From: Skylar A Gering Date: Fri, 29 Sep 2023 14:49:27 -0700 Subject: [PATCH] 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 measure of randomness for difference in radius between points (between 0 and 1) 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