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}, StaticArraysCore.SMatrix{4, 4, Int64, 16}}([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.0and 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.235803568569594
47.179415507626665
47.0932077352213
46.98908834882418
46.87896544589783
46.774747123914494
46.68834148033668
46.63165661263463
46.61391092263422
⋮
39.09483624644603
39.09666708631454
39.098177071127225
39.09891621696852
39.09843453992283
39.09628205607458
39.09200878150819
39.08516473230809
39.07529992455869Plotting 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}, StaticArraysCore.SMatrix{4, 4, Int64, 16}}([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.646999e6To 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.1086Plotting 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)