Smooth approximation of point data


I’m writing CFD code to model the tides in estuaries. Let’s say I have some measurement data of the bathymetry on an irregular grid.

What is the best way in NGSolve to create a smooth (C1 or higher) CoefficientFunction/GridFunction that interpolates this data?

I have seen several other posts:
It is suggested to use a linear approximation of the data, by using a H1 space of order 1. However, I need to compute the Hessian of my solution and I have found that it is very sensitive to the slope changes of the triangles.

An other suggestion is to use the VoxelCoefficient functions. It seems that VoxelCoefficients are piecewise constant or interpolate the data linearly only.

Lastly, BSpline CoefficientFunctions are suggested. However, can BSpline CFs be used to represent a surface?

Could you provide a MWE?