using Plots, GeoStats, InverseDistanceWeighting, Distances
# data
x = [50.0, 100.0, 200.0]
y = [-30.0, 30.0, 10.0]
z = [4.0, -1.0, 3.0]
# list of properties with coordinates
props = OrderedDict(:prop => z)
coord = [(lon,lat) for (lon,lat) in zip(x,y)]
# define spatial data
sdata = PointSetData(props, coord)
# the solution domain (basically a grid covering the entire globe)
dims = (180, 89)
start = (1.0, -89.0)
finish = (359.0, 89.0)
sdomain = RegularGrid(start, finish, dims=dims)
# define estimation problem for any data column(s) (e.g. :precipitation)
problem = EstimationProblem(sdata, sdomain, :prop)
# axes for plots
xs = range(start[1], stop=finish[1], length=dims[1])
ys = range(start[2], stop=finish[2], length=dims[2])
# Plotting solution with data on top
plts = Any[]
tit(x) = string(x.distance, ", neighbours=", x.neighbors)
for NT in [NamedTuple(),
           (distance=Haversine(1.0),),
           (distance=Chebyshev(),),
           (neighbors=1,)]
    solver = InvDistWeight(:variable => NT)
    solution = GeoStats.solve(problem, solver)
    zs = permutedims(reshape(solution.mean[:prop], size(solution.domain)), (2,1))
    plt = contourf(xs, ys, zs)
    scatter!(plt, sdata, title=tit(solver.vparams[:variable]), titlefontsize=10)
    push!(plts, plt)
end
plot(plts..., layout=(2,2))