diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
index 99a19e5..0d25df0 100644
--- a/.github/ISSUE_TEMPLATE/bug_report.md
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -13,7 +13,7 @@ A clear and concise description of what the bug is.
 A [minimal reproducible example](https://stackoverflow.com/help/minimal-reproducible-example) that is enough to show what the problem is
 
 ~~~ julia
-using SimpleSDMLayers
+using SpatialBoundaries
 
 # Add your code here
 ~~~
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
index 2ffcbbb..9691f60 100644
--- a/.github/ISSUE_TEMPLATE/feature_request.md
+++ b/.github/ISSUE_TEMPLATE/feature_request.md
@@ -15,7 +15,7 @@ A clear and concise description of what you want to happen.
 Ideally, this can take the form of code you would like to write:
 
 ~~~ julia
-using SimpleSDMLayers
+using SpatialBoundaries
 
 # Write your dream code here
 ~~~
diff --git a/Project.toml b/Project.toml
index 2b0eb98..717dfa7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SpatialBoundaries"
 uuid = "8d2ba62a-3d23-4a2b-b692-6b63e9267be2"
 authors = ["Tanya Strydom <tanya.strydom@umontreal.ca>", "Timothée Poisot <timothee.poisot@umontreal.ca>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -10,6 +10,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986"
 
 [compat]
+Statistics = "1"
 StatsBase = "0.34"
 VoronoiDelaunay = "0.4"
 julia = "1.9"
diff --git a/src/lib/overallmean.jl b/src/lib/overallmean.jl
index c963fc6..21b670d 100644
--- a/src/lib/overallmean.jl
+++ b/src/lib/overallmean.jl
@@ -18,10 +18,10 @@ function Statistics.mean(w::Vector{T}) where {T <: Womble}
         di = filter(!isnan, [deg2rad(womble.θ[_idx]) for womble in w])
         if !isempty(ch)
             m[_idx] = mean(ch)
-            α[_idx] = atan(mean(sin.(di)), mean(cos.(di)))
+            α[_idx] = atan(sum(sin.(di)), sum(cos.(di)))
         end
     end
-    average_direction = rad2deg.(α) .+ 180.0
+    average_direction = rad2deg.(α .- π)
     if w[1] isa LatticeWomble
         return LatticeWomble(m, average_direction, w[1].x, w[1].y)
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index b0ae9a2..a1c2a95 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,12 +1,14 @@
 global anyerrors = false
 
 tests = [
-    "mock test" => "00_allgood.jl",
+    "Mock test" => "00_allgood.jl",
     "Linear gradient" => "01_gradient.jl",
     "Boundaries detection" => "02_boundaries.jl",
     "Triangulation wombling" => "03_delaunay.jl",
     "Direction of change" => "04_direction.jl",
     "Overall mean wombling" => "05_mean.jl",
+    "Correct values are returned" => "06_values.jl",
+    "Overall mean correct values are returned" => "07_overallmean.jl"
 ]
 
 for test in tests
diff --git a/test/units/06_values.jl b/test/units/06_values.jl
new file mode 100644
index 0000000..9f006bc
--- /dev/null
+++ b/test/units/06_values.jl
@@ -0,0 +1,21 @@
+module SBTestGradient
+
+using SpatialBoundaries
+using NeutralLandscapes
+using Test
+
+angles = rand(100) .* 360
+
+for angle in angles
+    landscape = rand(PlanarGradient(angle), 20, 20)
+    x = collect(LinRange(0.2, 1.8, size(landscape, 1)))
+    y = collect(LinRange(0.2, 1.8, size(landscape, 2)))
+
+    W = wombling(x, y, landscape)
+
+    for θ in W.θ
+        @test deg2rad(abs(θ  - angle)) ≈ π atol = 0.01
+    end
+end
+
+end
diff --git a/test/units/07_overallmean.jl b/test/units/07_overallmean.jl
new file mode 100644
index 0000000..73855b7
--- /dev/null
+++ b/test/units/07_overallmean.jl
@@ -0,0 +1,32 @@
+module SBTestGradient
+
+using SpatialBoundaries
+using NeutralLandscapes
+using Test
+
+angles = rand(100) .* 360
+
+for replicate in Base.OneTo(10)
+    θ1 = rand(angles)
+    θ2 = rand(angles)
+    L1 = rand(PlanarGradient(θ1), 20, 20)
+    L2 = rand(PlanarGradient(θ2), size(L1)...)
+    x = collect(LinRange(0.2, 1.8, size(L1, 1)))
+    y = collect(LinRange(0.2, 1.8, size(L1, 2)))
+
+    W = mean([wombling(x, y, L) for L in [L1, L2]])
+
+    α1 = deg2rad(θ1)
+    α2 = deg2rad(θ2)
+    target = atan(sin(α1)+sin(α2), cos(α1)+cos(α2))
+    observed = deg2rad(rand(W.θ))
+
+    diff = (target-observed) - 2π
+    if abs(diff) > 0.01
+        @test diff ≈ -2π atol = 0.001
+    else
+        @test diff ≈ 0.0 atol = 0.001
+    end
+end
+
+end