Computationally efficient statistical models are needed to understand behaviour of spatial phenomena observed at huge number of locations all over the globe. Spatial modelling for worldwide geostatistical data gives rise to issues on how to build valid covariance functions accounting for geodesic distances between locations, and can be computationally demanding.  In this work we undertake an efficient non parametric approach to model and estimate the spatial field in these situations. We extend the Bayesian P-spline approach for smoothing data on a sphere, taking the sphere a surrogate for the globe. The key idea in traditional P-splines is to model the spatial field as the sum of scaled B-spline basis functions, each basis being defined on a regular grid of knots.  We extend this idea for data on a sphere by adopting a geodesic grid system, giving a set of quasi-equally spaced knots on the sphere. The example motivating our research is taken from a real climate study, where the goal is to identify the Intertropical Convergence Zone using high-resolution remote sensing data.