map.xyz(): interpolation of XYZ data and projection onto a map

     I am still struggling to get a grasp of R's mapping capabilities. Part of my frustration lies in the fact that I often work on areas near the poles, which complicates interpolation across the 180 degree line. For smaller areas, interpolation can be done using the interp() function in the package akima. I have taken the results from interp and projected the image onto a map. You will need the akima, maps, and mapproj packages and the functions new.lon.lat(), earth.dist(), and pos2coord().


As an example I have mapped the distance from Mecca, Saudi Arabia:

The function...

#interpolation of xyz data and projection onto a map
map.xyz <- function(xyz_data, file_name="map.xyz.png",
 
projection="stereographic", orientation=NULL,
 
par=NULL, fill=FALSE, nside=40, breaks=100, res=100,
 
xlim = NULL, ylim = NULL
 
){
 
 
 if(length(which(rowSums(is.na(xyz_data)) != 0)) > 0){
 
  xyz_data <- xyz_data[-which(rowSums(is.na(xyz_data)) != 0),]
 
 }
 
 
 require(akima)
 
 require(maps)
 
 require(mapproj)
 
 
 
 temp<-interp(
 
  xyz_data[,1], xyz_data[,2], xyz_data[,3],  
 
  xo=seq(min(xyz_data[,1]), max(xyz_data[,1]), length = nside), 
 
  yo=seq(min(xyz_data[,2]), max(xyz_data[,2]), length = nside, extrap=TRUE)
 
 )
 
 
 
 polys<-vector("list", length(as.vector(temp$z)))
 
 for(i in 1:length(polys)){
 
  lonx <- pos2coord(pos=i, dim.mat=dim(temp$z))[1]
 
  laty <- pos2coord(pos=i, dim.mat=dim(temp$z))[2]
 
  ifelse(laty < length(temp$y), neigh_y<-c(laty+1,lonx), neigh_y<-c(laty-1,lonx))
 
  ifelse(lonx < length(temp$x), neigh_x<-c(laty,lonx+1), neigh_x<-c(laty,lonx-1)) 
 
  dist_y <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_y[2]], temp$y[neigh_y[1]])
 
  dist_x <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_x[2]], temp$y[neigh_x[1]])
 
  s1 = new.lon.lat(temp$x[lonx], temp$y[laty], 180, dist_y/2)
 
  s3 = new.lon.lat(temp$x[lonx], temp$y[laty], 0, dist_y/2)
 
  s2 = new.lon.lat(temp$x[lonx], temp$y[laty], 270, dist_x/2)
 
  s4 = new.lon.lat(temp$x[lonx], temp$y[laty], 90, dist_x/2)
 
  polys[[i]] = cbind(c(s2[1], s2[1], s4[1], s4[1]), c(s1[2], s3[2], s3[2], s1[2]))
 
 }
 
 
 
 M.range=range(temp$z, na.rm=TRUE)
 
 M.breaks <- pretty(M.range, n=breaks)
 
 M.cols=colorRampPalette(c("blue","cyan", "yellow", "red"),space="rgb")
 
 colorlut <- M.cols(length(M.breaks)) # color lookup table
 
 colorvalues <- colorlut[((as.vector(temp$z)-M.breaks[1])/(range(M.breaks)[2]-range(M.breaks)[1])*length(M.breaks))+1] # assign colors to heights for each point
 
 if(is.null(xlim)){xlim <- range(xyz_data[,1], na.rm=TRUE)}
 
 if(is.null(ylim)){ylim <- range(xyz_data[,2], na.rm=TRUE)}
 
 png(filename=file_name, res=res)
 
 map("world",projection=projection, orientation=orientation, par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 for(i in 1:length(polys)){
 
  polygon(mapproject(x=polys[[i]][,1], y=polys[[i]][,2]), col=colorvalues[i], border=NA)
 
  print(i)
 
 }
 
 map("world", add=TRUE, projection="", par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 map.grid(c(-180, 180, -80, 80), nx=10, ny=18, labels=FALSE, col="grey")
 
 dev.off()
 
}