Update to Meshes.jl v0.43

.github/workflows/FormatPR.yml

@@ -2,7 +2,7 @@ name: FormatPR
 on:
   push:
     branches:
-      - master
+      - main
 jobs:
   build:
     runs-on: ubuntu-latest

Project.toml

@@ -29,11 +29,11 @@ GeoStatsFunctionsMakieExt = "Makie"
 [compat]
 Bessels = "0.2"
 Distances = "0.10"
-GeoTables = "1.19"
+GeoTables = "1.21"
 InteractiveUtils = "1.9"
 LinearAlgebra = "1.9"
 Makie = "0.20"
-Meshes = "0.42"
+Meshes = "0.43"
 NearestNeighbors = "0.4"
 Optim = "1.7"
 Printf = "1.9"

src/algorithms/ballsearch.jl

@@ -3,17 +3,20 @@
 # ------------------------------------------------------------------
 
 """
-    BallSearchAccum(maxlag, nlags, distance)
+    BallSearchAccum(nlags, maxlag, distance)
 
 Accumulate pairs of points in geospatial data with
 nearest neighbors inside metric ball.
 """
-struct BallSearchAccum{T,D} <: VariogramAccumAlgo
+struct BallSearchAccum{ℒ<:Len,D} <: VariogramAccumAlgo
   nlags::Int
-  maxlag::T
+  maxlag::ℒ
   distance::D
+  BallSearchAccum(nlags, maxlag::ℒ, distance::D) where {ℒ<:Len,D} = new{float(ℒ),D}(nlags, maxlag, distance)
 end
 
+BallSearchAccum(nlags, maxlag, distance) = BallSearchAccum(nlags, addunit(maxlag, u"m"), distance)
+
 function neighfun(algo::BallSearchAccum, pset)
   ball = MetricBall(algo.maxlag, algo.distance)
   searcher = BallSearch(pset, ball)

src/algorithms/fullsearch.jl

@@ -3,17 +3,20 @@
 # ------------------------------------------------------------------
 
 """
-    FullSearchAccum(maxlag, nlags, distance)
+    FullSearchAccum(nlags, maxlag, distance)
 
 Accumulate pairs of points in geospatial data with
 exhaustive (or full) search.
 """
-struct FullSearchAccum{T,D} <: VariogramAccumAlgo
+struct FullSearchAccum{ℒ<:Len,D} <: VariogramAccumAlgo
   nlags::Int
-  maxlag::T
+  maxlag::ℒ
   distance::D
+  FullSearchAccum(nlags, maxlag::ℒ, distance::D) where {ℒ<:Len,D} = new{float(ℒ),D}(nlags, maxlag, distance)
 end
 
+FullSearchAccum(nlags, maxlag, distance) = FullSearchAccum(nlags, addunit(maxlag, u"m"), distance)
+
 neighfun(::FullSearchAccum, pset) = j -> (j + 1):nelements(pset)
 
 skipfun(::FullSearchAccum) = (i, j) -> false

src/empirical.jl

@@ -60,7 +60,8 @@ function accumulate(data, var₁, var₂, estim::VariogramEstimator, algo::Vario
   V = returntype(estim, z₁, z₂)
 
   # lag sums and counts
-  Σx = zeros(nlags)
+  ℒ = Meshes.lentype(𝒫)
+  Σx = zeros(ℒ, nlags)
   Σy = zeros(V, nlags)
   ns = zeros(Int, nlags)
 
@@ -78,7 +79,7 @@ function accumulate(data, var₁, var₂, estim::VariogramEstimator, algo::Vario
       z₂ᵢ = z₂[i]
 
       # evaluate geospatial lag
-      h = evaluate(distance, coordinates(pᵢ), coordinates(pⱼ))
+      h = uevaluate(distance, pᵢ, pⱼ)
 
       # early exit if out of range
       exit(h) && continue

src/empirical/partition.jl

@@ -26,29 +26,29 @@ function EmpiricalVariogram(partition::Partition, var₁, var₂=var₁; kwargs.
 end
 
 """
-    DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6, [parameters])
+    DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6u"m", [parameters])
 
 Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in
-geospatial `data` along a given `direction` with band tolerance `dtol`.
+geospatial `data` along a given `direction` with band tolerance `dtol` in length units.
 
 Optionally, forward `parameters` for the underlying [`EmpiricalVariogram`](@ref).
 """
-function DirectionalVariogram(dir, data::AbstractGeoTable, var₁, var₂=var₁; dtol=1e-6, kwargs...)
+function DirectionalVariogram(dir, data::AbstractGeoTable, var₁, var₂=var₁; dtol=1e-6u"m", kwargs...)
   rng = MersenneTwister(123)
   Π = partition(rng, data, DirectionPartition(dir; tol=dtol))
   EmpiricalVariogram(Π, var₁, var₂; kwargs...)
 end
 
 """
-    PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6, [parameters])
+    PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6u"m", [parameters])
 
 Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in
 geospatial `data` along a plane perpendicular to a `normal` direction with plane
-tolerance `ntol`.
+tolerance `ntol` in length units.
 
 Optionally, forward `parameters` for the underlying [`EmpiricalVariogram`](@ref).
 """
-function PlanarVariogram(normal, data::AbstractGeoTable, var₁, var₂=var₁; ntol=1e-6, kwargs...)
+function PlanarVariogram(normal, data::AbstractGeoTable, var₁, var₂=var₁; ntol=1e-6u"m", kwargs...)
   rng = MersenneTwister(123)
   Π = partition(rng, data, PlanePartition(normal; tol=ntol))
   EmpiricalVariogram(Π, var₁, var₂; kwargs...)

src/empirical/variogram.jl

@@ -12,7 +12,7 @@ geospatial `data`.
 ## Parameters
 
   * nlags     - number of lags (default to `20`)
-  * maxlag    - maximum lag (default to 1/10 of maximum lag of data)
+  * maxlag    - maximum lag in length units (default to 1/10 of maximum lag of data)
   * distance  - custom distance function (default to `Euclidean` distance)
   * estimator - variogram estimator (default to `:matheron` estimator)
   * algorithm - accumulation algorithm (default to `:ball`)
@@ -44,8 +44,8 @@ See also: [`DirectionalVariogram`](@ref), [`PlanarVariogram`](@ref),
 * Hoffimann, J and Zadrozny, B. 2019. [Efficient variography with partition variograms]
   (https://www.sciencedirect.com/science/article/pii/S0098300419302936)
 """
-struct EmpiricalVariogram{V,D,E}
-  abscissa::Vector{Float64}
+struct EmpiricalVariogram{ℒ<:Len,V,D,E}
+  abscissa::Vector{ℒ}
   ordinate::Vector{V}
   counts::Vector{Int}
   distance::D
@@ -73,7 +73,7 @@ function EmpiricalVariogram(
   # sanity checks
   @assert nelem > 1 "variogram requires at least 2 elements"
   @assert nlags > 0 "number of lags must be positive"
-  @assert maxlag > 0 "maximum lag distance must be positive"
+  @assert maxlag > zero(maxlag) "maximum lag distance must be positive"
   @assert estimator ∈ (:matheron, :cressie) "invalid variogram estimator"
   @assert algorithm ∈ (:full, :ball) "invalid accumulation algorithm"
 
@@ -82,7 +82,7 @@ function EmpiricalVariogram(
 
   # ball search with NearestNeighbors.jl requires AbstractFloat and MinkowskiMetric
   # https://github.com/KristofferC/NearestNeighbors.jl/issues/13
-  isfloat = coordtype(𝒟) <: AbstractFloat
+  isfloat = Unitful.numtype(Meshes.lentype(𝒟)) <: AbstractFloat
   isminkowski = distance isa MinkowskiMetric
 
   # warn users requesting :ball option with invalid parameters

src/empirical/varioplane.jl

@@ -5,14 +5,14 @@
 """
     EmpiricalVarioplane(data, var₁, var₂=var₁;
                         normal=spheredir(0,0),
-                        nangs=50, ptol=0.5, dtol=0.5,
+                        nangs=50, ptol=0.5u"m", dtol=0.5u"m",
                         [parameters])
 
 Given a `normal` direction, estimate the (cross-)variogram of variables
 `var₁` and `var₂` along all directions in the corresponding plane of variation.
 
-Optionally, specify the tolerance `ptol` for the plane partition, the tolerance
-`dtol` for the direction partition, the number of angles `nangs` in the plane,
+Optionally, specify the tolerance `ptol` in length units for the plane partition, the tolerance
+`dtol` in length units for the direction partition, the number of angles `nangs` in the plane,
 and forward the `parameters` to the underlying [`EmpiricalVariogram`](@ref).
 """
 struct EmpiricalVarioplane{T,V}
@@ -26,8 +26,8 @@ function EmpiricalVarioplane(
   var₂=var₁;
   normal=spheredir(0, 0),
   nangs=50,
-  ptol=0.5,
-  dtol=0.5,
+  ptol=0.5u"m",
+  dtol=0.5u"m",
   kwargs...
 )
   # sanity checks

src/fitting.jl

@@ -129,44 +129,54 @@ function fit_impl(
   y = y[n .> 0]
   n = n[n .> 0]
 
-  # strip units if necessary
-  𝓊 = unit(first(y))
-  y = ustrip.(y)
+  # strip units of coordinates
+  ux = unit(eltype(x))
+  uy = unit(eltype(y))
+  x′ = ustrip.(x)
+  y′ = ustrip.(y)
+
+  # strip units of kwargs
+  range′ = isnothing(range) ? range : ustrip(ux, range)
+  sill′ = isnothing(sill) ? sill : ustrip(uy, sill)
+  nugget′ = isnothing(nugget) ? nugget : ustrip(uy, nugget)
+  maxrange′ = isnothing(maxrange) ? maxrange : ustrip(ux, maxrange)
+  maxsill′ = isnothing(maxsill) ? maxsill : ustrip(uy, maxsill)
+  maxnugget′ = isnothing(maxnugget) ? maxnugget : ustrip(uy, maxnugget)
 
   # evaluate weights
   f = algo.weightfun
-  w = isnothing(f) ? n / sum(n) : map(f, x)
+  w = isnothing(f) ? n / sum(n) : map(xᵢ -> ustrip(f(xᵢ)), x)
 
   # objective function
   function J(θ)
     γ = V(ball(θ[1]), sill=θ[2], nugget=θ[3])
-    sum(i -> w[i] * (γ(x[i]) - y[i])^2, eachindex(w, x, y))
+    sum(i -> w[i] * (γ(x′[i]) - y′[i])^2, eachindex(w, x′, y′))
   end
 
   # linear constraint (sill ≥ nugget)
   L(θ) = θ[2] ≥ θ[3] ? 0.0 : θ[3] - θ[2]
 
   # penalty for linear constraint (J + λL)
-  λ = sum(yᵢ -> yᵢ^2, y)
+  λ = sum(yᵢ -> yᵢ^2, y′)
 
   # maximum range, sill and nugget
-  xmax = maximum(x)
-  ymax = maximum(y)
-  rmax = isnothing(maxrange) ? xmax : maxrange
-  smax = isnothing(maxsill) ? ymax : maxsill
-  nmax = isnothing(maxnugget) ? ymax : maxnugget
+  xmax = maximum(x′)
+  ymax = maximum(y′)
+  rmax = isnothing(maxrange′) ? xmax : maxrange′
+  smax = isnothing(maxsill′) ? ymax : maxsill′
+  nmax = isnothing(maxnugget′) ? ymax : maxnugget′
 
   # initial guess
-  rₒ = isnothing(range) ? rmax / 3 : range
-  sₒ = isnothing(sill) ? 0.95 * smax : sill
-  nₒ = isnothing(nugget) ? 1e-6 : nugget
+  rₒ = isnothing(range′) ? rmax / 3 : range′
+  sₒ = isnothing(sill′) ? 0.95 * smax : sill′
+  nₒ = isnothing(nugget′) ? 1e-6 : nugget′
   θₒ = [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 + δ)
+  rₗ, rᵤ = isnothing(range′) ? (zero(rmax), rmax) : (range′ - δ, range′ + δ)
+  sₗ, sᵤ = isnothing(sill′) ? (zero(smax), smax) : (sill′ - δ, sill′ + δ)
+  nₗ, nᵤ = isnothing(nugget′) ? (zero(nmax), nmax) : (nugget′ - δ, nugget′ + δ)
   l = [rₗ, sₗ, nₗ]
   u = [rᵤ, sᵤ, nᵤ]
 
@@ -176,7 +186,7 @@ function fit_impl(
   θ = Optim.minimizer(sol)
 
   # optimal variogram (with units)
-  γ = V(ball(θ[1]), sill=θ[2] * 𝓊, nugget=θ[3] * 𝓊)
+  γ = V(ball(θ[1] * ux), sill=θ[2] * uy, nugget=θ[3] * uy)
 
   γ, ϵ
 end

src/sampling.jl

@@ -40,8 +40,8 @@ end
 function _spacing(γ::Variogram, g::Geometry)
   r = range(γ)
   s = sides(boundingbox(g))
-  l = minimum(filter(>(0), s))
-  r > 0 ? min(r, l) / 3 : l / 3
+  l = minimum(filter(sᵢ -> sᵢ > zero(sᵢ), s))
+  r > zero(r) ? min(r, l) / 3 : l / 3
 end
 
 function _dims(γ::Variogram, g::Geometry)

src/utils.jl

@@ -2,6 +2,26 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
+const Len{T} = Quantity{T,u"𝐋"}
+
+addunit(x::Number, u) = x * u
+addunit(::Quantity, _) = throw(ArgumentError("invalid units, please check the documentation"))
+
+unitless(a::Number, b::Quantity) = a, ustrip(b)
+function unitless(a::Quantity, b::Quantity)
+  u = Unitful.promote_unit(unit(a), unit(b))
+  ustrip(u, a), ustrip(u, b)
+end
+
+function uevaluate(distance, p₁, p₂)
+  u₁ = unit(Meshes.lentype(p₁))
+  u₂ = unit(Meshes.lentype(p₂))
+  u = Unitful.promote_unit(u₁, u₂)
+  v₁ = ustrip.(u, to(p₁))
+  v₂ = ustrip.(u, to(p₂))
+  evaluate(distance, v₁, v₂) * u
+end
+
 """
     spheredir(θ, φ)
 

src/variogram.jl

@@ -105,9 +105,7 @@ Evaluate the variogram at points `u` and `v`.
 """
 function (γ::Variogram)(u::Point, v::Point)
   d = metric(γ.ball)
-  x = coordinates(u)
-  y = coordinates(v)
-  h = evaluate(d, x, y)
+  h = uevaluate(d, u, v)
   γ(h)
 end
 

src/variogram/circular.jl

@@ -23,6 +23,7 @@ function (γ::CircularVariogram)(h)
   r = radius(γ.ball)
   s = γ.sill
   n = γ.nugget
-  v = h ≤ r ? 1 - (2 / π) * acos(h / r) + (2h / (π * r)) * sqrt(1 - (h^2 / r^2)) : one(h)
-  (s - n) * v + (h > zero(h)) * n
+  h′, r′ = unitless(h, r)
+  v = h′ ≤ r′ ? 1 - (2 / π) * acos(h′ / r′) + (2h′ / (π * r′)) * sqrt(1 - (h′^2 / r′^2)) : one(h′)
+  (s - n) * v + (h′ > zero(h′)) * n
 end

src/variogram/cubic.jl

@@ -23,17 +23,19 @@ variotype(::CubicVariogram) = CubicVariogram
 
 isstationary(::Type{<:CubicVariogram}) = true
 
-function (γ::CubicVariogram)(h::T) where {T}
+function (γ::CubicVariogram)(h)
   r = radius(γ.ball)
   s = γ.sill
   n = γ.nugget
+  h′, r′ = unitless(h, r)
 
   # constants
+  T = typeof(h′)
   c1 = T(35) / T(4)
   c2 = T(7) / T(2)
   c3 = T(3) / T(4)
-  s1 = 7 * (h / r)^2 - c1 * (h / r)^3 + c2 * (h / r)^5 - c3 * (h / r)^7
+  s1 = 7 * (h′ / r′)^2 - c1 * (h′ / r′)^3 + c2 * (h′ / r′)^5 - c3 * (h′ / r′)^7
   s2 = T(1)
 
-  (h < r) * (s - n) * s1 + (h ≥ r) * (s - n) * s2 + (h > 0) * n
+  (h′ < r′) * (s - n) * s1 + (h′ ≥ r′) * (s - n) * s2 + (h′ > 0) * n
 end

src/variogram/exponential.jl

@@ -28,5 +28,6 @@ function (γ::ExponentialVariogram)(h)
   r = radius(γ.ball)
   s = γ.sill
   n = γ.nugget
-  (s - n) * (1 - exp(-3(h / r))) + (h > 0) * n
+  h′, r′ = unitless(h, r)
+  (s - n) * (1 - exp(-3(h′ / r′))) + (h′ > 0) * n
 end

src/variogram/gaussian.jl

@@ -29,5 +29,6 @@ function (γ::GaussianVariogram)(h)
   r = radius(γ.ball)
   s = γ.sill
   n = γ.nugget + typeof(s)(1e-6)
-  (s - n) * (1 - exp(-3(h / r)^2)) + (h > 0) * n
+  h′, r′ = unitless(h, r)
+  (s - n) * (1 - exp(-3(h′ / r′)^2)) + (h′ > 0) * n
 end

src/variogram/matern.jl

@@ -26,19 +26,20 @@ variotype(::MaternVariogram) = MaternVariogram
 
 isstationary(::Type{<:MaternVariogram}) = true
 
-function (γ::MaternVariogram)(h::T) where {T}
+function (γ::MaternVariogram)(h)
   r = radius(γ.ball)
   s = γ.sill
   n = γ.nugget
   ν = γ.order
+  h′, r′ = unitless(h, r)
 
   # shift lag by machine precision to
   # avoid explosion at the origin
-  h′ = h + eps(T)
+  h′ += eps(typeof(h′))
 
-  δ = √(2ν) * 3(h′ / r)
+  δ = √(2ν) * 3(h′ / r′)
   Β = besselk(ν, δ)
   Γ = gamma(ν)
 
-  (s - n) * (1 - 2^(1 - ν) / Γ * δ^ν * Β) + (h > 0) * n
+  (s - n) * (1 - 2^(1 - ν) / Γ * δ^ν * Β) + (h′ > 0) * n
 end