Interpolate 2d function

Specific Domains Numerics
, ,
1

I want to interpolate 2d data. I have a function f(x, y) that is sampled as follows. For any y value in some array y_arr, I have the a set of points x_arr. Note that x_arr is different for any y, and even the boundaries of the x values are not the same.

I would like to perform interpolation of this 2d function.

For 1d interpolations, I am using BSplineKit, but from its documentation, I could not understand if it is possible to use it for 2d.

Would love to get any suggestions!

2

The README in Dierckx package seems to promise 2D interpolation on an irregular grid:
https://github.com/kbarbary/Dierckx.jl

3

True. Though I am having problems with this package.

Also, recently I saw that it is not actively maintained, so, therefore, I am looking for another solution.

4

You might try

in particular if you know beforehand the shape of the domain where you want to interpolate and which does not need to be a rectangle.

2 Likes
5

While Dierckx isn’t actively maintained, the existing features should work correctly, including 2D interpolation. New features may be slow to get added, that’s what the maintenance comment might have been about.

2 Likes
6

Thanks.

I tried to use Dierckx and I get this error:

No more knots can be added because the additional knot would
(quasi) coincide with an old one: s too small or too large a weight to
an inaccurate data point.  The weighted least-squares spline
corresponds to the current set of knots.

I tried to play with s, and it did not help. Not sure what the problem is here.
I am fine also with the simplest linear interpolation between points.

7

Multidimensional interpolations, using tensor product B-splines, have been in the plans for BSplineKit.jl for a while. There’s even a slighly outdated open PR with a possible implementation, but I’m not very happy with the proposed interface and I haven’t had the time to do something better.

If it helps, here it’s a bit of code for manually doing 2D interpolations with BSplineKit.jl:

using BSplineKit
using LinearAlgebra

# Target function
f(x, y) = 2.0 + sinpi(x) * cospi(y)

xs = 0:0.1:1
ys = 0:0.1:2
fdata = f.(xs, ys')

ord = BSplineOrder(4)  # cubic splines

# Create B-spline knots based on interpolation points (uses an internal function)
ts_x = SplineInterpolations.make_knots(xs, ord, nothing)
ts_y = SplineInterpolations.make_knots(ys, ord, nothing)

# Create B-spline bases
Bx = BSplineBasis(ord, ts_x; augment = Val(false))
By = BSplineBasis(ord, ts_y; augment = Val(false))

# Create and factorise interpolation matrices
Cx = lu!(collocation_matrix(Bx, xs))
Cy = lu!(collocation_matrix(By, ys))

# 2D B-spline coefficients (output)
coefs = similar(fdata)

# Solve linear systems
for j ∈ eachindex(ys)
    @views ldiv!(coefs[:, j], Cx, fdata[:, j])
end
for i ∈ eachindex(xs)
    @views ldiv!(Cy, coefs[i, :])
end

# Evaluate 2D (tensor product) spline at point (x, y).
function eval_spline2D(coefs::AbstractMatrix, (Bx, By), (x, y))
    i, bx = Bx(x)  # evaluate B-spline basis at point x
    j, by = By(y)
    kx = order(Bx)
    ky = order(By)
    val = zero(eltype(coefs))  # spline evaluated at (x, y)
    for δj ∈ 1:ky, δi ∈ 1:kx
        coef = coefs[i - δi + 1, j - δj + 1]
        bi = bx[δi]
        bj = by[δj]
        val += coef * bi * bj
    end
    val
end

# Verification: evaluate 2D spline at data points
for m ∈ eachindex(ys), n ∈ eachindex(xs)
    x = xs[n]
    y = ys[m]
    val = eval_spline2D(coefs, (Bx, By), (x, y))
    fval = fdata[n, m]
    @assert val ≈ fval  # check that the value is correct
end
2 Likes
8

Thanks @jipolanco !

This code seems like it is designed for an even grid, isn’t it?
I have a large set of x and y, and the values z(x,y). So this is an irregular grid. Is there a way to perform the 2d interpolation with the code you attached?

9

The example uses a regular grid, but it should work without a problem on irregular ones. (Just define xs and ys as vectors with sorted irregular values.)

1 Like
10

But how can it be sorted? I essentially have a set (x_i, y_i) and the values z_i = z(x_i, y_i).

11

Maybe I misunderstood your initial post. If you have a curvilinear grid (your x_i coordinates vary with y), then I’m not sure I can help…

1 Like
12

Linking also this other thread, which contains examples of GMT.jl and DIVAnd.jl solutions.

1 Like
13

Thanks! I will try those!

14

Also, you should check out ScatteredInterpolation.jl which seems to be exactly what you’re looking for.

1 Like