This document is a work in progress – give us feedback

1 Before we start

1.1 Install R Packages

raster install.packages("raster")
sp install.packages("sp")
RColorBrewer install.packages("RColorBrewer")
tidyverse install.packages("tidyverse")
maptools install.packages("maptools")
dismo install.packages("dismo")

Check if the package is installed first.

"raster" %in% installed.packages()[, "Package"]

## [1] TRUE

2 Rasters

In this exercise, we are going to use raster data that represent ‘bioclimatic variables’ from the WorldClim database (http://www.worldclim.org, Hijmans et al., 2004)

2.1 Downloading raster climate data

Worldclim database http://worldclim.org/bioclim

BIO1 = Annual Mean Temperature
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO12 = Annual Precipitation
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter

# using the raster package
bio_ly <- getData("worldclim",var="bio",res=10, path="./data/")

## Group of raster layers
class(bio_ly)

## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

names(bio_ly)

##  [1] "bio1"  "bio2"  "bio3"  "bio4"  "bio5"  "bio6"  "bio7"  "bio8" 
##  [9] "bio9"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
## [17] "bio17" "bio18" "bio19"

tempRast <- bio_ly$bio1/10
precipRast <- bio_ly$bio12

tempRast

## class       : RasterLayer 
## dimensions  : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : bio1 
## values      : -26.9, 31.4  (min, max)

2.2 Exploring a raster

Spatial resolution = The size of each cell (usually meters)
Spatial extent = X, Y coordinates of the corners of the raster in the geographic space
Projection = Spatial reference system describing how the data are “flattened” in into a 2D space
Coordinate Reference System = X,Y coordinate space associated with the projection.

“If you have the same dataset saved in two different projections, these two files won’t line up correctly”

2.3 Plotting raster

plot(tempRast)

## using the ssplot function from the raster package
spplot(tempRast)

2.3.1 Changing breaks and colours

## Changing breaks
temphist<-hist(tempRast, breaks=4)

temphist

## $breaks
## [1] -40 -20   0  20  40
## 
## $counts
## [1]  11693 181932 203695 187201
## 
## $density
## [1] 0.001000221 0.015562486 0.017424096 0.016013197
## 
## $mids
## [1] -30 -10  10  30
## 
## $xname
## [1] "v"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

# Using the raster object with different breaks
plot(tempRast, 
     breaks = temphist$breaks,
     col = brewer.pal(4, "Spectral"))

2.4 Challenge!

Using the annual mean precipitation raster, create a plot that customises the color map with 5 breaks.

2.5 Cropping Rasters

## Croping rasters based on extent
#Define the extent of the crop by clicking on the plot
# cropbox1 <- drawExtent()

cropbox1<-c(-124.7258, -66.94989, 24.49813, 49.38436)

#crop the raster, then plot the new cropped raster
Tempcrop1 <- crop(tempRast, cropbox1)
Precicrop1 <- crop(precipRast, cropbox1)

#Plot the cropped extent
plot(Tempcrop1)

2.6 Layering rasters

# create a color ramp of grey colors
plot(Tempcrop1,
    col=grey(1:100/100))

## Including both temperature and precipitation
plot(Tempcrop1,
    col=grey(1:100/100), 
    legend=F)

plot(Precicrop1,
     breaks = c(0,1000,2000,3000,4000),
     col = rev(brewer.pal(4, "Spectral")),
     alpha=0.4,
     add=T)

# include state polygons borders
# This can take a while
# USA<-getData('GADM', country='USA', level=1, path="./data/")
# plot(USA,add=TRUE)

2.7 Challenge!

Crop both the temperature and precipitation rasters for Australia
Create a plot using both cropped climatic layers

2.8 Raster calculations

#Calculate NPP using the Miami model.
npp.mia.t <- 3000/(1+exp(1.315-0.119*tempRast))
npp.mia.p <- 3000*(1-exp(-0.000664*precipRast))

nppRast<-min(stack(npp.mia.t,npp.mia.p))

plot(nppRast, main="NPP")

2.9 Extract values from a raster

2.9.1 From x,y locations

2.9.1.1 Getting records of species occurrences

With the climatic raster previously downloaded, we can determine the main climatic conditions that a particular species is adapted for.

In this exercise, we are going to use occurrence data from the Global Biodiversity Information Facility (GBIF, http://gbif.org), which have been compiled for million of species worldwide. We can access the data directly using the function gbif() from the R package dismo using the scientific name of the species of interest.

cuteSleepy<-gbif("Bradypus","variegatus")

## 1027 records found

## 0-300-600-900-1027 records downloaded

Before we extract the climatic values from our coordinates, we need to do some cleaning

cuteSleepy<-
  cuteSleepy %>% 
  filter(!is.na(lon)&!is.na(lat))

## Get rid of duplicate occurrences

dups=duplicated(cuteSleepy[, c("lon", "lat")])
cuteSleepy <-cuteSleepy[!dups, ]

Now let’s look at how the occurrences are distributed. We are going to make the map extent 1 degree around our observations.

data(wrld_simpl)

plot(wrld_simpl, xlim=c(min(cuteSleepy$lon)-1,max(cuteSleepy$lon)+1), ylim=c(min(cuteSleepy$lat)-1,max(cuteSleepy$lat)+1), axes=TRUE, col="light yellow")

points(cuteSleepy$lon, cuteSleepy$lat, col="orange", pch=16, cex=0.75)

2.9.1.2 Extract climatic variables

To extract the climatic values where our points occurr, we need first to convert the x,y point locations to SpatialPointsDataFrames. To do this, we are going to use the coordinates() function from the R package sp

This will allow us to work with the data frame as spatial data along with other spatial data – like rasters.

##Convert ocurrence data to spatial points 
cuteSleepy_geo<-cuteSleepy
coordinates(cuteSleepy_geo)<-~lon + lat

We also need to make sure that our spatial point data and the raster have the same projection and coordinate reference system. We are going to use the function proj4string() functions from the raster package. This function retrieves the projection attributes of a spatial object. In this case, we are assigning the projection system from the raster tempRast to the spatial point data.

proj4string(cuteSleepy_geo)<-proj4string(tempRast)

Now that we have a spatial object, we can display our spatial points directly to the map. We can also crop the map directly using the extent of our spatial pooints as a reference.

## Plotting climatic variable and coordinates
new_map<-crop(tempRast,extent(cuteSleepy_geo))

plot(new_map)
points(cuteSleepy_geo,pch=20,col=alpha("blue",0.5))

## Extracting data
cuteSleepy_geo$temp<-raster::extract(tempRast,cuteSleepy_geo)
hist(cuteSleepy_geo$temp)

2.10 Challenge!

Extract the values of precipitation and productivity using the cuteSleepy_geo spatial point object.
Plot the distribution of these values using histograms

Introduction to Spatial Analysis

Susy Echeverría-Londoño and Drew Kerkhoff