球面上的地理数据插值 [英] Interpolation of geodata on surface of a sphere

问题描述

我有一个包含经/纬度坐标的数据集，并且每个地理位置都有一个对应的0/1值(4至200个以上的数据点).现在，我想对空隙进行插值，并根据插值结果向地球表面添加颜色.我遇到的主要问题是对全球"进行插值，因为当前我是在飞机上进行的，这显然行不通.

I have a dataset with lat/lon coordinates and a corresponding 0/1 value for each geolocation (4 to 200+ datapoints). Now, I want to interpolate the voids and add colors to the surface of the globe based on the interpolation results. The main problem I have is to interpolate "around the globe", because currently I do in a plane, which obviously does not work.

我的数据

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填充直接在纬度/经度坐标上完成二元样条插值.这是可行的，但没有考虑纬度/经度坐标位于球体上.

Currently, the bivariate spline interpolation is done directly on the lat/lon coordinates using the akima papckage. This works, but does not take into account that the lat/lon coordinates lie on a sphere.

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)

显然，这不适用于球体，因为左侧与右侧不匹配.在球体上，插值应该是无缝的.

Obviously this does not work for a sphere, as the left side does not match the right side. On a sphere,the interpolation should be seamless.

我可以使用哪些方法对R中的球面进行插值?

推荐答案

您可以自己计算点与网格之间的距离，然后使用自己的插值.例如，下面是数据示例中的反距离插值.

You can calculate distances between points and grid yourself and then use your own interpolation. For instance, below is a inverse distance interpolation on your data example.

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使用大圆距.这意味着两点之间的计算距离最短.

Calculate distances between points and raster

Function spDists in library sp use the Great Circle Distance when coordinates are not projected. This means that the distance calculated between two points is the shortest.

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

使用反距离插值在球上插值

在这里，我建议使用功率为2(idp=2)的经典逆距离插值进行插值.如果要使用其他幂或线性插值，或者要对有限数量的邻居进行插值，则可以修改自己的计算.

Interpolate on a sphere using inverse distance interpolation

Here I propose to interpolate using the classical inverse distance interpolation with power 2 (idp=2). You can modify yourself the calculation if you want other power or linear interpolation, or if you want to interpolate with a limited number of neighbors.

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

这篇关于球面上的地理数据插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！