diff --git a/src/fitting.jl b/src/fitting.jl
index 498a897..6cb304f 100644
--- a/src/fitting.jl
+++ b/src/fitting.jl
@@ -130,12 +130,13 @@ function fit_impl(
   n = n[n .> 0]
 
   # strip units if necessary
-  𝓊 = unit(first(y))
+  ux = unit(first(x))
+  uy = unit(first(y))
   y = ustrip.(y)
 
   # evaluate weights
   f = algo.weightfun
-  w = isnothing(f) ? n / sum(n) : map(xᵢ -> f(ustrip(xᵢ)), x)
+  w = isnothing(f) ? n / sum(n) : map(xᵢ -> ustrip(f(xᵢ)), x)
 
   # objective function
   function J(θ)
@@ -160,15 +161,15 @@ function fit_impl(
   rₒ = isnothing(range) ? rmax / 3 : range
   sₒ = isnothing(sill) ? 0.95 * smax : sill
   nₒ = isnothing(nugget) ? 1e-6 : nugget
-  θₒ = [ustrip(rₒ), sₒ, nₒ]
+  θₒ = [ustrip(ux, rₒ), sₒ, nₒ]
 
   # box constraints
   δ = 1e-8
-  rₗ, rᵤ = isnothing(range) ? (0.0, rmax) : (range - δ, range + δ)
-  sₗ, sᵤ = isnothing(sill) ? (0.0, smax) : (sill - δ, sill + δ)
-  nₗ, nᵤ = isnothing(nugget) ? (0.0, nmax) : (nugget - δ, nugget + δ)
-  l = [ustrip(rₗ), sₗ, nₗ]
-  u = [ustrip(rᵤ), sᵤ, nᵤ]
+  rₗ, rᵤ = isnothing(range) ? (zero(rmax), rmax) : (range - δ * unit(range), range + δ * unit(range))
+  sₗ, sᵤ = isnothing(sill) ? (zero(smax), smax) : (sill - δ, sill + δ)
+  nₗ, nᵤ = isnothing(nugget) ? (zero(nmax), nmax) : (nugget - δ, nugget + δ)
+  l = [ustrip(ux, rₗ), sₗ, nₗ]
+  u = [ustrip(ux, rᵤ), sᵤ, nᵤ]
 
   # solve optimization problem
   sol = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ)
@@ -176,7 +177,7 @@ function fit_impl(
   θ = Optim.minimizer(sol)
 
   # optimal variogram (with units)
-  γ = V(ball(θ[1]), sill=θ[2] * 𝓊, nugget=θ[3] * 𝓊)
+  γ = V(ball(θ[1]), sill=θ[2] * uy, nugget=θ[3] * uy)
 
   γ, ϵ
 end
diff --git a/src/sampling.jl b/src/sampling.jl
index 8df7ee8..b4e42ac 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -40,7 +40,7 @@ end
 function _spacing(γ::Variogram, g::Geometry)
   r = range(γ)
   s = sides(boundingbox(g))
-  l = minimum(filter(>(zero(eltype(s))), s))
+  l = minimum(filter(sᵢ -> sᵢ > zero(sᵢ), s))
   r > zero(r) ? min(r, l) / 3 : l / 3
 end
 
diff --git a/src/variogram.jl b/src/variogram.jl
index 8b5d36f..ce3501c 100644
--- a/src/variogram.jl
+++ b/src/variogram.jl
@@ -137,8 +137,6 @@ function (γ::Variogram)(U::Geometry, V::Geometry)
   mean(γ(u, v) for u in us, v in vs)
 end
 
-(γ::Variogram)(h) = γ(addunit(h, u"m"))
-
 """
     pairwise(γ, domain)
 
diff --git a/test/empirical.jl b/test/empirical.jl
index bbd0164..fdd47f8 100644
--- a/test/empirical.jl
+++ b/test/empirical.jl
@@ -101,8 +101,20 @@
     data = georef((z=img,))
     γ = EmpiricalVarioplane(data, :z, maxlag=50.0)
     @test sprint(show, γ) == "EmpiricalVarioplane"
-    @test sprint(show, MIME"text/plain"(), γ) ==
-          "EmpiricalVarioplane\n  N° pairs\n  └─0.00° → 372500\n  └─3.67° → 304782\n  └─7.35° → 298306\n  └─11.02° → 297432\n  └─14.69° → 297243\n  ⋮\n  └─165.31° → 293643\n  └─168.98° → 295850\n  └─172.65° → 296931\n  └─176.33° → 306528\n  └─180.00° → 372500"
+    @test sprint(show, MIME"text/plain"(), γ) == """
+    EmpiricalVarioplane
+      N° pairs
+      └─0.00° → 372500
+      └─3.67° → 304782
+      └─7.35° → 298306
+      └─11.02° → 297432
+      └─14.69° → 297243
+      ⋮
+      └─165.31° → 293643
+      └─168.98° → 295850
+      └─172.65° → 296931
+      └─176.33° → 306528
+      └─180.00° → 372500"""
   end
 
   @testset "Directional" begin
diff --git a/test/fitting.jl b/test/fitting.jl
index a8d2ba7..878b9ed 100644
--- a/test/fitting.jl
+++ b/test/fitting.jl
@@ -14,28 +14,28 @@
   @test isapprox(sill(γ₄), 0.054, atol=1e-3)
 
   # fix parameters
-  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0)
+  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m")
   @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m")
   γ = GeoStatsFunctions.fit(GaussianVariogram, g, sill=0.07)
   @test isapprox(sill(γ), 0.07, atol=1e-3)
   γ = GeoStatsFunctions.fit(GaussianVariogram, g, nugget=0.05)
   @test isapprox(nugget(γ), 0.05, atol=1e-3)
-  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, sill=0.07)
+  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", sill=0.07)
   @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m")
   @test isapprox(sill(γ), 0.07, atol=1e-3)
-  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, nugget=0.05)
+  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", nugget=0.05)
   @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m")
   @test isapprox(nugget(γ), 0.05, atol=1e-3)
   γ = GeoStatsFunctions.fit(GaussianVariogram, g, sill=0.07, nugget=0.05)
   @test isapprox(sill(γ), 0.07, atol=1e-3)
   @test isapprox(nugget(γ), 0.05, atol=1e-3)
-  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, sill=0.07, nugget=0.05)
+  γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", sill=0.07, nugget=0.05)
   @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m")
   @test isapprox(sill(γ), 0.07, atol=1e-3)
   @test isapprox(nugget(γ), 0.05, atol=1e-3)
 
   # fix maximum parameters
-  γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxrange=5.0)
+  γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxrange=5.0u"m")
   @test isapprox(range(γ), 5.0u"m", atol=1e-3u"m")
   γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxsill=0.04)
   @test isapprox(sill(γ), 0.04, atol=1e-3)
diff --git a/test/nesting.jl b/test/nesting.jl
index e2d0922..8a59b4c 100644
--- a/test/nesting.jl
+++ b/test/nesting.jl
@@ -53,7 +53,7 @@
   γ = C₁ * GaussianVariogram(range=1.0) + C₂ * SphericalVariogram(range=2.0)
   @test range(γ) ≈ 2.0u"m"
   @test sill(γ) ≈ C₁ .+ C₂
-  @test γ(10.0) ≈ sill(γ)
+  @test γ(10.0u"m") ≈ sill(γ)
   @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ)
   @test isstationary(γ)
 
@@ -62,7 +62,7 @@
   γ = C * GaussianVariogram() + C * ExponentialVariogram() + C * CubicVariogram()
   @test range(γ) ≈ 1.0u"m"
   @test sill(γ) ≈ [3.0 0.0; 0.0 3.0]
-  @test γ(10.0) ≈ sill(γ)
+  @test γ(10.0u"m") ≈ sill(γ)
   @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ)
   @test isstationary(γ)
 
diff --git a/test/variogram.jl b/test/variogram.jl
index bdce539..a907aa9 100644
--- a/test/variogram.jl
+++ b/test/variogram.jl
@@ -1,6 +1,6 @@
 @testset "Variogram" begin
   rng = StableRNG(123)
-  h = range(0, stop=10, length=50)
+  h = range(0, stop=10, length=50) * u"m"
   x, y = rand(rng, Point{3}), rand(rng, Point{3})
 
   # stationary variogram models
@@ -61,12 +61,12 @@
 
   # some variograms are non-decreasing
   for γ in (γnd ∪ [sum(γnd)])
-    @test all(γ.(h) .≤ γ.(h .+ 1))
+    @test all(γ.(h) .≤ γ.(h .+ 1u"m"))
   end
 
   # variograms are valid at the origin
   for γ in (γs ∪ γn ∪ γnd)
-    @test !isnan(γ(0.0)) && !isinf(γ(0.0))
+    @test !isnan(γ(0.0u"m")) && !isinf(γ(0.0u"m"))
   end
 
   # practical ranges
@@ -97,8 +97,8 @@
   # nugget regularization
   γ = GaussianVariogram(range=20.0, nugget=0.1)
   C = sill(γ) .- GeoStatsFunctions.pairwise(γ, pset)
-  @test γ(0) == 0
-  @test γ(1e-6) > 0
+  @test γ(0u"m") == 0
+  @test γ(1e-6u"m") > 0
   @test cond(C) < 100.0
 
   # sill and nugget in single precision
@@ -125,7 +125,7 @@
   # unitful non-stationary types
   γn = [PowerVariogram(scaling=1.0u"K^2")]
   for γ in γs
-    @test unit(γ(1.0)) == u"K^2"
+    @test unit(γ(1.0u"m")) == u"K^2"
   end
 
   𝒟 = PointSet(Matrix(1.0I, 3, 3))
@@ -188,9 +188,6 @@
   g = GeoStatsFunctions.scale(γ, 2)
   @test g == γ
 
-  # error: h must have length units
-  @test_throws ArgumentError γ(0.0u"s")
-
   # shows
   γ = CircularVariogram()
   @test sprint(show, γ) == "CircularVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)"