B-spline interpolation function

If the B-spline structure is defined, it can be used to define the B-spline interpolation function. In this chapter, we will continue with the examples from the B-spline structure section for the one and two dimensional case.

One dimensional case

In the B-spline structure section, we began with the example file rhine_data_cubic-nak.jl from the examples folder where we already defined the B-spline structure for the cubic B-spline interpolation with not-a-knot end condition and smoothing.

# Include packages
using TrixiBottomTopography
using CairoMakie

# Define data path
root_dir = pkgdir(TrixiBottomTopography)
data = joinpath(root_dir, "examples", "data", "rhine_data_1d_20_x.txt")

# Define B-spline structure
spline_struct = CubicBSpline(data; end_condition = "not-a-knot", smoothing_factor = 999)

CubicBSpline{Vector{Float64}, Float64, Vector{Float64}, Matrix{Int64}}([357000.0, 357020.0, 357040.0, 357060.0, 357080.0, 357100.0, 357120.0, 357140.0, 357160.0, 357180.0  …  357800.0, 357820.0, 357840.0, 357860.0, 357880.0, 357900.0, 357920.0, 357940.0, 357960.0, 357980.0], 20.0, [46.91426684256018, 47.397904854124576, 46.99689666444544, 46.50889411070535, 46.73154903008692, 46.862924837363984, 46.74903591362994, 46.40702226591983, 46.31904390610859, 46.654018220049146  …  39.07284424693978, 39.1231081223708, 39.13738000648763, 39.136760208084766, 39.10325245367299, 39.08822484618267, 39.09410545889054, 39.1030246596446, 39.08484073202447, 39.009411959609665], [-1 3 -3 1; 3 -6 3 0; -3 0 3 0; 1 4 1 0])

To define the B-spline interpolation function for a variable x, we use the spline_interpolation function. For the one dimensional case, this function is implemented as spline_interpolation(linear_struct, x) and spline_interpolation(cubic_struct, x).

# Define B-spline interpolation function
spline_func(x) = spline_interpolation(spline_struct, x)

This defines the cubic B-spline interpolation function with not-a-knot end condition and smoothing with respect to variable x from the previously created spline_struct. If we want to visualize the interpolation function with 100 interpolation points, we define the following:

# Define interpolation points
n = 200
x_int_pts = Vector(LinRange(spline_struct.x[1], spline_struct.x[end], n))

200-element Vector{Float64}:
 357000.0
 357004.92462311557
 357009.84924623114
 357014.77386934677
 357019.6984924623
 357024.6231155779
 357029.5477386934
 357034.47236180905
 357039.3969849247
 357044.3216080402
      ⋮
 357940.6030150754
 357945.52763819095
 357950.4522613065
 357955.3768844221
 357960.3015075377
 357965.2261306533
 357970.15075376886
 357975.07537688443
 357980.0

and evaluate to obtain the corresponding y values by:

# Get interpolated values
y_int_pts = spline_func.(x_int_pts)

200-element Vector{Float64}:
 47.25046382058398
 47.2358035685696
 47.179415507626665
 47.0932077352213
 46.98908834882418
 46.87896544589783
 46.774747123914494
 46.68834148033668
 46.63165661263463
 46.61391092263422
  ⋮
 39.09483624644603
 39.09666708631454
 39.098177071127225
 39.098916216968526
 39.09843453992283
 39.09628205607458
 39.09200878150819
 39.08516473230809
 39.07529992455869

Plotting the interpolated points can be done via the command below and produces a one-dimensional plot

plot_topography(x_int_pts, y_int_pts; xlabel = "ETRS89 East", ylabel = "DHHN2016 Height")

Alternatively, one can plot the interpolated bottom topography together with the interpolation knots as follows

# Get the original interpolation knots
x_knots = spline_struct.x
y_knots = spline_func.(x_knots)

plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, x_knots, y_knots;
                                         xlabel = "ETRS89 East", ylabel = "DHHN2016 Height")

Two dimensional case

In the B-spline structure section, we examined the example file rhine_data_bicubic-nak.jl from the examples folder where we already created the B-spline structure for the bicubic B-spline interpolation with not-a-knot end condition and smoothing.

# Include packages
using TrixiBottomTopography
using CairoMakie

# Define data path
root_dir = pkgdir(TrixiBottomTopography)
data = joinpath(root_dir, "examples", "data", "rhine_data_2d_20.txt")

# Define B-spline structure
spline_struct = BicubicBSpline(data; end_condition = "not-a-knot", smoothing_factor = 9999)

BicubicBSpline{Vector{Float64}, Vector{Float64}, Float64, Matrix{Float64}, Matrix{Int64}}([357000.0, 357020.0, 357040.0, 357060.0, 357080.0, 357100.0, 357120.0, 357140.0, 357160.0, 357180.0  …  357800.0, 357820.0, 357840.0, 357860.0, 357880.0, 357900.0, 357920.0, 357940.0, 357960.0, 357980.0], [5.646019e6, 5.646039e6, 5.646059e6, 5.646079e6, 5.646099e6, 5.646119e6, 5.646139e6, 5.646159e6, 5.646179e6, 5.646199e6  …  5.646819e6, 5.646839e6, 5.646859e6, 5.646879e6, 5.646899e6, 5.646919e6, 5.646939e6, 5.646959e6, 5.646979e6, 5.646999e6], 20.0, [46.600451320002854 47.282399695552705 … 47.09971997322787 47.10735082392629; 45.82016087061047 46.502109246160636 … 46.963741530201055 46.97137238089921; … ; 45.99186284178017 45.686484464847375 … 39.098457400873656 39.23988800954714; 46.243985904572256 45.938607527639306 … 39.06537978820518 39.206810396878865], [-1 3 -3 1; 3 -6 3 0; -3 0 3 0; 1 4 1 0])

To define the B-spline interpolation function for variables x and y, we use the spline_interpolation function for the two dimensional case. This functionality is implemented as spline_interpolation(bilinear_struct, x, y) and spline_interpolation(bicubic_struct, x, y).

# Define B-spline interpolation function
spline_func(x,y) = spline_interpolation(spline_struct, x, y)

This defines the bicubic B-spline interpolation function with not-a-knot end condition and smoothing for variables x and y due to the previously constructed spline_struct. If we want to visualize the bicubic interpolation function with 100 interpolation points in each spatial direction, we define:

# Define interpolation points
n = 100
x_int_pts = Vector(LinRange(spline_struct.x[1], spline_struct.x[end], n))
y_int_pts = Vector(LinRange(spline_struct.y[1], spline_struct.y[end], n))

100-element Vector{Float64}:
 5.646019e6
 5.646028898989899e6
 5.646038797979798e6
 5.646048696969697e6
 5.646058595959596e6
 5.646068494949495e6
 5.646078393939394e6
 5.646088292929293e6
 5.646098191919192e6
 5.646108090909091e6
 ⋮
 5.646919808080808e6
 5.646929707070707e6
 5.646939606060606e6
 5.646949505050505e6
 5.646959404040404e6
 5.646969303030303e6
 5.646979202020202e6
 5.646989101010101e6
 5.646999e6

To fill a matrix z_int_pts, which contains the corresponding z values for x_int_pts and y_int_pts, we use the helper function evaluate_bicubicspline_interpolant.

# Get interpolated matrix
z_int_pts = evaluate_bicubicspline_interpolant(spline_func, x_int_pts, y_int_pts)

100×100 Matrix{Float64}:
 46.4652  46.0495  45.59    45.0608  …  45.5853  45.6011  45.6461  45.719
 46.7486  46.3415  45.8956  45.3881     45.556   45.5411  45.557   45.6029
 46.9957  46.6121  46.1894  45.7059     45.5651  45.5146  45.4976  45.5147
 47.1991  46.8501  46.4565  46.0011     45.6092  45.5245  45.4746  45.4658
 47.3511  47.0408  46.6819  46.2619     45.6843  45.5739  45.4969  45.468
 47.4475  47.1729  46.8535  46.4774  …  45.79    45.6678  45.5745  45.5312
 47.4999  47.2541  46.9726  46.642      45.9414  45.8197  45.717   45.6593
 47.5231  47.2966  47.0444  46.7527     46.1512  46.0398  45.9296  45.8516
 47.5254  47.3085  47.0747  46.8123     46.3974  46.3073  46.1955  46.0947
 47.5134  47.2972  47.0707  46.8264     46.6459  46.5899  46.4879  46.369
  ⋮                                  ⋱                             
 46.9407  46.8666  46.7907  46.7222     39.06    39.0671  39.071   39.072
 46.8836  46.8109  46.7421  46.682      39.0394  39.049   39.0552  39.0585
 46.848   46.7785  46.7166  46.6635     39.019   39.031   39.0395  39.0445
 46.8304  46.7651  46.7088  46.6602     39.0024  39.0165  39.027   39.0328
 46.8323  46.7707  46.7164  46.6678  …  38.9936  39.0082  39.0196  39.0259
 46.854   46.7937  46.7365  46.6822     38.9969  39.0093  39.02    39.0263
 46.8873  46.8252  46.7617  46.6978     39.0192  39.0257  39.033   39.0377
 46.922   46.8547  46.7832  46.7084     39.0675  39.0638  39.0642  39.0638
 46.9479  46.8731  46.7926  46.7075     39.1495  39.1302  39.1184  39.1086

Plotting the interpolated values gives the representation directly below the plot_topography command

plot_topography(x_int_pts, y_int_pts, z_int_pts;
                xlabel="ETRS89\n East",
                ylabel="ETRS89\n North",
                zlabel="DHHN2016\n Height",
                azimuth_angle = 54 * pi / 180,
                elevation_angle = 27 * pi / 180)

Alternatively, one can plot the interpolated two-dimensional bottom topography together with the interpolation knots as follows

# Get the original interpolation knots
x_knots = spline_struct.x
y_knots = spline_struct.y
z_knots = evaluate_bicubicspline_interpolant(spline_func, x_knots, y_knots)

plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, z_int_pts,
                                         x_knots, y_knots, z_knots;
                                         xlabel="ETRS89\n East",
                                         ylabel="ETRS89\n North",
                                         zlabel="DHHN2016\n Height",
                                         azimuth_angle = 54 * pi / 180,
                                         elevation_angle = 27 * pi / 180)