A ridiculous proof of concept: xyz interpolation

Ridiculous Orb

This is really the last one on this theme for a while... I had alluded to a combination of methods regarding xyz interpolation at the end of my last post and wanted to demonstrate this in a final example.

The ridiculousness that you see above involved two interpolation steps. First, a thin plate spline interpolation ("Tps" function of the fields package) is applied to the original random xyz field of distance to Mecca. This fitted model is then used to predict values at a new grid of 2° resolution. Finally, in order to avoid plotting polygons for each grid (which can be slow for fine grids), I obtain their projected coordinates with the mapproject function. Using these projected coordinates and their respective z values, a second interpolation is done with the "interp" function of the akima package onto a relatively fine grid of 1000x1000 positions. The result is a smooth field that can then be overlayed on the map using the "image" function (very fast).

So you may ask - When is this even necessary? I would say that it really only makes sense for projecting a filled.contour-type plot for relatively sparse geographic data. Be warned - for large amounts of xyz data, the interpolation algorithms can take a long time.

A couple of functions, found within this blog, are needed to reproduce the plot (earth.dist, color.palette).

the code to reproduce the figure...

#required packages
require(akima)
require(maps)
require(mapproj)
require(fields)
 
#required functions (from "www.menugget.blogspot.com")
source("earth.dist.R")
source("color.palette.R")
 
#create xyz data - distance to Mecca
set.seed(1)
mecca <- data.frame(lon=39.826138, lat=21.422508)
x <- runif(1000, min=-180, max=180)
y <- runif(1000, min=-90, max=90)
z <- earth.dist(x, y, mecca$lon, mecca$lat)
 
#Method 4. Thin-spline fit onto new grid and smooth interpolate
coord <- data.frame(lon=x,lat=y)
# a useful function from field.pdf manual:
dircos<- function(x1){
coslat1 <- cos((x1[, 2] * pi)/180)
sinlat1 <- sin((x1[, 2] * pi)/180)
coslon1 <- cos((x1[, 1] * pi)/180)
sinlon1 <- sin((x1[, 1] * pi)/180)
cbind(coslon1*coslat1, sinlon1*coslat1, sinlat1)}
# fit in 3-d to direction cosines
fit <- Tps(dircos(coord),z)
#create grid and polygons for projection
grd.res <- 2 #grid resolution (degrees)
xcenter <- seq(-180+grd.res/2, 180-grd.res/2, by=grd.res)
ycenter <- seq(-90+grd.res/2, 90-grd.res/2, by=grd.res)
tmp <- list(lon=xcenter, lat=ycenter)
grd <- make.surface.grid(tmp)
#predict fitted model to grid
pred <- predict(fit, dircos(grd))
 
#Plot device settings
HEIGHTS <- c(1,4,1)
WIDTHS <- c(1,4,1)
OMI <- c(0,0,0,0)
#open a device
png("weird_interpolation.png", width=sum(WIDTHS)+sum(OMI[c(2,4)]), height=sum(HEIGHTS)+sum(OMI[c(1,3)]), units="in", res=400)
layout(matrix(c(1,1,1,1,2,1,1,1,1),nrow=3,ncol=3), width=WIDTHS, height=HEIGHTS, respect=TRUE)
layout.show(2)
par(mai=c(0,0,0,0), bg=rgb(0,0,0.1), omi=OMI)
#plot some stars
plot(1, ylim=c(-1,1), xlim=c(-1,1), xaxt="n", yaxt="n", t="n", xaxs="i", yaxs="i")
xstar <- runif(1000, min=-1, max=1)
ystar <- runif(1000, min=-1, max=1)
cexstar <- rnorm(1000, mean=0.3, sd=0.2)
points(xstar, ystar, pch=24, cex=cexstar, col="white", bg="white", lwd=0)
points(xstar, ystar, pch=25, cex=cexstar, col="white", bg="white", lwd=0)
#map setup
breaks <- pretty(z, n=1000)
XLIM=c(-180,180)
YLIM=c(-90,90)
PROJ="orthographic"
ORIENT=c(60,0,0)
PAR=NULL
pal <- color.palette(rev(c("purple","black","yellow")), n.steps.between=c(2,2)) 
#plot blank map
map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA)
#project grid
p <- mapproject(x=grd[,1], y=grd[,2])
#select grid points that mapped correctly
incl <- which(!is.na(p$x))
#build xyz object (projected coordinates)
pxyz <- data.frame(x=p$x[incl], y=p$y[incl],  z=pred[incl])
#interpolate between projected coordinates (fine grid 1000x1000)
INTERP <- interp(pxyz$x, pxyz$y, pxyz$z, 
 xo=seq(min(pxyz$x), max(pxyz$x), length = 1000),
 yo=seq(min(pxyz$y), max(pxyz$y), length = 1000)
)
#overlay image of interpolation
image(INTERP, col=pal(length(breaks)-1), breaks=breaks, add=TRUE)
#turn off device
dev.off()

Created by Pretty R at inside-R.org