 # 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],
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.

###### Citation:

For attribution, please cite this work as:
 Rochette S. (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/. 

BibTex citation:
@misc{Roche2017Spati,
author = {Rochette S},
title = {Spatial interpolation on Earth as a 3D sphere},
url = {https://statnmap.com/2017-11-01-spatial-interpolation-on-earth-as-a-3d-sphere/},
year = {2017}
}