地球表面的地理數據插值

Question

我有一個包含經緯度座標和相應的每個地理位置（4到200+個數據點）的0/1值的數據集。現在，我想插入空隙並根據插值結果向地球表面添加顏色。我所面臨的主要問題是插入「全球各地」，因爲目前我在一架飛機上做，這顯然不起作用。

我的數據

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 

可視化的數據點

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) 

二元插值平面

目前，雙變量樣條插值是使用akima papckage直接在緯度/經度座標上完成的。這可行，但不考慮緯度/經度座標位於球體上。

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) 

顯然，這並不爲球體工作，因爲左側不右側匹配。在一個球體上，插值應該是無縫的。

什麼樣的方法我可以用不插在R A球？

Answer 1

您可以自己計算點和網格之間的距離，然後使用自己的插值。例如，下面是您的數據示例中的反距離插值。

生成數據

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 raster image 
r <- raster(nrows=ny, ncols=nx, crs='+proj=longlat', 
      xmn=-180, xmx=180, ymn=-90, ymx=90) 

計算時座標不被投射在庫sp使用大圓距離點和光柵

功能spDists之間的距離。這意味着兩點之間計算的距離最短。

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

插值使用逆距離內插

在這裏，我建議使用具有功率2（idp=2）經典的逆距離內插來內插一個球體。如果需要其他功率或線性插值，或者如果要插入有限數量的鄰居，則可以修改自己的計算。

# 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 

然後畫出結果

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