Hi @CosmoProf, you have two approaches to your problem:
Approach 1 (high-level API)
If you are doing interpolation in a geospatial application, then it is a good idea to use the data types in the GeoStats.jl stack to create a high-level problem with all the necessary elements and then call the Kriging solver to solve the problem.
First step consists of georeferencing your data:
using GeoStats
geodata = georef((y=[1.23,5.35,8.64,9.23,7.76],), [1.0 2.0 3.0 4.0 5.0])
5 PointSet{Float64,1}
variables
└─y (Float64)
The second step consists of creating the domain where you want to do the interpolation:
geodomain = RegularGrid((1.0,), (5.0,), dims=(100,))
100 RegularGrid{Float64,1}
origin: (1.0,)
spacing: (0.04040404040404041,)
Finally, you define your problem by specifying which variables you want to interpolate:
problem = EstimationProblem(geodata, geodomain, :y)
1D EstimationProblem
data: 5 SpatialData{Float64,1}
domain: 100 RegularGrid{Float64,1}
variables: y (Float64)
You can solve the problem with the Kriging solver, and specify parameters for each variable:
solver = Kriging(:y => (variogram=GaussianVariogram(),))
solution = solve(problem, solver)
1D EstimationSolution
domain: 100 RegularGrid{Float64,1}
variables: y
You can plot the solution and the data to make sure the result meets yours expectations:
using Plots
plot(solution)
plot!(geodata)
Notice that this approach is particularly useful when you have interpolation in 2D and 3D domains. In your case, the second approach may be simpler.
Approach 2 (low-level API)
If all you need is access to a Kriging interpolator on a set of points, you can use the traditional fit/predict workflow in statistics. You don’t need to load the entire GeoStats.jl stack for that.
First you create your interpolator object:
using KrigingEstimators
using Variography
OK = OrdinaryKriging(GaussianVariogram())
OrdinaryKriging{GaussianVariogram{Float64,Distances.Euclidean}}(GaussianVariogram{Float64,Distances.Euclidean}
range: Float64 1.0
sill: Float64 1.0
nugget: Float64 0.0
distance: Distances.Euclidean
)
Second, you fit it to the data:
krig = fit(OK, [1.0 2.0 3.0 4.0 5.0], [1.23,5.35,8.64,9.23,7.76])
And finally you predict at new locations:
μ, σ² = predict(krig, [1.0])
(1.23, 0.0)
The second result is the conditional variance, which should be zero at interpolation points.
