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

}```