球体表面地理数据的插值
Interpolation of geodata on surface of a sphere
我有一个数据集，其中包含 lat/lon 坐标和每个地理位置对应的 0/1 值（4 到 200+ 个数据点）。现在，我想根据插值结果对空隙进行插值并为地球表面添加颜色。我的主要问题是插值 "around the globe"，因为目前我是在飞机上做的，这显然行不通。
我的数据
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 包直接在 lat/lon 坐标上完成的。这可行，但没有考虑到 lat/lon 坐标位于球体上。
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 中的球体进行插值？
您可以自己计算点与网格之间的距离，然后使用您自己的插值法。例如，下面是您的数据示例的反距离插值。
生成数据
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)
我有一个数据集，其中包含 lat/lon 坐标和每个地理位置对应的 0/1 值（4 到 200+ 个数据点）。现在，我想根据插值结果对空隙进行插值并为地球表面添加颜色。我的主要问题是插值 "around the globe"，因为目前我是在飞机上做的，这显然行不通。
我的数据
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 包直接在 lat/lon 坐标上完成的。这可行，但没有考虑到 lat/lon 坐标位于球体上。
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 中的球体进行插值？
您可以自己计算点与网格之间的距离，然后使用您自己的插值法。例如，下面是您的数据示例的反距离插值。
生成数据
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)