Keywords. : Interpolation of scattered data, sphere, Shepard operator, least squares approximation, thin plate spline, inverse multiquadric, spiral points, error estimations.

MSC 2020 Subject Classification. : 41A05, 41A25, 41A80, 65D05, 65D15, 65D12, 33C55.
Abstract. We construct some smooth functions defined over a sphere that interpolate large sets of scattered data, using some modified Shepard methods, the least squares thin-plate spline and the inverse multiquadric functions. We illustrate the benefits of our methods in numerical examples for some test functions and some real data applications.

1 Preliminaries

We consider here the problem of interpolating large scattered data sets lying on a sphere by using a Shepard type method and some radial basis functions. This interpolation problem is important as it appears in solving some problems related to some physical phenomenons, such as temperature, rainfall, pressure, ozone or gravitational forces, measured at various points on the surface of the Earth; they could also be applied in modelling closed surfaces in CAGD, see, e.g., [ 7 ] . As it is mentioned in [ 9 ] and [ 25 ] , these type of data fitting problems, where the underlying domain is the sphere, arises in many areas, including, e.g., geophysics and meteorology, as, in general, the sphere is taken as a model of the Earth. Many authors have investigated the approximation of functions on the sphere by means of polynomials or radial basis functions (see, e.g, [ 1 ] , [ 2 ] , [ 7 ] - [ 12 ] , [ 17 ] , [ 19 ] , [ 25 ] ), the main motivation being the need to approximate geophysical quantities.

The Shepard method, introduced in [ 24 ] , is one of the best suited methods for approximating large sets of data. It has the advantages of a small storage requirement and an easy generalization to additional independent variables, but it suffers from no good reproduction quality, low accuracy and a high computational cost relative to some alternative methods [ 22 ] , these being the reasons for finding new methods that improve it (see, e.g., [ 3 ] - [ 6 ] , [ 13 ] , [ 20 ] , [ 27 ] ). The radial basis functions are suitable tools for scattered data, and for data that varies rapidly over short distances.

Our purpose in this paper is to introduce some combined Shepard methods using the least squares thin-plate spline and the inverse multiquadric functions, so our method involves also least squares approximation method, not just radial basis functions as, for example, in [ 9 ] , [ 12 ] and [ 25 ] .

Combined polynomial and radial basis function approximations have often been studied in the context of radial basis functions constructed from conditionally positive definite kernels, in which case a polynomial part is needed to make the theory work. Here, we restrict our attention to the case of (conditionally) strictly positive definite kernels, including a polynomial component, motivated by the fact that approximations of this kind offer real advantages.

The numerical tests on different types of data prove the efficiency of the method. Moreover, a physical phenomenon is investigated, namely temperature prediction on the Earth’s surface, and the results show that this method is a powerful instrument for solving various problems that model real life phenomena.

2 Combined spherical Shepard method

Following the idea and some notations from [ 7 ] , we introduce several modifications of the spherical Shepard method using spherical radial basis functions.

Let \(S\ \)be the unit sphere in \(\mathbb {R}^{3}\) centered at the origin. Let us consider the set of distinct nodes \(\mathbf{x}_{i} = (x_i, y_i, z_i),i=1,...,n\), lying on \(S\) and the corresponding function values \(f_{i},i=1,...,n\), with \(f: S \to \mathbb {R}\). For \(\mathbf{x}=(x,y,z) \in S\) the modified spherical Shepard operator is defined as (see, e.g., [ 7 ] )

with

and \(R_{j}^{w}\) is a radius of influence about the node \(j\) and the geodesic distance \(g\) between \(\mathbf{x}\) and \(\mathbf{x}_j\) defined as

where \(u_+ = \max \{ 0, u \} \) and \('' \cdot ''\) denotes the Euclidean inner product of the two vectors.

For a set of scattered data \(X=\{ \mathbf{x}_1,\mathbf{x}_2,...,\mathbf{x}_n\} \) the mesh norm \(h_{X}\) is defined by \(h_{X}=\sup \limits _{x\in S}g(x,X),\) where \(g\) denotes the geodesic distance.

If we denote by (see, e.g., [ 7 ] )

then we have the following properties of the weights:

1. the cardinality property: \(\overline{w}_j(\mathbf{x}_i) = \delta _{ij}\);

2. the partition of unity property: \(\sum \limits _{j=1}^{n}\overline{w}_j(\mathbf{x}) = 1\).

We recall some notions from [ 19 ] . A polynomial \(p:\mathbb {R}^3 \to \mathbb {R}\) of degree \(d,\; d \geq 0\) is homogeneous of degree \(d\) if \(p(t \mathbf{x}) = t^d p(\mathbf{x})\), for any \( \mathbf{x} \in \mathbb {R}^3 \) and \(t{\gt}0\). The polynomial \(p\) is harmonic if \(\Delta p(\mathbf{x}) = 0\), for any \(\mathbf{x} \in \mathbb {R}^3\), where \(\Delta \) is the Laplace operator, i.e.,

So, if \(\mathcal{P}_d\) is the set of polynomials of degree \(d\) in \(\mathbb {R}^3\) that are harmonic and homogeneous of order \(d\), we can define the set

of spherical harmonics of exact order \(d\). The dimension of \(\mathcal{H}_{d}^{*}(S)\) is

Denoting the space of spherical harmonics of maximum order \(d\) by \( \mathcal{H}_{d}(S) \), we have

It follows that the dimension of \( \mathcal{H}_{d}(S) \) is

Now we introduce the new local Shepard interpolants,

with the local zonal basis function given by

and

with

considering the thin-plate spline spherical radial basis function (see, e.g., [ 2 ] )

with the local zonal basis function given by

and

with

considering the inverse multiquadric spherical radial basis function (see, e.g., [ 2 ] )

Since,

we get

3 Test results

• Numerical results on test functions

We use the following test functions (see, e.g., [ 15 ] , [ 16 ] , [ 21 ] , [ 22 ] ) to illustrate the theoretical results:

We consider two cases for the choice of nodes and points on which the operators are constructed.

For the first case, we evaluate the operators on 1000 random points on the unit sphere \(S\), using sets of 500, 1000 and 2000 spiral nodes, respectively. These nodes are constructed based on a method proposed in [ 23 ] . To describe them, we use their spherical coordinates

with

In Tables 1 and 2 we list the maximum absolute errors (MAEs) and the root mean square errors (RMSEs) for the thin-plate spline case, given in (2.9) and for the inverse multiquadric case, given in (2.11).

The shape parameter \(c\) in the inverse multiquadric case, we have computed as \(1/(0.815d),\) where \(d = \frac{1}{n} \sum \limits _{i=1}^n d_i\), with \(d_i\) being the distance from the node \(i\) to its closest neighbour, following the idea proposed by Hardy in [ 18 ] . We also used the method proposed in [ 14 , Program 2 ] , but the results are comparable.

We consider also the case of zonal basis functions combined with spherical harmonics. For \(S_2\) and for \(S_4\), given in (2.10) and (2.12), we take the particular cases

Tables 3 and 4 contain the approximation errors for \(S_2\) and \(S_4\).

Numerical experiments have shown that in our cases, a good value for the number of closest neighbours of a node is \(n_Z=17\), but its choice is not an easy task, since there are many factors that could influence the numerical results, as it is also emphasized in [ 10 ] .

We also show the computational times (CPU) of the methods considered above. Comparing the approximation results, one can see that in most of the cases, the thin-plate spline function is slightly more accurate, although good results were obtained in the case of inverse multiquadric function as well. For both zonal basis functions, the addition of a spherical harmonic produces an improved accuracy in the majority of cases, the computational times being quite equivalent.

In the second case, for a fair distribution of points on the sphere, as proposed in [ 7 ] , we evaluate the operator on a set of 600 spiral points, using sets of 1000, 2000, 5000 Halton nodes, respectively. The Halton nodes are constructed following the method proposed in [ 26 ] . The approximation errors for the operators \(S_1\), \(S_3\) are displayed in Tables 5 and 6 and the results for the particular case of operators combined with a spherical harmonic \(S_2\) and \(S_4\) are given in Tables 7 and 8. The same advantages as in the first case can be observed.

To obtain a real visualisation effect of the numerical results, in Figures 4, 8 we plot the functions values \(f_i\), \(i=1,2\), and the values of the corresponding Shepard interpolants \(S_1f_i\) and \(S_3f_i\), \(i=1,2\), considering 1000 Halton nodes and 20000 spiral points on the unit sphere. In Table 9 we list the MAEs and the RMSEs in this case, for all the functions \(f_i,\; i=1,...,4\).

Finally, we compared our results with the ones obtained in the case of the modified spherical Shepard operator \(S\), defined in (2.1). In Table 10 we have the MAEs and RMSEs for the case of spiral nodes and 1000 random points and in Table 11 we have the errors for the case of Halton nodes and 600 spiral points. Comparing the results, one can observe the benefits of choosing the combined operators with the zonal basis functions, since in the majority of cases, the interpolation results are more accurate.

The experiments were performed in Matlab.

• Numerical results in a real data application

We consider an example of approximation for some real data, namely the monthly-mean temperatures on the Globe in January 2010 and June 2010. The set of data was selected from [ 28 ] . For our numerical tests we considered \(1073\) nodes and we reconstructed the temperature values for \(21449\) points on the Globe, after we projected the latitude and longitude of each point onto the unit sphere \(S\). A representation for the sets of nodes and points is given in Figure 11.

To handle the global part, we considered the case of zonal basis functions combined with a spherical harmonic, \(S_2\) and \(S_4\). For the inverse multiquadric case, we computed the parameter \(c\) as in the previous numerical tests (see, e.g., [ 18 ] ). The value obtained was \(c=35.26\). In Table 12 there are given the mean absolute errors (MAEs) and the root mean square errors (RMSEs). Figures 15, 19 present the graphical results for the temperatures in January 2010 and June 2010, respectively. Studying the results obtained for this example, there are observed and confirmed the advantages of using the previously introduced combined Shepard operators to real data applications.

1

Allasia, G., Cavoretto, R., De Rossi, A.: Hermite-Birkhoff interpolation on scattered data on the sphere and other manifolds. Appl. Math. Comput. 318, 35–50 (2018)

2

Baxter, B.J.C., Hubbert, S.: Radial basis functions for the sphere. In: Recent Progress in Multivariate Approximation, Witten-Bommerholz, 2000. Internat. Ser. Numer. Math., vol. 137, pp. 33–47. Birkhäuser, Basel (2001)

3

Cătinaş, T.: The combined Shepard-Abel-Goncharov univariate operator. Rev. Anal. Numér. Théor. Approx. 32, 11–20 (2003)

4

Cătinaş, T.: The combined Shepard-Lidstone bivariate operator. In: Trends and Applications in Constructive Approximation. International Series of Numerical Mathematics (de Bruin, M.G. et al. (eds.)), 151, pp. 77–89, Springer Group-Birkhäuser Verlag (2005)

5

Cătinaş, T.: Bivariate interpolation by combined Shepard operators. In: Proceedings of \(17^{th}\) IMACS World Congress, Scientific Computation, Applied Mathematics and Simulation (Borne, P., Benrejeb, M., Dangoumau, N., Lorimier, L. (eds.)), 7 pp., Paris, July 11–15 (2005)

6

Cătinaş, T.: The bivariate Shepard operator of Bernoulli type. Calcolo 44(4), 189–202 (2007)

7

De Rossi, A.: Spherical interpolation of large scattered data sets using zonal basis functions. In: Mathematical Methods for Curves and Surfaces (Daehlen, M., Morken, K., Schumaker, L. (eds.)), pp. 125–134, Trømso 2004. Nashboro Press (2005)

8

De Rossi, A.: Hybrid spherical approximation. arXiv:1404.1475 (2014)

9

Cavoretto, R., De Rossi, A.: Fast and accurate interpolation of large scattered data sets on the sphere. J. Comput. Appl. Math. 234, 1505–1521 (2010)

10

Cavoretto, R., De Rossi, A.: Numerical comparison of different weights in Shepard’s interpolants on the sphere. Appl. Math. Sci. 4, 3425–3435 (2010)

11

Cavoretto, R., De Rossi, A.: Spherical interpolation using the partition of unity method: An efficient and flexible algorithm. Appl. Math. Lett. 25, 1251–1256 (2012)

12

Cavoretto, R., De Rossi, A.: Achieving accuracy and efficiency in spherical modelling of real data. Math. Methods Appl. Sci. 37, 1449–1459 (2014)

13

Farwig, R.: Rate of convergence of Shepard’s global interpolation formula. Math. Comp. 46, 577–590 (1986)

14

Fasshauer, G.E., Zhang J.G.: On choosing “optimal” shape parameters for RBF approximation. Numer. Algor. 45, 345–368 (2007)

15

Franke, R.: Scattered data interpolation: tests of some methods. Math. Comp. 38, 181–200 (1982)

16

Franke, R., Nielson, G.: Smooth interpolation of large sets of scattered data. Int. J. Numer. Meths. Engrg. 15, 1691–1704 (1980)

17

Le Gia, Q.T., Sloan, I. H., Wendland, H.: Multiscale analysis in Sobolev spaces on the sphere. SIAM J. Numer. Anal. 48(6), 2065–2090 (2010)

18

Hardy, R.L.: Multiquadric equations of topography and other irregular surfaces. J. Geophys. Res. 76, 1905–1915 (1971)

19

Hubbert, S., Le Gia, Q.T., Morton, T.: Spherical Radial Basis Functions, Theory and Applications. Springer (2015)

20

Lazzaro, D., Montefusco, L.B.: Radial basis functions for multivariate interpolation of large scattered data sets. J. Comput. Appl. Math. 140, 521–536 (2002)

21

Renka, R.J., Cline, A.K.: A triangle-based \(C^{1}\) interpolation method. Rocky Mountain J. Math. 14, 223–237 (1984)

22

Renka, R.J.: Multivariate interpolation of large sets of scattered data. ACM Trans. Math. Software 14, 139–148 (1988)

23

Saff, E.B., Kuijlaars, A.B.J.: Distributing many points on a sphere. Math. Intelligencer. 19, 5–11 (1997)

24

Shepard, D.: A two dimensional interpolation function for irregularly spaced data. In: Proceedings of the 1968 23rd ACM National Conference, pp. 517–523. ACM (1968)

25

Sloan, I.H., Sommariva, A.: Approximation on the sphere using radial basis functions plus polynomials. Adv. Comput. Math., 29, 147–177 (2008)

26

Wong, T.T., Luk, W.S., Heng, P.A.: Sampling with Hammersley and Halton Points. J. Graph. Tools, 2(2), 9-24 (1997)

27

Zuppa, C.: Error estimates for moving least square approximations. Bull. Braz. Math. Soc. (N.S), 34(2), 231–249 (2003)

28

https://www.kaggle.com/datasets/shishu1421/global-temperature?select=air_temp.2010https://www.kaggle.com/datasets/shishu1421/global-temperature?select=air_temp.2010