Interpolation de géodonnées à la surface d'une sphère

Question

J'ai un jeu de données avec coordonnées lat/lon et une valeur 0/1 correspondante pour chaque géolocalisation (4 à 200+ points de données). Maintenant, je veux interpoler les vides et ajouter des couleurs à la surface du globe en fonction des résultats de l'interpolation. Le principal problème que j'ai est d'interpoler "autour du globe", car actuellement je fais dans un avion, ce qui évidemment ne fonctionne pas.

Mes données

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) 

Bivariée interpoler dans le plan

Actuellement, l'interpolation spline bivariée est effectuée directement sur les coordonnées lat/lon en utilisant le document akima. Cela fonctionne, mais ne tient pas compte du fait que les coordonnées lat/lon se trouvent sur une sphère.

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) 

Il est évident que cela ne fonctionne pas pour une sphère, comme le côté gauche ne correspond pas au côté droit. Sur une sphère, l'interpolation devrait être transparente.

Quelles approches puis-je utiliser n'interpoler sur une sphère en R?

Answer 1

Vous pouvez calculer vous-même les distances entre les points et la grille, puis utiliser votre propre interpolation. Par exemple, ci-dessous est une interpolation de distance inverse sur votre exemple de données.

générer des données

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 

Créer raster pour l'interpolation de grille

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

distances Calculer entre les points et raster

Fonction spDists dans la bibliothèque sp utiliser la distance orthodromique lorsque les coordonnées ne sont pas projetées . Cela signifie que la distance calculée entre deux points est la plus courte.

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

Interpoler sur une sphère en utilisant l'interpolation de la distance inverse

Ici, je vous propose d'interpoler en utilisant l'interpolation de la distance inverse classique avec une puissance 2 (idp=2). Vous pouvez modifier vous-même le calcul si vous voulez une autre interpolation de puissance ou linéaire, ou si vous voulez interpoler avec un nombre limité de voisins.

# 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 

Ensuite, tracer les résultats

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