Lambert azimuthal equal area projection of the northern hemisphere, with species range color coded by country, and points for collection locations. See http://wp.me/p1Ye5e-a5 for relevant blog post.

ModernMapEurasia.R

#ModernMapEurasia.r
#Kathryn Turner, June 17, 2013
#With extensive help from Andy South

#to create a reprojected map of the Northern Hemisphere
#color code countries for range
#and add collection sites
#in a Lambert azimuthal equal area projection, with the North Pole in the center 

library(rgdal) # Commands for reprojecting the vector data.
library(rworldmap) # Recently updated mapping program.
library(rworldxtra) # Add-ons for rworldmap.

###set output
pdf("collectionMap_color.pdf", useDingbats=FALSE)

#Lambert azimuthal equal area projection
#from http://gis.stackexchange.com/questions/30054/defining-datum-for-lambert-azimuthal-equal-area-and-converting-to-geographic-coo

###the view of the map is centered on the North pole, but slightly offset
#using the exact location causes an error
projectionCRS <- CRS("+proj=laea +lon_0=0.001 +lat_0=89.999 +ellps=sphere") 
#the ellps 'sphere' has a radius of 6370997.0m
par(mai=c(0,0,0.2,0))
#To determines the width of the margins for the map; the "i" designations make the map go to the edge of the plot window.
#par(mai=c(0,0,0.2,0),xaxs="i",yaxs="i")
#including Antarctica causes an error, so I hope you don't need that
sPDF <- getMap()[-which(getMap()$ADMIN=='Antarctica')] 
sPDF <- spTransform(sPDF, CRS=projectionCRS)

###use setLims to reproject xlim and ylim, so that we can zoom in on the map
setLims <- TRUE #FALSE back to whole world
#setLims <- FALSE
if ( !setLims )
{
  xlim <- ylim <- NA
} else
{
  ### TRY FIDDLING WITH THESE LIMITS###
  #how this zooming works in non-intuitive, I suggest trial and error to find your numbers
  xlimUnproj <- c(-52,120)
  ylimUnproj <- c(10,30)
  sPointsLims <- data.frame(x=xlimUnproj, y=ylimUnproj)
  coordinates(sPointsLims) = c("x", "y")
  #set CRS (coordinate reference system) for the points
  #assuming WGS84
  proj4string(sPointsLims) <- CRS("+proj=longlat +ellps=WGS84")
  sPointsLims <- spTransform(sPointsLims, CRS=projectionCRS)
  xlim <- coordinates(sPointsLims)[,"x"]
  ylim <- coordinates(sPointsLims)[,"y"]  
}

###to color code countries
#to view a list of country name formats
sPDF$ADMIN
#setup a color code column filled with 1's
sPDF$colCode <- 1
#set codes for specified countries
sPDF$colCode[ which(sPDF$ADMIN %in% c("Canada","United States of America"))] <- 2
sPDF$colCode[ which(sPDF$ADMIN %in% c("Armenia","Azerbaijan", "Bulgaria", "Georgia", 
                                      "Greece", "Moldova", "Romania","Russia", "Turkey",
                                      "Ukraine", "Serbia"))] <- 3
sPDF$colCode[ which(sPDF$ADMIN %in% c("Poland", "Belarus", "Italy", "Syria", "Czech Republic",
                                      "Estonia", "Switzerland","Latvia","Lithuania", 
                                      "Slovenia", "Serbia","Austria","Belgium", "France",
                                      "Germany","Hungary","Luxembourg","Norway","Slovakia",
                                      "Spain", "United Kingdom", "Kazakhstan", "Turkmenistan", "China"))] <- 4

#create a color palette - note for each value not for each country
colourPalette <- c("lightgray","#F8766D","#00BFC4", "cadetblue1")
#colourPalette can be set to either a vector of colors or one of :white2Black black2White palette heat topo terrain rainbow negpos8 negpos9


mapCountryData(sPDF, nameColumnToPlot="colCode", mapTitle='Global range of C. diffusa', 
               colourPalette=colourPalette, borderCol ='gray24', addLegend = FALSE,
               xlim=xlim, ylim=ylim, catMethod=c(0,1,2,3,4))
#note that catMethod defines the breaks and values go in a category if they are <= upper end
#I specified 4 different colors, so catMethod = c(0,1,2,3,4)

###to plot collection locations from GPS coordinates
> head(popInv)
Latitude Longitude   Pop pch
1 34.85570 -110.2360 US024   1
2 39.22427 -103.1223 US023  19
3 39.69667 -104.8076 US007   1
4 39.70283 -105.3242 US006   1
5 40.12227 -101.2803 US014  19
6 40.37191 -104.4726 US026  19
#I have three different types of collections, specified by point type (pch column in dataframe of coordinates)

#coordinates must be projected to show up on map properly
#invasive populations
popInv <- read.csv("InvPopCoord.csv")
coordinates(popInv) = c("Longitude", "Latitude")
proj4string(popInv) <- CRS("+proj=longlat +ellps=WGS84")
sPointsDF <- spTransform(popInv, CRS=projectionCRS)
points(sPointsDF, pch=popInv$pch, cex=1)
#native populations
popNat <- read.csv("NatPopCoord.csv")
coordinates(popNat) = c("Longitude", "Latitude")
proj4string(popNat) <- CRS("+proj=longlat +ellps=WGS84")
sPointsDFNat <- spTransform(popNat, CRS=projectionCRS)
points(sPointsDFNat, pch=popNat$pch, cex=1) #pch2 for triangles

###Adding things to map
#to plot latitude lines on to map
llgridlines(sPDF, easts=c(-90,-180,0,90,180), norths=seq(0,90,by=15), 
            plotLabels=FALSE, ndiscr=1000) 
#ndiscr=number points in lines, more makes a smoother curved line

#to plot latitude and longitude labels on the map, again, coordinates must be projected
markings <- data.frame(Latitude=as.numeric(c(75,60,45,30,15,85,85)), Longitude=as.numeric(c(-45,-45,-45,-45,-45,0,180)),name=c("75", "60","45","30","15","0","180"))
coordinates(markings) = c("Longitude", "Latitude")
proj4string(markings) <- CRS("+proj=longlat +ellps=WGS84")
sPointsDFmark <- spTransform(markings, CRS=projectionCRS)
text(sPointsDFmark, labels = sPointsDFmark$name, cex=1)

#to plot the pole
# pole <- data.frame(x=0, y=90)
# coordinates(pole) = c("x", "y")
# proj4string(pole) <- CRS("+proj=longlat +ellps=WGS84")
# pole <- spTransform(pole, CRS=projectionCRS)
# points(pole, pch=8, cex=2, lwd=2)

#legends for point type and color code
legend("bottomleft", c("Germination trial", "Broad CG","Maternal effects CG"), 
       pch=c(1,19,15),  bg="white", title = "Sampled populations")
legend("topright", c("Invasive", "Native","Present, status unknown"), fill=c("#F8766D","#00BFC4", "cadetblue1"),
       title="Origin", bg="white")

#shameless plug !
mtext("map made using rworldmap", line=-1, side=1, adj=1, cex=0.6)


dev.off()