Wednesday, April 8, 2015

Map projection cheat sheet


Here's a cheat sheet for map projection settings for the mapproject function (mapproj package). I've repeated the same general procedure for each map: 1. create blank world map, 2. project colored polygons, 3. overlay country boundaries. Some specifications were needed for individual map projections (i.e. different orientation, lat/lon limits, additional parameters). Most projections work well as long as the limits and orientation are valid - in a few cases country boundaries ('square' and 'hex') were wrapped incorrectly. I was unable to get the two spheroid projections working due to the escape character in their names ('sp\_mercator', 'sp\_albers') (maybe someone has a suggestion here?).

In addition to the maps and mapproj packages, the sinkr package is also used for the calculation of global distance, value to color conversion, and polygon creation.


Script to reproduce (high resolution png and pdf versions):
# required packages =======================================
# library(devtools); install_github("marchtaylor/sinkr") # sinkr package installation
library(sinkr)
library(maps)
library(mapproj)
 
# make data ===============================================
deg <- 10
x <- seq(-180+deg/2, 180-deg/2, by=10)
y <- seq(-90+deg/2, 90-deg/2, by=10)
grd <- expand.grid(x=x, y=y)
grd$z <- earthDist(40+(38/60), 11+(10/60), grd$x, grd$y) # distance to Hadar, Ethiopia [km]
 
mat <- matrix(grd$z, length(x), length(y))
image(x,y, mat, col=jetPal(100))
POLY <- matrixPoly(x, y, z=mat)
PAL <- colorRampPalette(c(rgb(1,0.5,0.5), rgb(1,1,0.5), rgb(0.5,1,1), rgb(0.5,0.5,1)), bias=2)
COLLEV <- PAL(31)
ZLIM <- c(0, max(mat))
COL <- val2col(mat, zlim=ZLIM, col=COLLEV)
 
# create projection settings ==============================
PROJ <- c(
  "mercator","sinusoidal","cylequalarea","cylindrical","rectangular","gall",
  "mollweide","gilbert","azequidistant","azequalarea","gnomonic",
  "perspective","orthographic","stereographic","laue","fisheye","newyorker",
  "conic","simpleconic","lambert","albers","bonne","polyconic","aitoff","lagrange",
  "bicentric","elliptic","globular","vandergrinten","eisenlohr","guyou","square",
  "tetra","hex","harrison","trapezoidal","lune","mecca","homing"
)
PARAM <- vector(mode="list", length(PROJ))
ORIENT <- lapply(PARAM, FUN=function(x){c(90,0,0)})
XLIM <- lapply(PARAM, FUN=function(x){c(-180, 180)})
YLIM <- lapply(PARAM, FUN=function(x){c(-90, 90)})
 
tmp=which(PROJ=="cylequalarea"); PARAM[[tmp]]=0
tmp=which(PROJ=="rectangular"); PARAM[[tmp]]=0
tmp=which(PROJ=="gall"); PARAM[[tmp]]=0
tmp=which(PROJ=="perspective"); PARAM[[tmp]]=2; ORIENT[[tmp]]=c(45,0,0)
tmp=which(PROJ=="orthographic"); ORIENT[[tmp]]=c(45,0,0)
tmp=which(PROJ=="azequalarea")
tmp=which(PROJ=="gnomonic"); ORIENT[[tmp]]=c(0,0,0); XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="laue"); ORIENT[[tmp]]=c(45,0,0); XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="fisheye"); PARAM[[tmp]]=1.5
tmp=which(PROJ=="newyorker"); PARAM[[tmp]]=10
tmp=which(PROJ=="conic"); PARAM[[tmp]]=45
tmp=which(PROJ=="simpleconic"); PARAM[[tmp]]=c(30,-80)
tmp=which(PROJ=="lambert"); PARAM[[tmp]]=c(-45,45); XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="albers"); PARAM[[tmp]]=c(45,-80)
tmp=which(PROJ=="bonne"); PARAM[[tmp]]=0
tmp=which(PROJ=="bicentric"); PARAM[[tmp]]=0
tmp=which(PROJ=="elliptic"); PARAM[[tmp]]=0
tmp=which(PROJ=="tetra"); ORIENT[[tmp]]=c(45,0,0)
tmp=which(PROJ=="harrison"); PARAM[[tmp]]=c(2, 0); ORIENT[[tmp]]=c(45,0,0)
tmp=which(PROJ=="trapezoidal"); PARAM[[tmp]]=c(0,45); XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="lune"); ORIENT[[tmp]]=c(10,0,0); PARAM[[tmp]]=c(-10,10); XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="mecca"); PARAM[[tmp]]=0; XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
tmp=which(PROJ=="homing"); PARAM[[tmp]]=10; XLIM[[tmp]]=c(-50,50); YLIM[[tmp]]=c(-50,50)
 
 
# map projection testing ==================================
i <- tmp
map("world", proj=PROJ[[i]], par=PARAM[[i]], orient=ORIENT[[i]], type="n", xlim=XLIM[[i]], ylim=YLIM[[i]])
for(j in seq(POLY)){
  polygon(mapproject(POLY[[j]]), border="white", col=COL[j], lwd=0.5)
}
map("world", projection="", add=TRUE, wrap=TRUE)
box()
mtext(PROJ[[i]], line=0.1, side=3)
 
 
# create cheat sheet (run from "pdf(..." down for pdf version ============
png("mapproj_cheatsheet.png", width=16, height=10, units="in", res=200, type="cairo", family="Helvetica")
#pdf("mapproj_cheatsheet.pdf", width=16, height=10, family="Helvetica")
layout(matrix(1:40, 5, 8), widths=rep(1,8), heights=rep(1,5), respect=TRUE)
op <- par(cex=1, ps=10, mar=c(0.5,1,1,0.5))
for(i in seq(length(PROJ))){
  map("world", projection=PROJ[[i]], parameters=PARAM[[i]], orientation=ORIENT[[i]], type="n", xlim=XLIM[[i]], ylim=YLIM[[i]])
  usr <- par("usr")
  rect(usr[1], usr[3], usr[2], usr[4], col="grey95")
  for(j in seq(POLY)){
    polygon(mapproject(POLY[[j]]), border="white", col=COL[j], lwd=0.5)
  }
  map("world", projection="", add=TRUE, wrap=TRUE, col="grey50", lwd=0.2)
  box()
  mtext(PROJ[[i]], line=0.1, side=3, font=2)
  reltext(0,0.2,labels=
            paste(
              paste0("proj=", paste0("'", PROJ[[i]], "'")),
              paste0("par=c(", paste(PARAM[[i]], collapse=","), ")"),
              paste0("orient=c(", paste(ORIENT[[i]], collapse=","), ")"),
              paste0("xlim=c(", paste(XLIM[[i]], collapse=","), ")"),
              paste0("ylim=c(", paste(YLIM[[i]], collapse=","), ")"),
              sep="\n"
            ),
          pos=4, cex=0.8, font=2
  )
  print(i)
}
par(mar=c(1,2,1,7), las=2)
imageScale(c(mat), zlim=ZLIM, col=COLLEV, axis.pos=4)
box()
par(las=0)
mtext("Distance to \nHadar, Ethiopia [km]", side=4, line=4, font=2)
par(op)
dev.off()
Created by Pretty R at inside-R.org



2 comments:

  1. Hi,
    This is really useful and great to see.. but I can't quite get this to work for how I'd like it to work. Basically I've found that increasingly people are reporting there plots on a Mollweide projection and I'd like to know how to do this in R if you have data say from a climate model. From what I can see the majority of solutions involve using the raster package but I was wondering if there is a more generic solution that you can suggest?
    Best,
    Alex

    ReplyDelete
    Replies
    1. Hi Alex, I wish I knew the raster package better, but I haven't invested the time yet. I have usually employed the strategy of producing polygons of my spatial grid, and projecting these on a map with mapproj (e.g. in Ex. 2 here: http://menugget.blogspot.de/2012/04/create-polygons-from-matrix.html).

      Delete