## 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

- Drape the globe with predictions using
- Add a map of the world on the 3d sphere
- Retrieve geographical coordinates and transform as cartesian with function
`sph2car`

- Retrieve geographical coordinates and transform as cartesian with function
- 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],
radius = rep(0.2e6, nrow(obs.cart)),
color = obs.cart[,4] + 2)
# Black background
rgl.bg(color = c("black"))
# Add Earth rotation axe
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.