# Spatial interpolation on Earth as a 3D sphere

## Spatial interpolation at Earth scale

Earth map is usually presented flat. Geographical coordinates usually go from -180° to +180°. Spatial interpolations using coordinates may be tricky as -180° is equal to +180°. Here I propose a way to realize spatial interpolations on Earth as a sphere and then map the outputs in 3D using `rgl`. The complete R-script extracted from this post is here on my github.

## Random dataset

Let’s create a random dataset with observations of presences and absences over the globe (This is simulated with a little spatial auto-correlation North / South). We do not account for possible difference between terrestrial and oceans areas.
I use a ‘modulo’ trick to restrict coordinates between -180° and +180°. This is just for me to keep a trace of it (otherwise I could have used `runif` to simulate data).

``````# library(sp)
library(rgdal)
library(sf)
library(raster)
library(ggplot2)

# Simulate dataset with a little spatial auto-correlation
set.seed(42)
n <- 30
obs <- rbind(data.frame(lon = rnorm(n, 0, 180),
lat = rnorm(n, 50, 35),
value = 0),
data.frame(lon = rnorm(n, 180, 180),
lat = rnorm(n, -50, 35),
value = 1))
# Modulo trick
obs\$lon <- obs\$lon %% 360 -180
obs\$lat <- obs\$lat %% 180 -90

# Transform points as spatial points with sf ----
obs_sf <- st_as_sf(obs, coords = c("lon", "lat"),
crs = 4326)

# Plot on a worldmap
worldmap <- borders("world",
colour = "#fefefe",
fill = "#808080"
)

# Plot over a worldmap
ggplot() + worldmap +
geom_point(data = obs,
aes(x = lon, y = lat,
colour = factor(value))) +
coord_quickmap()``````

## Inverse distance interpolation

A simple interpolation is the inverse distance interpolation. We can calculate distances between observations and world map using great circle.

• Create an empty `raster` that will be used as a target for interpolations.
• Transform points as spatial points with library `sf`
• Calculate (great circle) distances between points and raster cells
• Interpolate using inverse distance interpolation
``````# Create an empty world raster ----
ny <- 41
nx <- 80
r <- raster(
nrows = ny, ncols = nx,
crs = '+proj=longlat',
xmn = -180, xmx = 180,
ymn = -90, ymx = 90
)

# Transform raster as spatial points with sf
r_sf <- st_as_sf(as.data.frame(coordinates(r)),
coords = c("x", "y"),
crs = 4326)

# Distance between points and raster ----
obs.r.dists <- st_distance(obs_sf, r_sf)
obs.r.dists <- unclass(obs.r.dists)

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

# Fill in raster for predictions
r.pred <- r
values(r.pred) <- z

# Plot prediction raster
worldmap_predict <- borders("world",
colour = "#05050541",
fill = NA,
size = 0.5
)

rasterVis::gplot(r.pred) +
geom_tile(aes(fill = value)) +
scale_fill_gradient2(low = 'red', mid = "yellow",
high = 'green', midpoint = 0.5) +
geom_point(data = obs,
aes(x = lon, y = lat, colour = factor(value))) +
worldmap_predict +
guides(colour = FALSE) +
coord_quickmap()``````

## Kriging with user defined distances

To interpolate the dataset, you may want to use kriging. This interpolation relies on a variogram allowing to make use of the spatial structure of the data (contrary to idw which is fixed by the user).
For interpolations with custom distances and obstacles, I created the (totally in-development) `GeoDist` R-library. Here we do not have obstacles, but we have a user-defined distance matrix based on great circle distances that we want to use to fit a variogram and interpolate on the map. In GeoDist, I modified some functions of `geoR` so that it can fit a variogram and interpolate using a distance matrix defined by the user.

We calculate the variogram according to our own great circle distance matrix with `GeoDist::variog.dist`.

``````library(GeoDist)
# Transform data as geodata for geoR library
datageo <- as.geodata(obs, data.col = "value")
# Calculate distances between observations
obs.obs.dists <- st_distance(obs_sf)
obs.obs.dists <- unclass(obs.obs.dists)
# Create variogram with custom distances
data.v <- variog.dist(
datageo,
trend = "cte",
dist.mat = obs.obs.dists,
max.dist = max(obs.obs.dists),
breaks = seq(0, max(obs.obs.dists), length.out = 10))
# Fit variogram
data.vfit <- geoR::variofit(
data.v,
cov.model = "exponential")
# Plot
plot(data.v); lines(data.vfit, col = "blue")``````

We use the variogram to predict on the entire map, also using custom distances with `GeoDist::krige.conv.dist`.

``````# Distances to raster
obs.r.dists <- st_distance(obs_sf, r_sf)
obs.r.dists <- unclass(obs.r.dists)
# Distances between raster locations
r.r.dists <- st_distance(r_sf, r_sf)
r.r.dists <- unclass(r.r.dists)
# Krige with custom distances
r.krige <- krige.conv.dist(
geodata = datageo,
locations = coordinates(r),
krige = krige.control(obj.model = data.vfit),
dist.mat = obs.obs.dists,
loc.dist = obs.r.dists,
loc.loc.dist = r.r.dists)
# Fill in raster with predictions
r.pred <- r
values(r.pred) <- r.krige\$predict
# Plot prediction raster
rasterVis::gplot(r.pred) +
geom_tile(aes(fill = value)) +
scale_fill_gradient2(low = 'red', mid = "yellow",
high = 'green', midpoint = 0.5) +
geom_point(data = obs,
aes(x = lon, y = lat, colour = factor(value))) +
worldmap_predict +
guides(colour = FALSE) +
coord_quickmap()``````

## Fit a `gam` model on 3D-coordinates

If you want to model the global distribution, with a GAM for instance, using latitude and longitude, you will face the discontinuity in the longitude coordinates. In that case, it may be preferable to use 3D-cartesian coordinates of the globe. In `GeoDist`, I included the function `sph2car` to transform geographical coordinates into cartesian. This is a modification of `sphereplot::sph2car`, which allows to include the ellipsoid flattening (`geosphere::refEllipsoids()`).

``````# Transform data in 3D ----
library(ggplot2)
library(maptools)
library(raster)
library(rgl)

# Transform observation in cartesian coords
obs.cart <- data.frame(sph2car(obs[,1:2]), value = obs\$value)
# Approximation with a gam model
library(mgcv)
gam1 <- gam(value ~ te(x, y, z, k = 3),
data = obs.cart,
family = binomial)
# Predict on Raster
r.cart <- data.frame(sph2car(coordinates(r)))
pred <- predict(gam1, r.cart, type = "response")
r.pred <- r
values(r.pred) <- c(pred)

# Plot prediction raster
rasterVis::gplot(r.pred) +
geom_tile(aes(fill = value)) +
scale_fill_gradient2(low = 'red', mid = "yellow",
high = 'green', midpoint = 0.5) +
geom_point(data = obs,
aes(x = lon, y = lat, colour = factor(value))) +
worldmap_predict +
guides(colour = FALSE) +
coord_quickmap()``````

## Plot predictions on a globe with `rgl`

There are different ways to get a projection allowing to show Earth as a globe and to create rotation gif. In my case, I like using `rgl` library because the output in interactive and I can easily travel the world with the mouse.

Projection of our predicted map within a `rgl` 3D plot requires a few steps:

• Build the 3D representation of the prediction raster with triangulation
• [Edit: 2017-11-06] Directly create 3D triangulation using `geometry::convhulln` (Thanks to @mdsumner)
• [removed] Triangulate both halves of the globe separately using library `deldir`
• [removed] Combine both halves
• Create color for each triangle facet according to predictions
• Draw the 3D rgl map
• Drape the globe with predictions using `triangles3d`
• Add observations as `spheres3d`
• Add Earth rotation axis for fun
• Add a map of the world on the 3d sphere
• Retrieve geographical coordinates and transform as cartesian with function `sph2car`
• Use `transform3d` and `rotationMatrix` functions to rotate Earth and take snapshots
• Create a gif
``````# Plot in 3d using rgl ----
library(geometry)
library(dplyr)
library(rgl)

# ---- Edit - 2017-11-06 ----
# Triangulate entire globe directly with geometry ----
# Thanks to Michael Sumner @mdsumner
tri3d <- geometry::convhulln(r.cart)
r.cart.tri <- r.cart[t(tri3d), ] %>%
mutate(n = c(t(tri3d)))

# Edit: Kept here as reminder but simpler with "geometry"
if (FALSE) {
# _With deldir = 2D triangulation ----
library(deldir)
# Triangulate top half of the globe
r.cart.top <- as.data.frame(r.cart) %>% as.tbl() %>%
mutate(n = 1:n()) %>%
filter(z >= 0)
r.cart.top.del <- deldir(as.data.frame(r.cart.top[,1:2]))
r.cart.top.tri <- do.call(rbind, triang.list(r.cart.top.del))
r.cart.top.tri\$n <- r.cart.top\$n[r.cart.top.tri\$ptNum]

# Triangulate bottom half of the globe
r.cart.bottom <- as.data.frame(r.cart) %>% as.tbl() %>%
mutate(n = 1:n()) %>%
filter(z <= 0)
r.cart.bottom.del <- deldir(as.data.frame(r.cart.bottom[,1:2]))
r.cart.bottom.tri <- do.call(rbind, triang.list(r.cart.bottom.del))
r.cart.bottom.tri\$n <- r.cart.bottom\$n[r.cart.bottom.tri\$ptNum]

# Combine top and bottom
r.cart.tri <- rbind(r.cart.top.tri, r.cart.bottom.tri) %>%
mutate(z = r.cart[.\$n,"z"])
}
# ---- End of Edit ----

# Define a vector of colors for predictions
n.break <- 20
colors <- alpha(colorRampPalette(c("red", "yellow", "green"))(n.break), .4)
brk <- seq(0, 1, len = n.break + 1)
pred.col <- colors[as.numeric(as.character(
cut(pred, breaks = brk, include.lowest = TRUE, labels = 1:n.break)))]

# Print in 3d
triangles3d(r.cart.tri\$x,
r.cart.tri\$y,
r.cart.tri\$z,
col = pred.col[r.cart.tri\$n],
alpha = 0.9,
specular = "black")
# Observation with radius in scale of coordinates
spheres3d(obs.cart[,1:3],
color = obs.cart[,4] + 2)
# Black background
rgl.bg(color = c("black"))
segments3d(x = c(0,0),
y = c(0,0),
z = c(min(r.cart.tri\$z)*1.2,
max(r.cart.tri\$z)*1.2),
col = "white",
lwd = 2)

# Map world on the 3d sphere ----
library(rworldmap)
data(countriesCoarse)

for (i in 1:nrow(countriesCoarse)) { # i <- 1
Pols <- countriesCoarse@polygons[[i]]
for (j in 1:length(Pols)) { # j <- 1
lines3d(data.frame(sph2car(
countriesCoarse@polygons[[i]]@Polygons[[j]]@coords
)), col = "white", lwd = 2)
}
}

# Rotate and save as gif
# extraWD: directory where to save img of rotation
nb.img <- 45
angle.rad <- seq(0, 2*pi, length.out = nb.img)
# NorthPole on top, Europe-Africa in front.
uM0 <- rotationMatrix(-pi/2, 1, 0, 0) %>%
transform3d(rotationMatrix(-2, 0, 0, 1)) %>%
transform3d(rotationMatrix(-pi/12, 1, 0, 0))
# Change viewpoint
rgl.viewpoint(theta = 0, phi = 0, fov = 0, zoom = 0.7,
userMatrix = uM0)

for (i in 1:nb.img) {
# Calculate matrix rotation
uMi <- transform3d(uM0, rotationMatrix(-angle.rad[i], 0, 0, 1))
# Change viewpoint
rgl.viewpoint(theta = 0, phi = 0, fov = 0, zoom = 0.7,
userMatrix = uMi)
# Save image
filename <- paste0(extraWD, "/gif/pic", formatC(i, digits = 1, flag = "0"), ".png")
rgl.snapshot(filename)
}
# Create gif
system(glue::glue("convert -delay 10 {extraWD}/gif/*.png -loop 0 ../../static{StaticImgWD}/Globe3D_rgl.gif"))``````

The complete R-script extracted from this post is here on my github.

###### Citation:

``` Rochette Sébastien. (2017, Nov. 01). "Spatial interpolation on Earth as a 3D sphere". Retrieved from https://statnmap.com/2017-11-01-spatial-interpolation-on-earth-as-a-3d-sphere/. ```
``````@misc{Roche2017Spati,