Having recently received a comment on a post regarding geographic xyz data interpolation, I decided to return to my original "xyz.map" function and open it up for easier interpretation. This should make the method easier to adapt and follow.

The above graph shows the distance to Mecca as interpolated from 1000 randomly generated lat/lon data using the "interp" function of the akima package. Several functions, found within this blog, are needed to reproduce the plot (pos2coord, earth.dist, new.lon.lat, color.palette, val2col, image.scale). One thing you will notice is the strip of uninterpolated area within the stereographic projection. This is a problem that I have yet to resolve and has to do with the fact that the interpolation is not considering the connection along the 180° longitude line. This will probably require some other type of interpolation based on geographic distances rather than Cartesian coordinates.

**R code to produce the above graph...**

#required packages require(akima) require(maps) require(mapproj) #required functions (from "www.menugget.blogspot.com") source("pos2coord.R") source("earth.dist.R") source("new.lon.lat.R") source("color.palette.R") source("val2col.R") source("image.scale.R") ###create xyz data set.seed(1) mecca <- data.frame(lon=39.826138, lat=21.422508) x <- runif(10000, min=-180, max=180) y <- runif(10000, min=-90, max=90) z <- earth.dist(x, y, mecca$lon, mecca$lat) xyz_data <- cbind(x,y,z) ###interpolate data INTERP <- interp(x, y, z, xo=seq(-180, 180, length = 100), yo=seq(-90, 90, length = 100)) ###basic plot image(INTERP, col=rainbow(100)) map("world",add=TRUE, col=1) ###projection plots #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]] = cbind(c(s2[1], s2[1], s4[1], s4[1]), c(s1[2], s3[2], s3[2], s1[2])) } #values to colors pal <- color.palette(rev(c("purple","black","red","yellow")), n.steps.between=c(2,5,20)) COL <- val2col(INTERP$z, col=pal(100)) #plot parameters WIDTHS=c(4,4,1) HEIGHTS=c(4) OMI=c(0.1,0.1,0.1,0.1) #define plot device and layout png("xyz_2_projections.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:3),nrow=1,ncol=3), width=WIDTHS, height=HEIGHTS, respect=TRUE) layout.show(3) #1st plot - perspective projection XLIM=c(-180,180) YLIM=c(-90,90) PROJ="perspective" ORIENT=c(mecca$lat,mecca$lon,0) PAR=2 par(mai=c(0.1,0.1,0.5,0.1)) map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA) for(i in 1:length(polys)){ polygon(mapproject(x=polys[[i]][,1], y=polys[[i]][,2]), 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("Perspective projection", side=3, line=1) #2nd plot - stereographic projection XLIM=c(-180,180) YLIM=c(-90,90) PROJ="stereographic" ORIENT=c(mecca$lat,mecca$lon,0) PAR=NULL par(mai=c(0.1,0.1,0.5,0.1)) map("world", projection=PROJ, orientation=ORIENT, par=PAR, xlim=XLIM, ylim=YLIM, col=NA) for(i in 1:length(polys)){ polygon(mapproject(x=polys[[i]][,1], y=polys[[i]][,2]), 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("Stereographic projection", side=3, line=1) #3rd plot - color scale par(mai=c(0.1,0.1,0.5,0.7)) image.scale(INTERP$z, col=pal(100), horiz=FALSE, xaxt="n", yaxt="n", xlab="", ylab="") box() axis(4) mtext("Distance to Mecca [km]", side=4, line=3) #close device dev.off()

The projection function failed for some input data.

ReplyDeleteYes, I think for some projections, the country polygons do not always plot perfectly. However, in this example, you do not really see a dramatic effect. I have come across some tricks to avoid these problems in other blogs - will try to address this at some point in the future.

ReplyDeleteI obtained the right result! Thank you very much. :)

ReplyDeleteI like the hole at the 180th meridian. Yes, it is a data artifact, but does provide a nice guide to where things are on the map.

ReplyDeleteWell, there's that ... :-)

ReplyDeleteYou can also add labels to the grid lines in the map.grid function.

in:

ReplyDelete#define plot device and layout

the first code line(png)failed,why?

I obtained this error:

ReplyDeleteIn map("world", projection = PROJ, orientation = ORIENT, par = PAR, :

projection failed for some data

waht could be the reason?.Thanks

Hi ANGEL666 - As I mentioned in an earlier comment, some country polygons might be cut off or not plotted correctly by some projections. You don't see too much of a noticeable affect here though. I'm not sure why your "png" line failed. Was there a warning? Maybe, make sure that you have closed all open devices before starting the plot next time (by running "dev.off()" until no devices remain). - Marc

ReplyDelete