In this paper a local approximation method on the sphere is presented. As
interpolation scheme we consider a partition of unity method, such as the
modified spherical Shepard's method, which uses zonal basis functions (ZBFs)
plus spherical harmonics as local approximants. Moreover, a spherical zone
algorithm is efficiently implemented, which works well also when the amount of
data is very large, since it is based on an optimized searching procedure.
Numerical results show good accuracy of the method, also on real geomagnetic
data