地球表面的地理数据插值

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)