i have dataset lat/lon coordinates , corresponding 0/1 value each geolocation (4 200+ datapoints). now, want interpolate voids , add colors surface of globe based on interpolation results. main problem have interpolate "around globe", because in plane, not work.
my 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 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, = 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) 
bivariate interpolate in plane
currently, bivariate spline interpolation done directly on lat/lon coordinates using akima papckage. works, not take account lat/lon coordinates lie on 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 out of range interpolations values z[z < 0] <- 0 z[z > 1] <- 1 grd <- expand.grid(lon = seq(-180,180, = 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 not work sphere, left side not match right side. on sphere,the interpolation should seamless.
what approaches can use interpolate on sphere in r?
you can calculate distances between points , grid , use own interpolation. instance, below inverse distance interpolation on data example.
generate data
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 grid interpolation
## create raster image r <- raster(nrows=ny, ncols=nx, crs='+proj=longlat', xmn=-180, xmx=180, ymn=-90, ymx=90) calculate distances between points , raster
function spdists in library sp use great circle distance when coordinates not projected. means distance calculated between 2 points shortest.
# distance between points , raster s.r.dists <- spdists(x = coordinates(s), y = coordinates(r), longlat = true) interpolate on sphere using inverse distance interpolation
here propose interpolate using classical inverse distance interpolation power 2 (idp=2). can modify calculation if want other power or linear interpolation, or if want interpolate 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 then plot results
# 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) 
Comments
Post a Comment