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

Hi,

ReplyDeleteI 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!

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!

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

ReplyDelete