Confusion on Creating Gridded Interpolations

No worries! Note that, as they say on the ReadMe,

The method was initially developped with ocean data in mind, but it can be applied to any field where localized observations have to be used to produce gridded fields which are “smooth”.

1 Like

It is unclear what structure your grid has (regular or scattered). If scattered, you need a really large basis to interpolate 250k points — unless you are sure you need that many degrees of freedom, perhaps you are better served with eg a least squares approximation using a relatively low-rank basis.

Let me add that these issues have been carefully addressed in GeoStats.jl, and that you can not only do smooth interpolations, but also simulate spatial processes in a grid without effort.

Notice that most interpolation packages assume that you can store the full grid of coordinates into memory, whereas in GeoStats.jl we never materialize the grid:

julia> using GeoStats
julia> @allocated RegularGrid(100000, 100000, 100000)
0

This is crucial in 3D interpolation for example where a simple 500x500x500 grid with x,y,z coordinates requires 1GB of memory. Because we are also modeling the subsurface with the package, you can always expect that our efforts will scale well.

Check the GeoStatsTutorials and let us know if a feature or method is missing.

7 Likes

As I mentioned in the first link you posted, ScatteredInterpolation.jl can work with a bunch of limitations. Comparing to what exists in Python and MATLAB, there is no simple linear or nearest-neighbor scattered interpolation, and because of the underlying function calls the Distance package, for large grid you will fall into out of memory error. It would be really nice to have such functionalities, and maybe me or someone interested should contribute to that package.

@henry2004y the out of memory error is not caused by the dependency on Distances.jl it is about how it is used internally in the package. We never construct full pairwise distance matrices in GeoStats.jl exactly because of that.

1 Like

@TheCedarPrince I highly recommend learning more about geostatistics if you are doing work related to geospatial data. All these interpolation, simulation, etc problems are addressed properly by the theory, and more contributors are always welcome.

5 Likes

Thanks for pointing out this package again! I think the reason I end up using the Python call original approach is simply that I was worried about the dependencies that come with the GeoStats.jl package. The name itself suggests that it may contain many things that is irrelevant to my purpose. It would be wonderful if these basic functionalities exist in a more lower level package like ScatteredInterpolation.jl.

@henry2004y you mean the submodules of the package like GeoEstimation.jl, KrigingEstimators.jl etc? The package is modular if you read the source code of GeoStats.jl you will see it just reloads submodules in other repos.

But I agree with you it is a big project, and that dependencies can be an issue when someone is only interested in simple interpolations without all the GIS complications.

The main advantage of sticking with GeoStats.jl is when you have actual maps with different coordinates, different spatial objects playing with each other, etc. And you want to replicate spatial processes with some given spatial structure.

Actually I just realize that it comes as submodules :sweat_smile: However, the package name itself is still not informative enough for me. As someone who is looking for a proper solution for plotting, it is really unintuitive to think at first that the desired function lies underneath the GeoStats.jl package.

I get it. And if the goal is just plotting, then as I said, GeoStats.jl too big of a hammer.

I attempted to get this example to work for my project (data in csv format and loaded into DataFrames). The working MWE from above example and sample plots are below.

#https://discourse.julialang.org/t/confusion-on-creating-gridded-interpolations/50205/15

using ScatteredInterpolation
using GLMakie

xs = [187.0, 313.0, 157.0, 343.0, 132.0, 368.0, 157.0, 343.0, 187.0, 313.0, 85.0, 415.0, 46.0, 454.0, 85.0, 415.0, 250.0, 250.0, 57.0, 443.0, 129.0, 371.0, 164.0, 336.0, 250.0, 250.0, 250.0, 180.0, 320.0, 118.0, 382.0]
ys = [57.0, 57.0, 134.0, 134.0, 250.0, 250.0, 366.0, 366.0, 443.0, 443.0, 129.0, 129.0, 250.0, 250.0, 371.0, 371.0, 132.0, 368.0, 313.0, 313.0, 415.0, 415.0, 403.0, 403.0, 417.0, 454.0, 475.0, 466.0, 466.0, 433.0, 433.0] 
samples = [-19.995152, -46.38638, -28.329832, 81.12224, -67.002556, 21.798035, 32.538525, 38.022156, -14.682513, -56.849937, -20.680933, 95.41259, 22.357513, 36.404854, -37.628937, -3.4352758, -25.271269, 36.05853, 26.666668, -29.088951, -7.4670925, -47.769047, -29.178823, -43.958958, 29.638771, -10.448129, 5.71033, -25.612432, 37.268364, -22.657557, -2.0058606]

points = hcat(xs, ys)'

len = 1001
x2 = range(0, 500, length=len)
y2 = range(0, 500, length=len)

points2 = hcat(x2, y2)'
#@show points2

x3 =[]
y3 = []
for x in x2
    for y in y2
        push!(x3, x)
        push!(y3, y)
    end
end
#x3 = vec(x3)
#y3 = vec(y3)
#@show x3
#@show y3
points2 = hcat(vec(x3), vec(y3))'
#@show points2

#z2 = [itp(x,y) for y in y2, x in x2]
itp = ScatteredInterpolation.interpolate(Multiquadratic(), points, samples)
#interpolated = ScatteredInterpolation.evaluate(itp, [186.0, 60.0]) # evaluate it for how ever many coordinates that you want
interpolated = ScatteredInterpolation.evaluate(itp, points2) # evaluate it for how ever many coordinates that you want

gridded = reshape(interpolated, len^2)

fig = Figure(resolution=(1500, 800), figure_padding = 100)
ax1 = Axis3(fig[1,1])
ax2 = Axis3(fig[1,2])
#@show interpolated
#@show gridded
x3 = collect(Iterators.flatten(x3))
y3 = collect(Iterators.flatten(y3))

meshscatter!(ax1, xs, ys, samples, markersize = 5)
surface!(ax2, x3, y3, interpolated, markersize=1)

fig

The MWE is kind of hacky at the moment - appreciate your feedback on code optimisation for the following sections or any codes that can be optimised.

len = 1001
x2 = range(0, 500, length=len)
y2 = range(0, 500, length=len)

points2 = hcat(x2, y2)'
#@show points2

x3 =[]
y3 = []
for x in x2
    for y in y2
        push!(x3, x)
        push!(y3, y)
    end
end
#x3 = vec(x3)
#y3 = vec(y3)
#@show x3
#@show y3
points2 = hcat(vec(x3), vec(y3))'

and

x3 = collect(Iterators.flatten(x3))
y3 = collect(Iterators.flatten(y3))

Thanks.

@Humphrey_Lee, you can try this:

using ScatteredInterpolation, GLMakie

# INPUT:
xs = [187., 313., 157., 343., 132., 368., 157., 343., 187., 313., 85., 415., 46., 454., 85., 415., 250., 250., 57., 443., 129., 371., 164., 336., 250., 250., 250., 180., 320., 118., 382.]
ys = [57., 57., 134., 134., 250., 250., 366., 366., 443., 443., 129., 129., 250., 250., 371., 371., 132., 368., 313., 313., 415., 415., 403., 403., 417., 454., 475., 466., 466., 433., 433.] 
samples = [-19.9, -46.3, -28.3, 81.1, -67.0, 21.7, 32.5, 38.2, -14.6, -56.8, -20.6, 95.4, 22.3, 36.4, -37.6, -3.4, -25.2, 36.5, 26.6, -29.8, -7.4, -47.7, -29.1, -43.9, 29.6, -10.4, 5.7, -25.6, 37.2, -22.6, -2.]

# DATA CONDITIONING:
points = hcat(xs, ys)'
nx = 1001
ny = 2000
x2 = LinRange(0, 500, nx)
y2 = LinRange(0, 1000, ny)
xy = Iterators.product(x2, y2)
xx = getindex.(xy,1)
yy = getindex.(xy,2)
points2 = hcat(vec(xx), vec(yy))'

# INTERPOLATION:
itp = ScatteredInterpolation.interpolate(Multiquadratic(), points, samples)
interpolated = ScatteredInterpolation.evaluate(itp, points2)

# PLOTTING:
fig = Figure(resolution=(1500, 800), figure_padding = 100)
ax1 = Axis3(fig[1,1])
ax2 = Axis3(fig[1,2])
meshscatter!(ax1, xs, ys, samples, markersize = 5)
surface!(ax2, x2, y2, reshape(interpolated,nx,ny))
3 Likes

Copying @joa-quim, because GMT has probably one of the most used algorithms in the world of geosciences for gridding scattered X, Y, Z(x,y) data points like in the above example.

The reference for the algorithm used by GMT’s surface function is:

Smith, W. H. F, and P. Wessel, 1990, Gridding with continuous curvature splines in tension, Geophysics , 55, 293-305.

1 Like

Awesome @rafael.guerra!. Like your codes. Just wondering if there are more elegant ways.
Noticed that you are a GP. I’m a PE (PP to be more specific).

How does GMT.jl compares against GeoStats.jl?

The two packages have different purposes but for this specific gridding problem, GMT surface solution (using splines under tension) would provide a deterministic solution while the GeoStats package would provide a probabilistic one (by exploiting any spatial correlations present by using krigging for instance).

See an interesting discussion on this topic here.

1 Like

Thanks for the summary.

I’ve tried the code as being described here to plot some Topographic EEG/MEG plot using Julia. The results are not quite satisfying, however.

The methods are described as follows in MatLab: griddata .

I’m wondering whether it results from the method difference to interpolate values?

@likanzhan, you may check other options in this other thread, in particular GMT.jl and DIVAnd.jl.

NB:
GMT’s splines in tension method also produces “surfaces that always pass through the data points

1 Like