Interpolating geodata on a sphere surface
I have a dataset with lat / lon coordinates and a corresponding 0/1 value for each geolocation (4 to 200+ data points). Now I want to interpolate voids and add colors to the surface of the globe based on the results of the interpolation. The main problem is to interpolate "around the world" because I am currently doing in a plane, which obviously does not work.
My details
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
Render 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)
2D in-plane interpolation
Currently, 2D spline interpolation is done directly in lat / lon coordinates using akima
papckage. This works, but does not take into account that the lat / lon coordinates lie on the 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 doesn't work for a sphere, since the left side doesn't match the right side. On a sphere, interpolation should be smooth.
What approaches can I use to interpolate on a sphere in R?
source to share
You can calculate the distances between points and the mesh yourself and then use your own interpolation. For example, below is the reciprocal distance interpolation in your sample data.
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 a raster for mesh 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
The function spDists
in the library sp
uses 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 interpolation using classic inverse distance interpolation with a power of 2 ( idp=2
). You can change your calculation if you want a different cardinality 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 write down 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)
source to share