Monday, May 30, 2011

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

3 comments:

  1. Hi,

    I was wondering if this function would work for plotting an xyz array or does it only work on lists? Basically I use NetCDF data a lot, and all the data I work on is array based. What I would like to do is some stereographic projections but haven't been able to work out how to do them in R. Your plots (and the ones for the MCA) look great and I would really like to be able to do similar!

    ReplyDelete
  2. Hi Alex - Not sure If I completely understand your question. The "xyz_data" argument is actually a matrix with columns x, y, and z. However, the "polys" object is a list - and probably needs to remain so in order to define the coordinates of each polygon to be plotted on the projection. The above example uses a "perspective" projection I believe. The map.xyz function assumes a default of "stereographic". I hope this helps - thanks for your comment!

    ReplyDelete
  3. I have posted another example that more clearly demonstrates this method here: http://menugget.blogspot.com/2012/02/having-recently-received-comment-on.html

    ReplyDelete