The article treats the spherical interpolation and approximation, i.e., a way of measured data processing in case when the scalar data sampling is performed at nodes on the unit 2D sphere surface in the 3D Euclidean space (in general, on the (d−1)-dimensional sphere surface in the d-dimensional space). We use the spherical radial basis function interpolation with an inverse multiquadric basis function and a second degree polynomial in Cartesian coordinates as a trend.The formulae of this type can be useful in the interpretation of various physical measurements and have wide applications in geosciences, e.g., in the treatment of anisotropy of magnetic susceptibility data measured. We present the advantages of the formula chosen on numerical examples. However, in practical computation with high number of sampling nodes, the matrix of the system for determining interpolation coefficients may be ill-conditioned. Then the standard double precision LU factorization does not give a reliable solution due to the effect of round-off error accumulation and some more sophisticated means of solving the system are to be used.