diff --git a/Project.toml b/Project.toml
index 795a4fb..9dccec5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "LocalAnisotropies"
 uuid = "cb698b3a-2f1d-48b4-a1e9-ddc0b8f648e2"
 authors = ["Rafael Caixeta and contributors"]
-version = "0.5.4"
+version = "0.5.5"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -24,6 +24,7 @@ SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
@@ -35,16 +36,16 @@ LocalAnisotropiesMakieExt = "Makie"
 [compat]
 CSV = "0.10"
 Distances = "0.10"
-GeoStatsBase = "0.43"
-GeoStatsFunctions = "0.1, 0.2"
-GeoStatsModels = "0.2, 0.3"
-GeoStatsTransforms = "0.2, 0.3, 0.4, 0.5, 0.6"
+GeoStatsBase = "0.44"
+GeoStatsFunctions = "0.4"
+GeoStatsModels = "0.4"
+GeoStatsTransforms = "0.7"
 GeoTables = "1.20"
 Graphs = "1.4"
 ImageFiltering = "0.7"
 LinearAlgebra = "<0.0.1, 1"
-Makie = "0.20"
-Meshes = "0.38, 0.39, 0.40, 0.41, 0.42"
+Makie = "0.21"
+Meshes = "0.45"
 MultivariateStats = "0.10"
 NearestNeighbors = "0.4"
 RecipesBase = "1.1"
@@ -54,6 +55,7 @@ SimpleWeightedGraphs = "1.2"
 StaticArrays = "1.2"
 StatsBase = "0.33, 0.34"
 Tables = "1.10"
+Unitful = "1.20"
 WriteVTK = "1.10"
 julia = "1.9"
 
diff --git a/ext/LocalAnisotropiesMakieExt.jl b/ext/LocalAnisotropiesMakieExt.jl
index de66f8d..fdb786d 100644
--- a/ext/LocalAnisotropiesMakieExt.jl
+++ b/ext/LocalAnisotropiesMakieExt.jl
@@ -20,7 +20,7 @@ Makie.plottype(::LocalAnisotropy, ::LocalAnisotropies.SpatialData) = LA{<:Tuple{
 
 function pre_ellipsoid(lpars, geom, f)
   # pre process
-  coords = LocalAnisotropies.coords(geom)
+  coords = LocalAnisotropies.coords_(geom)
   D3 = size(coords,1) >= 3
   x = coords[1,:]
   y = coords[2,:]
@@ -55,7 +55,7 @@ function Makie.plot!(plot::LA)
   s = Makie.@lift $out[4]
   q = Makie.@lift $out[5]
 
-  Makie.meshscatter!(plot, x, y, z, markersize=s, rotations=q,
+  Makie.meshscatter!(plot, x, y, z, markersize=s, rotation=q,
      alpha=plot[:alpha], color=plot[:color], colormap=plot[:colormap])
 end
 
diff --git a/src/LocalAnisotropies.jl b/src/LocalAnisotropies.jl
index 0024691..e480701 100644
--- a/src/LocalAnisotropies.jl
+++ b/src/LocalAnisotropies.jl
@@ -22,6 +22,7 @@ using SimpleWeightedGraphs
 using StaticArrays
 using StatsBase:Weights,quantile,mean!
 using Tables
+using Unitful:ustrip
 using GeoStatsFunctions
 using WriteVTK
 
diff --git a/src/estimators/variograms.jl b/src/estimators/variograms.jl
index b9be42d..69c1577 100644
--- a/src/estimators/variograms.jl
+++ b/src/estimators/variograms.jl
@@ -21,7 +21,7 @@ function mwvario(model, localpar)
   # get reference pars. and apply local anisotropy to given structures
   p = structures(model.γ)
   γs = map(p[3]) do γ
-    Qs = Q ./ radii(γ.ball)' .^ 2
+    Qs = ustrip(Q ./ collect(radii(γ.ball))' .^ 2)
     γ = @set γ.ball = MetricBall(1.0, Mahalanobis(Symmetric(Qs)))
   end
   γl = NuggetEffect(p[1]) + sum(c*γ for (c, γ) in zip(p[2], γs))
@@ -58,7 +58,7 @@ function kccov(γ::Variogram, xi, xj, Qi::AbstractMatrix, Qj::AbstractMatrix)
   ## modify variogram with average anisotropy matrix
   p = structures(γ)
   γs = map(p[3]) do γx
-    Qs = Qij ./ radii(γx.ball)' .^ 2
+    Qs = ustrip(Qij ./ collect(radii(γx.ball))' .^ 2)
     γx = @set γx.ball = MetricBall(1.0, Mahalanobis(Symmetric(Qs)))
   end
   γl = NuggetEffect(p[1]) + sum(c*γx for (c, γx) in zip(p[2], γs))
diff --git a/src/parametrization.jl b/src/parametrization.jl
index af9723f..4c3a120 100644
--- a/src/parametrization.jl
+++ b/src/parametrization.jl
@@ -115,9 +115,9 @@ spars(D::LocalGeoData, i::Int) = spars(D)[i-nvals(D)]
 spars(D::LocalGeoData, i::AbstractVector{Int}) = view(spars(D), i .- nvals(D))
 nall(D::LocalGeoData) = sdata(D) ? nvals(D)+snvals(D) : nvals(D)
 centro(D::LocalGeoData, i::Int) = i<=nvals(D) ? centro(obj(D),i) : centro(sobj(D),i-nvals(D))
-coords(D::LocalGeoData, i::Int) = coordinates(centro(D,i))
-coords(D::SpatialData) = reduce(hcat, [coordinates(centro(D,x)) for x in 1:nvals(D)])
-coords(D::SpatialData, i::AbstractVector{Int}) = reduce(hcat, [coordinates(centro(D,x)) for x in i])
+coords_(D::LocalGeoData, i::Int) = ustrip(to(centro(D,i)))
+coords_(D::SpatialData) = reduce(hcat, [ustrip(to(centro(D,x))) for x in 1:nvals(D)])
+coords_(D::SpatialData, i::AbstractVector{Int}) = reduce(hcat, [ustrip(to(centro(D,x))) for x in i])
 
 rotation(D::LocalGeoData) = D.localaniso.rotation
 srotation(D::LocalGeoData) = view(rotation(D),spars(D))
diff --git a/src/parametrization/geometric.jl b/src/parametrization/geometric.jl
index c6be6f6..b8b10db 100644
--- a/src/parametrization/geometric.jl
+++ b/src/parametrization/geometric.jl
@@ -26,7 +26,7 @@ lpars = idwpars(prelpars, searcher, grid) # interpolate local pars to a grid
 function localanisotropies(::Type{Geometrical}, searcher::NeighborSearchMethod;
 	simplify::Bool=true)
 	D = searcher.domain
-	X = coords(D)
+	X = coords_(D)
 	N, len = size(X)
 
     quat = Array{Quaternion}(undef,len)
@@ -60,7 +60,7 @@ end
 
 function pca(X, simplify)
 	N = size(X,1)
-	M = MultivariateStats.fit(PCA, X, maxoutdim=N, pratio=1)
+	M = MultivariateStats.fit(PCA, ustrip(X), maxoutdim=N, pratio=1)
 	λ = principalvars(M) ./ principalvars(M)[1]
 	nv = length(λ)
 	v = N == 3 ? projection(M) : vcat(projection(M),[0 0 1][1:nv])
diff --git a/src/parametrization/interpolation.jl b/src/parametrization/interpolation.jl
index 9abbcc1..a5b024d 100644
--- a/src/parametrization/interpolation.jl
+++ b/src/parametrization/interpolation.jl
@@ -72,7 +72,7 @@ function interpolate(lpars, searcher::NeighborSearchMethod, domain=nothing;
 
     Threads.@threads for i in 1:len
         ic  = domain==nothing ? centro(D,i) : centro(domain,i)
-		icoords  = coordinates(ic)
+		icoords  = ustrip(to(ic))
 		neighids = search(ic, searcher)
 
 		if length(neighids) == 0
@@ -85,7 +85,7 @@ function interpolate(lpars, searcher::NeighborSearchMethod, domain=nothing;
 				m[:,i] .= mapreduce(x->quantile(view(mi,x,:),0.5), vcat, 1:N)
 			end
 		else
-			xcoords = coords(D, neighids)
+			xcoords = coords_(D, neighids)
 			prewgts = 1 ./ (eps() .+ Distances.colwise(metric, icoords, xcoords)) .^ power
 			weights = prewgts ./ sum(prewgts)
 			quat[i] = quatavg(rotation(lpars,neighids), weights)
diff --git a/src/parametrization/recipes.jl b/src/parametrization/recipes.jl
index 8f22577..e2a3593 100644
--- a/src/parametrization/recipes.jl
+++ b/src/parametrization/recipes.jl
@@ -9,8 +9,8 @@
   quats = lpars.rotation
   u = [quat_to_dcm(quats[x])[1,1] for x in 1:length(quats)]
   v = [quat_to_dcm(quats[x])[1,2] for x in 1:length(quats)]
-  x = [coordinates(centro(D,x))[1] for x in 1:nvals(D)]
-  y = [coordinates(centro(D,x))[2] for x in 1:nvals(D)]
+  x = [to(centro(D,x))[1] for x in 1:nvals(D)]
+  y = [to(centro(D,x))[2] for x in 1:nvals(D)]
   c = lpars.magnitude[iaxis(axis),:]
 
   seriestype --> :quiver
diff --git a/src/parametrization/searchers.jl b/src/parametrization/searchers.jl
index 571e3bc..50693eb 100644
--- a/src/parametrization/searchers.jl
+++ b/src/parametrization/searchers.jl
@@ -13,8 +13,8 @@ end
 
 # nearest grid data id to some other sample data
 function grid2hd_ids(pdata,pdomain)
-  hd   = [coordinates(centro(pdata, i)) for i in 1:nvals(pdata)]
-  grid = [coordinates(centro(pdomain, i)) for i in 1:nvals(pdomain)]
+  hd   = [ustrip(to(centro(pdata, i))) for i in 1:nvals(pdata)]
+  grid = [ustrip(to(centro(pdomain, i))) for i in 1:nvals(pdomain)]
 
   tree = KDTree(grid)
   idxs, dists = nn(tree, hd)
diff --git a/src/parametrization/utilities.jl b/src/parametrization/utilities.jl
index 0fbd1ee..fa57b64 100644
--- a/src/parametrization/utilities.jl
+++ b/src/parametrization/utilities.jl
@@ -112,7 +112,7 @@ function adjust_rake!(lpar::LocalAnisotropy, az::AbstractVector)
 		p2 = Plane(Point(0,0,0),n2)
 
 		dcm = collect(dcm)
-		icross = coordinates(intersection(p1,p2).geom.b)
+		icross = ustrip(to(intersection(p1,p2).geom.b))
 		dcm[1,:] .= icross
 		dcm[2,:] .= cross(dcm[1,:],dcm[3,:])
 		det(dcm) < 0 && (dcm = Diagonal([-1,1,1]) * dcm)
@@ -195,7 +195,7 @@ function to_vtk(vtkfile, coords::AbstractArray, lpars::LocalAnisotropy;
 end
 
 function to_vtk(vtkfile, D::SpatialData, lpars::LocalAnisotropy; kwargs...)
-	coords = reduce(hcat, [coordinates(centro(D,x)) for x in 1:nvals(D)])
+	coords = reduce(hcat, [ustrip(to(centro(D,x))) for x in 1:nvals(D)])
 	to_vtk(vtkfile, coords, lpars; kwargs...)
 end
 
diff --git a/src/parametrization/variograms.jl b/src/parametrization/variograms.jl
index 6ace41c..603b433 100644
--- a/src/parametrization/variograms.jl
+++ b/src/parametrization/variograms.jl
@@ -30,13 +30,13 @@ function pseudolocalpartition(obj, lpars, axis, tol, maxratio1, maxratio2)
   for i in 1:n
     maxratio1 != Inf && ratio1(lpars,i) > maxratio1 && continue
     maxratio2 != Inf && ratio2(lpars,i) > maxratio2 && continue
-    x = coordinates(centro(obj, i))
+    x = to(centro(obj, i))
     v = rotmat(lpars, i)[iaxis(axis),1:dims]
     p = DirectionPartition(Tuple(v), tol=tol)
     s = Int[i]
     for j in 1:n
       i == j && continue
-      y = coordinates(centro(obj, j))
+      y = to(centro(obj, j))
       p(x, y) && push!(s, j)
     end
     push!(subs, s)
diff --git a/src/spacetransforms/deformation.jl b/src/spacetransforms/deformation.jl
index d1145c3..c77d74d 100644
--- a/src/spacetransforms/deformation.jl
+++ b/src/spacetransforms/deformation.jl
@@ -157,6 +157,6 @@ first three dimensions in order to help plotting of the data.
 """
 function to_3d(s)
 	dom = domain(s)
-	c = reduce(hcat,[coordinates(centro(dom,x))[1:3] for x in 1:nvals(dom)])
+	c = reduce(hcat,[to(centro(dom,x))[1:3] for x in 1:nvals(dom)])
 	georef(values(s),PointSet(c))
 end
diff --git a/src/spacetransforms/metrics.jl b/src/spacetransforms/metrics.jl
index db2cda2..c2c440f 100644
--- a/src/spacetransforms/metrics.jl
+++ b/src/spacetransforms/metrics.jl
@@ -3,12 +3,12 @@
 # ------------------------------------------------------------------
 
 function colwise(D::LocalGeoData, ::Type{AnisoDistance}, i::Int, xj::Vector{Int})
-	Qi, ix = qmat(rotation(D,i),magnitude(D,i)), coords(D,i)
+	Qi, ix = qmat(rotation(D,i),magnitude(D,i)), coords_(D,i)
 	d = Vector{Float64}()
 	for j in xj
 		Qj = qmat(rotation(D,j),magnitude(D,j))
 		Qij = (Qi+Qj)/2
-		jx = coords(D,j)
+		jx = coords_(D,j)
 		push!(d,Distances.evaluate(Mahalanobis(Symmetric(Qij)), ix, jx))
 	end
 	d
@@ -33,7 +33,7 @@ end
 
 function evaluate(D::LocalGeoData, ::Type{AnisoDistance}, i::Int, j::Int)
 	Qi, Qj = [qmat(rotation(D,x),magnitude(D,x)) for x in (i,j)]
-	ix, jx = [coords(D,x) for x in (i,j)]
+	ix, jx = [coords_(D,x) for x in (i,j)]
 	Qij = (Qi+Qj)/2
 	Distances.evaluate(Mahalanobis(Symmetric(Qij)), ix, jx)
 end