Making fancy relief maps with GBIF records and sampling localities for a rangewide population genomics project

Craig F. Barrett Department of Biology West Virginia University

This code shows you how generate fancy relief maps, gather occurrence data from GBIF, and include a .csv file of latitudes, longitudes, and taxa

e.g. Here is what the CSV file with collections infor looks like:

accession	tax	lon	lat
103a_NM	vree	-105.7481667	32.94283333
103a_NM	vree	-105.7481667	32.94283333
103a_NM	vree	-105.7481667	32.94283333
etc	etc	etc	etc

### Get appropriate libraries
# if needed, install.packages(c("raster","dismo","rgbif","dplyr"))
library("raster") # for downloading relief map data
library("dismo") # for downloading gbif data -- approach 1
library("rgbif") # for downloading gbif data -- approach 2
library("dplyr")   # for subsetting the data by taxon

### Get data for each species with "dismo" package
all<-gbif("Corallorhiza striata")

3226 records found
0-300-600-900-1200-1500-1800-2100-2400-2700-3000-3226 records downloaded

str<-gbif("Corallorhiza striata var. striata")

632 records found
0-300-600-632 records downloaded

vree<-gbif("Corallorhiza striata var. vreelandii")

235 records found
0-235 records downloaded

inv<-gbif("Corallorhiza striata var. involuta")

60 records found
0-60 records downloaded

ben<-gbif("Corallorhiza bentleyi")

34 records found
0-34 records downloaded

# Get altitude, spole, and aspect layers
alt.usa <- getData('alt', country = c("USA"), level = 1)

trying URL 'https://biogeo.ucdavis.edu/data/diva/msk_alt/USA_msk_alt.zip'
Content type 'application/zip' length 20073726 bytes (19.1 MB)
downloaded 19.1 MB

Discarded datum Unknown based on WGS84 ellipsoid in Proj4 definitionDiscarded datum Unknown based on WGS84 ellipsoid in Proj4 definitionDiscarded datum Unknown based on WGS84 ellipsoid in Proj4 definitionDiscarded datum Unknown based on WGS84 ellipsoid in Proj4 definitionreturning a list of RasterLayer objects

slope.usa <- terrain(alt.usa[[1]], opt='slope') # only first layer (country level)
aspect.usa <- terrain(alt.usa[[1]], opt='aspect')
alt.mex <- getData('alt', country = c("Mexico"), level = 1)

Discarded datum Unknown based on WGS84 ellipsoid in Proj4 definition

slope.mex <- terrain(alt.mex[[1]], opt='slope') # only first layer (country level)
aspect.mex <- terrain(alt.mex[[1]], opt='aspect')
alt.can <- getData('alt', country = c("Canada"), level = 1)

Discarded datum Unknown based on WGS84 ellipsoid in Proj4 definition

slope.can <- terrain(alt.can[[1]], opt='slope') # only first layer (country level)
aspect.can <- terrain(alt.can[[1]], opt='aspect')

# combine layers
hill.usa <- hillShade(slope.usa, aspect.usa, angle = 40, direction = 270)
hill.can <- hillShade(slope.can, aspect.can, angle = 40, direction = 270)
hill.mex <- hillShade(slope.mex, aspect.mex, angle = 40, direction = 270)
# Combine USA, Canada, Mexico data
hill.na<-merge(hill.usa, hill.can, hill.mex)

## Read in OUR data for collection localities
samples<-read.csv("latlon.csv")
## Use dplyr to subset based on taxon
stri<-samples %>% filter(tax == "stri")
vree<-samples %>% filter(tax == "vree")
snv<-samples %>% filter(tax == "snv")
cra<-samples %>% filter(tax == "cra")

## plot the points
# Plot the base map
plot(hill.na, col=grey(0:100/100), xlab="Latitude", ylab="Longitude", cex=2.0, legend=FALSE, xlim = c(-135, -50), ylim = c(15, 60)) +
## for the C striata GBIF records
points(all$lon, all$lat, pch = 21, col = "black", cex = 1.5)

integer(0)

## For our collections
points(stri$lon, stri$lat, pch = 19, col = "red", cex = 1.5)

points(vree$lon, vree$lat, pch = 19, col = "green", cex = 1.5)
points(snv$lon, snv$lat, pch = 19, col = "blue", cex = 1.5)

points(cra$lon, cra$lat, pch = 19, col = "magenta", cex = 1.5)

### add custom legend ####
legend("bottomright", legend=c("Corallorhiza striata, no variety listed", "Corallorhiza striata var. striata", "Corallorhiza striata var. vreelandii", "Corallorhiza striata Sierra Nevada", "Corallorhiza striata Coast Ranges/Cascades"), col=c("black", "red", "green", "blue", "magenta"), pch=16, cex=0.8, pt.cex=2.0, inset=0.02)

Error in strwidth(legend, units = "user", cex = cex, font = text.font) : 
  plot.new has not been called yet