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

## No comments:

## Post a Comment