I have a dataset with lat/lon coordinates and a corresponding 0/1 value for each geolocation (4 to 200+ datapoints). Now, I want to interpolate the voids and add colors to the surface of the globe based on the interpolation results. The main problem I have is to interpolate "around the globe", because currently I do in a plane, which obviously does not work.

My data

set.seed(41)
n <- 5
s <- rbind(data.frame(lon = rnorm(n, 0, 180),
                      lat = rnorm(n, 90, 180),
                      value = 0),
           data.frame(lon = rnorm(n, 180, 180),
                      lat = rnorm(n, 90, 180),
                      value = 1))
s$lon <- s$lon %% 360 -180
s$lat <- s$lat %% 180 -90
s_old <- s

Visualize datapoints

library(sp)  
library(rgdal)
library(scales)
library(raster)
library(dplyr)

par(mfrow=c(2,1), mar=c(0,0,0,0))
grd <- expand.grid(lon = seq(-180,180, by = 20), 
                   lat = seq(-90, 90, by=10))
coordinates(grd) <- ~lon + lat
gridded(grd) <- TRUE
plot(grd, add=F, col=grey(.8))

coordinates(s) = ~lon + lat
points(s, col=s$value + 2, pch=16, cex=.6)

Bivariate interpolate in plane

Currently, the bivariate spline interpolation is done directly on the lat/lon coordinates using the akima papckage. This works, but does not take into account that the lat/lon coordinates lie on a sphere.

nx <- 361
ny <- 181
xo <- seq(-180, 179, len=nx)
yo <- seq(-90, 89, len=ny)
xy <- as.data.frame(coordinates(s))
int <- akima:::interp(x = xy$lon, y = xy$lat, z = s$value, 
                      extrap = T, 
                      xo = xo, yo = yo, 
                      nx = nx, ny=100, 
                      linear = F)
z <- int$z
# correct for out of range interpolations values
z[z < 0] <- 0
z[z > 1] <- 1

grd <- expand.grid(lon = seq(-180,180, by = 20), 
                   lat = seq(-90, 90, by=10))
coordinates(grd) <- ~lon + lat
gridded(grd) <- TRUE
plot(grd, add=F, col=grey(.8))

## create raster image
r <- raster(nrows=ny, ncols=nx, crs='+proj=longlat',
            xmn=-180, xmx=180, ymn=-90, ymx=90)
values(r) <- as.vector(z)  

# tweaking of color breaks
colors <- alpha(colorRampPalette(c("red", "yellow", "green"))(21), .4)
br <- seq(0.3, 0.7, len=20)
image(xo, yo, z, add = T, col = colors, breaks=c(-.1, br, 1.1))
points(s, col=s$value + 2, pch=16, cex=.6)

Obviously this does not work for a sphere, as the left side does not match the right side. On a sphere,the interpolation should be seamless.

What approaches can I use do interpolate on a sphere in R?

You can calculate distances between points and grid yourself and then use your own interpolation. For instance, below is a inverse distance interpolation on your data example.

Generate Data

library(sp)
library(rgdal)

# Data
set.seed(41)
n <- 5
s <- rbind(data.frame(lon = rnorm(n, 0, 180),
                      lat = rnorm(n, 90, 180),
                      value = 0),
           data.frame(lon = rnorm(n, 180, 180),
                      lat = rnorm(n, 90, 180),
                      value = 1))
s$lon <- s$lon %% 360 -180
s$lat <- s$lat %% 180 -90
s_old <- s

Create raster for grid interpolation

## create raster image
r <- raster(nrows=ny, ncols=nx, crs='+proj=longlat',
            xmn=-180, xmx=180, ymn=-90, ymx=90)

Calculate distances between points and raster

Function spDists in library sp use the Great Circle Distance when coordinates are not projected. This means that the distance calculated between two points is the shortest.

# Distance between points and raster
s.r.dists <- spDists(x = coordinates(s), y = coordinates(r), longlat = TRUE)

Interpolate on a sphere using inverse distance interpolation

Here I propose to interpolate using the classical inverse distance interpolation with power 2 (idp=2). You can modify yourself the calculation if you want other power or linear interpolation, or if you want to interpolate with a limited number of neighbors.

# Inverse distance interpolation using distances
# pred = 1/dist^idp
idp <- 2
inv.w <- (1/(s.r.dists^idp))
z <- (t(inv.w) %*% matrix(s$value)) / apply(inv.w, 2, sum)

r.pred <- r
values(r.pred) <- z

Then plot the results

# tweaking of color breaks
colors <- alpha(colorRampPalette(c("red", "yellow", "green"))(21), .4)
br <- seq(0.3, 0.7, len=20)
plot(r.pred, col = colors, breaks=c(-.1, br, 1.1), legend=F)
points(s, col=s$value + 2, pch=16, cex=.6)