XYZ geographic data interpolation, part 3

This will be probably be a final posting on interpolation of xyz data as I believe I have come to some conclusions to my original issues. I show three methods of xyz interpolation:
1. The quick and dirty method of interpolating projected xyz points (bi-linear)
2. Interpolation using Cartesian coordinates (bi-linear)
3. Interpolation using spherical coordinates and geographic distances (thin plate spline)


The 1st map shows the "quick and dirty" method. First xy coordinates are converted to coordinates of the projected map in the device . For xy coordinates that are visible, xyz data is interpolated with the interp function. By interpolating to a much finer grid, as defined by the interp arguments "xo" and "yo" of the akima package ("bi-linear" interpolation), you obtain a nice field by which to cover the map. Finally, contour lines and country borders are placed on top. Interestingly, this method deals a work-around to a problem that several folks have asked about in the blogosphere - Can the results of the function "filled.contour" be mapped to a different projection.

The 2nd map uses the same method from my previous post, yet I have now changed how I defined the arguments "xo" and "yo" in the "interp" function ("bi-linear" interpolation) . I had previously used the x and y limits of the new grid sides for the interpolation, which resulted in a missing strip of interpolated grids.  The correct method is to use the coordinates of the grid centers, which fixed most of the problem - only some grids near the pole remain uninterpolated due to lower coverage of the original xyz data.

The 3rd map illustrates the slowest of the three interpolations - yet the most accurate in that it considers geographic distances ("great circle" distances). Here I utilize a "thin plate spline" interpolation, from the fields package, that takes into account the geographic distance between xy coordinates and allows for extrapolation to defined grid centers. As recommended by the authors, interpolations across large geographic domains should first convert xy coordinates to three dimensions using their direction cosines (see page 150 of the fields reference manual). Their function for this transformation, "dircos", is included in the code below.

As a final note, the first method could probably be used to create a smoothed representation of the interpolated values from the other methods (i.e. first project the interpolated values, and then interpolate a finer field by which to cover the map) (see example here).


the code to produce the above figure...

#required packages
require(akima)
require(maps)
require(mapproj)
require(fields)
 
#required functions (from "www.menugget.blogspot.com")
source(paste(fun_path, "pos2coord.R", sep=""))
source(paste(fun_path, "earth.dist.R", sep=""))
source(paste(fun_path, "new.lon.lat.R", sep=""))
source(paste(fun_path, "color.palette.R", sep=""))
source(paste(fun_path, "val2col.R", sep=""))
source(paste(fun_path, "image.scale.R", sep=""))
 
#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)
 
#Plot device settings
HEIGHTS <- c(4)
WIDTHS <- c(4,4,4,1)
OMI <- c(0.1,0.1,1.5,1.5)
 
#open a device
png("interpolation_comparison.png", width=sum(WIDTHS)+sum(OMI[c(2,4)]), height=sum(HEIGHTS)+sum(OMI[c(1,3)]), units="in", res=400)
#x11(width=sum(WIDTHS)+sum(OMI[c(2,4)]), height=sum(HEIGHTS)+sum(OMI[c(1,3)]))
layout(matrix(c(1:4),nrow=1,ncol=4), width=WIDTHS, height=HEIGHTS, respect=TRUE)
layout.show(4)
 
###Interpolations
#map definition
par(mai=c(0.1,0.1,0.1,0.1))
breaks <- pretty(z, n=20)
XLIM=c(-180,180)
YLIM=c(-90,90)
PROJ="orthographic"
ORIENT=c(60,0,0)
PAR=NULL
pal <- color.palette(rev(c("purple","black","red","yellow")), n.steps.between=c(2,5,20)) 
 
#Plot 1 - Method 1 - project xyz data and interpolate
map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA)
xy <- data.frame(x=x,y=y)
p <- mapproject(xy)
incl <- which(!is.na(p$x))
pxyz <- data.frame(x=p$x[incl], y=p$y[incl],  z=z[incl])
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)
)
image(INTERP, col=pal(length(breaks)-1), breaks=breaks, add=TRUE)
map("world", add=TRUE, projection="", xlim=XLIM, ylim=YLIM, col="grey")
contour(INTERP, levels=breaks, add=TRUE, col="white", lwd=1.5)
mtext("Project and interpolate - \n'quick and dirty'", side=3, line=1)
 
#Plot 2 - Method 2 - Interpolate onto cartesian grid
#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)
#interpolate to new grid
INTERP <- interp(x, y, z, xo=xcenter, yo=ycenter)
#make polygons
polys<-vector("list", length(as.vector(INTERP$z)))
for(i in 1:length(polys)){
 lonx <- pos2coord(pos=i, dim.mat=dim(INTERP$z))[1]
 laty <- pos2coord(pos=i, dim.mat=dim(INTERP$z))[2]
 ifelse(laty < length(INTERP$y), neigh_y<-c(laty+1,lonx), neigh_y<-c(laty-1,lonx))
 ifelse(lonx < length(INTERP$x), neigh_x<-c(laty,lonx+1), neigh_x<-c(laty,lonx-1)) 
 dist_y <- earth.dist(INTERP$x[lonx], INTERP$y[laty], INTERP$x[neigh_y[2]], INTERP$y[neigh_y[1]])
 dist_x <- earth.dist(INTERP$x[lonx], INTERP$y[laty], INTERP$x[neigh_x[2]], INTERP$y[neigh_x[1]])
 s1 = new.lon.lat(INTERP$x[lonx], INTERP$y[laty], 180, dist_y/2)
 s3 = new.lon.lat(INTERP$x[lonx], INTERP$y[laty], 0, dist_y/2)
 s2 = new.lon.lat(INTERP$x[lonx], INTERP$y[laty], 270, dist_x/2)
 s4 = new.lon.lat(INTERP$x[lonx], INTERP$y[laty], 90, dist_x/2)
 polys[[i]] = list(x=c(s2[1], s2[1], s4[1], s4[1]), y=c(s1[2], s3[2], s3[2], s1[2]))
}
#color levels
COL <- val2col(INTERP$z, col=pal(length(breaks)-1), breaks=breaks)
#plot and project polygons
map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA)
for(i in 1:length(polys)){
 polygon(mapproject(polys[[i]]), col=COL[i], border=COL[i], lwd=0.1)
}
map("world", add=TRUE, projection="", xlim=XLIM, ylim=YLIM, col="grey")
map.grid(c(-180, 180, -80, 80), nx=12, ny=6, labels=FALSE, col="grey")
mtext("Bi-linear interpolation\nonto Cartesian grid", side=3, line=1)
 
#Method 3. Thin-spline fit onto new grid
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)
pred<- predict(fit, dircos(grd))
#make polygons for grid
polys <- vector(mode="list", nrow(grd))
for(i in seq(polys)){
 xci <- grd[i,1]
 yci <- grd[i,2]
 xs <- c(xci-grd.res/2, xci-grd.res/2, xci+grd.res/2, xci+grd.res/2)
 ys <- c(yci-grd.res/2, yci+grd.res/2, yci+grd.res/2, yci-grd.res/2)
 polys[[i]]$x <- xs
 polys[[i]]$y <- ys
}
#color levels
COL <- val2col(pred, col=pal(length(breaks)-1), breaks=breaks)
 #plot and project polygons
map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA)
for(i in 1:length(polys)){
 polygon(mapproject(polys[[i]]), col=COL[i], border=COL[i], lwd=0.1)
}
map("world", add=TRUE, projection="", xlim=XLIM, ylim=YLIM, col="grey")
map.grid(c(-180, 180, -80, 80), nx=12, ny=6, labels=FALSE, col="grey")
mtext("Thin plate spline interpolation\nusing great circle distances", side=3, line=1)
 
### Color scale
#Plot 4 - color scale
par(mai=c(0.1,0.1,0.1,0.7))
image.scale(1, col=pal(length(breaks)-1), breaks=breaks, horiz=FALSE, xaxt="n", yaxt="n", xlab="", ylab="")
box()
axis(4)
mtext("Distance to Mecca [km]", side=4, line=3)
 
#close device
dev.off()

Created by Pretty R at inside-R.org