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

@@ -135,7 +135,7 @@ function fit_impl(
 
   # evaluate weights
   f = algo.weightfun
-  w = isnothing(f) ? n / sum(n) : map(f, x)
+  w = isnothing(f) ? n / sum(n) : map(xᵢ -> f(ustrip(xᵢ)), x)
 
   # objective function
   function J(θ)
@@ -160,15 +160,15 @@ function fit_impl(
   rₒ = isnothing(range) ? rmax / 3 : range
   sₒ = isnothing(sill) ? 0.95 * smax : sill
   nₒ = isnothing(nugget) ? 1e-6 : nugget
-  θₒ = [rₒ, sₒ, nₒ]
+  θₒ = [ustrip(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 = [rₗ, sₗ, nₗ]
-  u = [rᵤ, sᵤ, nᵤ]
+  l = [ustrip(rₗ), sₗ, nₗ]
+  u = [ustrip(rᵤ), sᵤ, nᵤ]
 
   # solve optimization problem
   sol = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ)

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(>(zero(eltype(s))), s))
+  r > zero(r) ? min(r, l) / 3 : l / 3
 end
 
 function _dims(γ::Variogram, g::Geometry)

src/utils.jl

@@ -2,6 +2,22 @@
 # 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(xs::Quantity...) = ustrip.(promote(xs...))
+
+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
 
@@ -139,6 +137,8 @@ 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)
 

src/variogram/circular.jl

@@ -19,10 +19,11 @@ variotype(::CircularVariogram) = CircularVariogram
 
 isstationary(::Type{<:CircularVariogram}) = true
 
-function (γ::CircularVariogram)(h)
+function (γ::CircularVariogram)(h::Len)
   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::Len)
   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

@@ -24,9 +24,10 @@ variotype(::ExponentialVariogram) = ExponentialVariogram
 
 isstationary(::Type{<:ExponentialVariogram}) = true
 
-function (γ::ExponentialVariogram)(h)
+function (γ::ExponentialVariogram)(h::Len)
   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

@@ -23,11 +23,12 @@ variotype(::GaussianVariogram) = GaussianVariogram
 
 isstationary(::Type{<:GaussianVariogram}) = true
 
-function (γ::GaussianVariogram)(h)
+function (γ::GaussianVariogram)(h::Len)
   # add small eps to nugget
   # for numerical stability
   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::Len)
   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′ = 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

src/variogram/nugget.jl

@@ -21,10 +21,10 @@ isisotropic(::NuggetEffect) = true
 
 sill(γ::NuggetEffect) = γ.nugget
 
-Base.range(::NuggetEffect{T}) where {T} = zero(T)
+Base.range(::NuggetEffect) = 0u"m"
 
 scale(γ::NuggetEffect, ::Real) = γ
 
-(γ::NuggetEffect)(h) = (h > 0) * γ.nugget
+(γ::NuggetEffect)(h::Len) = (h > zero(h)) * γ.nugget
 
 (γ::NuggetEffect)(u::Point, v::Point) = ifelse(u == v, zero(γ.nugget), γ.nugget)