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.

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.

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.

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

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?

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.)