We may often be interested in assessing various covariates that may influence resource selection of our study animals. If we have a priori knowledge that elevation or slope may influence selection for or use of portions of the landscape then we need to create these layers for analysis. While this may not seem like a very complicated process because it is routinely done in ArcMap, those same available layers can be used and manipulated in R as in Chapter 1. We can then create slope, aspect, hillshade or other variables within R using concepts in Chapter 6 and extract those covariates for use in modeling all within the R environment.
8.3.1 Manipulating raster layers for inclusion in modeling procedures
- Exercise 8.3 - Download and extract zip folder into your preferred location
- Set working directory to the extracted folder in R under File - Change dir...
- First we need to load the packages needed for the exercise
library(sp)
library(lattice)
library(rgdal)#readOGR
library(rgeos)#gIntersection
library(raster)#to import rasters
library(adehabitatHR)
library(maptools)#readAsciiGrid - Now open the script "MD_DataPrep.R" and run code directly from the script
muleys <-read.csv("muleysexample.csv", header=T)
str(muleys)
#Remove outlier locations
newmuleys <-subset(muleys, muleys$Long > -110.90 & muleys$Lat > 37.80)
muleys <- newmuleys
newmuleys <-subset(muleys, muleys$Long < -107)
muleys <- newmuleys
#Make a spatial data frame of locations after removing outliers
coords<-data.frame(x = muleys$X, y = muleys$Y)
utm.crs<-"+proj=utm +zone=12 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
utm.spdf <- SpatialPointsDataFrame(coords= coords, data = muleys, proj4string = CRS(utm.crs))
Albers.crs <-CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
deer.spdf <-spTransform(utm.spdf, CRS=Albers.crs) - We can create a bounding box around locations and clip as above or using coordinates of box.
bbox(deer.spdf)
bb1 <- cbind(x=c(-1127964,-1127964,-1115562,-1115562,-1127964),
y=c(1718097,1724868,1724868,1718097,1718097))
AlbersSP <- SpatialPolygons(list(Polygons(list(Polygon(bb1)),"1")), proj4string=CRS(proj4string(deer.spdf)))
plot(AlbersSP)
points(deer.spdf, col="red")
#Let's buffer around the bounding box to be sure it encompasses all #locations
buffSP <- gBuffer(AlbersSP,width=1000)
plot(buffSP)
points(deer.spdf,col="red")
#NOTE: Be sure that the raster layers extend beyond AlbersSP to avoid #errors later. - Now it is time to import some raster layers of the covariates we are interested in for our analysis. Start with raster of vegetation from the 2012 NRCS Crop data that is a nice dataset that is crop specific for each year. Crop data can be found at the NRCS webpage Cropland Data Layer that can be accessed for each county of each state.
crops <-raster("crop12clip.tif")
plot(crops)
class(crops)
proj4string(crops)
crops
# Reclassify crops raster from above into 9 groups as in previous exercises.
# all values between 0 and 20 equal 1, etc.
m <- c(-Inf,0,NA,2, 7, 2, 20, 60, 3, 60, 70, 4, 110, 132, 5, 133, 150, 6, 151, 172, 7, 180, 183, 8, 189, 191, 9,192,Inf,NA)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rc <- reclassify(crops, rclmat)
plot(rc)
rc
as.matrix(table(values(rc)))#Confirms new number of vegetation categories
#Clip using buffSP polygon created earlier to reduce size of raster (if needed).
bbclip <- crop(rc, buffSP)
plot(bbclip)
points(deer.spdf,col="red")
plot(buffSP, add=T)
#Be sure all are in Albers projection before moving forward.
proj4string(bbclip)
proj4string(deer.spdf)
proj4string(buffSP) - We also want to look at elevation and covariates related to elevation (e.g., slope, aspect) from Section 6.2. These can be created directly in R using the terrain function in the raster package.
elevation <- raster("dem_albers.txt")
image(elevation, col=terrain.colors(10))
contour(elevation, add=TRUE)
#Create slope and aspect
slope = terrain(elevation,opt='slope', unit='degrees')
aspect = terrain(elevation,opt='aspect', unit='degrees')
elevation#NOTE the number of cells for all 3 layers
slope
aspect
#Clip the 3 layers within the buffSP around locations if needed
demclip <- crop(elevation, buffSP)
sloclip <- crop(slope, buffSP)
aspclip <- crop(aspect, buffSP) - Cast over all 4 layers to a Spatial Grid Data Frame to permit combining into one layer. One very important note here is that the 3 lines of code below may not work on a standard desktop PC because the memory will be maxed out. Fortunately, we have a Super Computer as we call it that was custom built by Jeff, our IT guy, with specs as follows:
############################################
Intel Core i7-3930K Sandy Bridge-E 3.2GHz LGA 2011 Six-Core Processor
G.SKILL Ripjaws Z Series 32GB (4 x 8GB) 240-Pin DDR3 SDRAM
Western Digital WD Black WD1002FAEX 1TB
7200 RPM 64MB Cache SATA 6.0Gb/s 3.5" Internal Hard Drive
######################################################################
nlcd <- as.data.frame(as(bbclip, "SpatialGridDataFrame"))
elev <- as.data.frame(as(demclip, "SpatialGridDataFrame"))
slo <- as.data.frame(as(sloclip, "SpatialGridDataFrame"))
asp <- as.data.frame(as(aspclip, "SpatialGridDataFrame"))
#Now check to be sure the number of cells in each layer are the same before proceeding the next step of combining layers.
str(elev)
str(slo)
str(asp)
#Combine elevation, slope and aspect into one layer.
layers = cbind(nlcd,elev, asp, slo)
head(layers2)
layers = layers[,-c(2,3, 5,6,8,9)]
names(layers) = c("nlcd","elevation""aspect", "slope",, "x", "y")
head(layers)
#Turn aspect into categorical variable is recommended
aspect_categorical = rep(NA, nrow(layers))
aspect_categorical[layers$aspect < 45 | layers$aspect >= 315] = "N"
aspect_categorical[layers$aspect >= 45 & layers$aspect < 135] = "E"
aspect_categorical[layers$aspect >= 135 & layers$aspect < 225] = "S"
aspect_categorical[layers$aspect >= 225 & layers$aspect < 315] = "W"
table(aspect_categorical)
table(is.na(aspect_categorical))
layers$aspect_categorical = aspect_categorical
head(layers)
write.table(layers,"layer1.txt",sep=",",col.names=TRUE, quote=FALSE)###############################################
#NOTE: Script may contains Demonstration code that will subset number of #locations to speed up processing of data during a course exercise. To #prevent this, skip this line of code: muleys <- #muleys[sample(nrow(muleys), 100),]
###############################################
- We can now begin the task of sampling each of our locations using the code below. This code was created by Ryan Nielsen of West Inc. and was very helpful in this exercise. Alternatively, we could have extracted each covariate layer by layer and included it in our dataset.
# Grab values for points created above
grab.values = function(layer, x, y){
# layer is data.frame of spatial layer, with values 'x', 'y', and ____?
# x is a vector
# y is a vector
if(length(x) != length(y)) stop("x and y lengths differ")
z = NULL
for(i in 1:length(x)){
dist = sqrt((layer$x - x[i])^2 + (layer$y-y[i])^2)
#Could adjust this line or add another line to calculate moving window or
#distance to nearest feature
z = rbind(z, layer[dist == min(dist),][1,])
}
return(z)
}
#Grab all values for used and available points based on combined layer data set that can take several minutes.
used = grab.values(layers, muleys$X, muleys$Y)
used$x = muleys$X
used$y = muleys$Y
used$animal_id = muleys$id
used$use = 1
used[1:10,] - We also need to get some measure of what is available for our mule deer population (2nd order selection) or for each mule deer (3rd order selection). We really do not understand the need for 2nd order selection unless you are looking at deer across different landscapes but hardly seems necessary for deer occupying similar areas such as our mule deer in southwestern Colorado. Below we will focus on 3rd order selection with used locations for each deer being compared to available locations randomly determined within each deer's MCP.
#Create MCP for all locations for each deer by ID (3nd order selection).
cp = mcp(deer.spdf[,2],percent=100)
as.data.frame(cp)
## Plot the home ranges and relocations.
plot(cp)
plot(deer.spdf, col=as.data.frame(deer.spdf)[,2], add=TRUE)
#Determine the habitat available using all code below
#First create random sample of points in each polygon
random <- sapply(slot(cp, 'polygons'), function(i) spsample(i, n=500, type='random', offset=c(0,0)))
plot(cp) ; points(random[[1]], col='red', pch=3, cex=.5)
#The number in the line of code above in double brackets changes polygons
#Stack into a single SpatialPoints object
random.merged <- do.call('rbind', random)
#Extract the original IDs
ids <- sapply(slot(cp, 'polygons'), function(i) slot(i, 'ID'))
#Determine the number of ACTUAL sample points generated for each polygon
newpts <- sapply(random, function(i) nrow(i@coords))
#Nice check of how many points generated per polygon
newpts
# generate a reconstituted vector of point IDs
pt_id <- rep(ids, newpts)
#Promote to SpatialPointsDataFrame
random.final <- SpatialPointsDataFrame(random.merged,
data=data.frame(poly_id=pt_id))
#Plot random final on MCPs
plot(cp) ; points(random.final, col=random.final$poly_id, pch=3, cex=0.5)
random.final
#Make random.final a data frame to extract raster covariates for each
random.df = as.data.frame(random.final,coords=coords)
names(random.df) = c("ID", "x", "y")
#Grab covariates as we did for mule deer locations above
available = grab.values(layers, random.df$x, random.df$y)
available$x = random.df$x
available$y = random.df$y
available$animal_id = pt_id
available$use = 0
available[1:10,] - Bind together mule deer locations with covariates extracted (used) and random locations within each polygon by deer ID (available) into a master dataset for modeling (data). The (use) column identifies 1 as (used) and 0 as (available)
data = rbind(available, used)
str(muleys)
#A quick check of the data to determine if correct number of records.
#100 locations used +
#100 locations available (2 animals X 50 random locations)
#= 100 #Confirmed in code below
str(data)
#''data.frame': 200 obs. of 9 variables:
# nlcd : num 7 6 7 7 7 7 7 7 7 6 ...
# elevation : int 2058 2058 2068 2068 2070 2072 2076 2062 2071 2071 ...
# aspect : num 105 278 105 80 135 ...
# slope : num 2.72 3.37 4.68 4.11 6.05 ...
# x : num -1127639 -1127610 -1127864 -1127805 -1127862 ...
# y : num 1724257 1724271 1724091 1724218 1724174 ...
# aspect_categorical: chr "E" "W" "E" "E" ...
# animal_id : chr "D12" "D12" "D12" "D12" ...
# use : num 0 0 0 0 0 0 0 0 0 0 ... - The above code is for 3rd order selection within home range of each deer. We could also look at 3rd order selection within a buffered circle around each mule deer location that is common in Discrete Choice Models. The code is similar except the initial steps of creating buffered polygons and obviously includes a lot more polygons than simply MCPs for each deer. Determining the daily distance moved was done in Chapter 3 but new code is available to estimate for each deer or all deer combined.
- First create buffered circles with radius of 500 m
#Be sure to include " _")
#Determine the number of ACTUAL sample points generated for each #polygon
buffpts <- sapply(ranbuff, function(i) nrow(i@coords))
buffpts[1:20] #Nice check of how many points generated per polygon
#Generate a reconstituted vector of point IDs
buffpt_id <- rep(buff_ids, buffpts)
#Promote to SpatialPointsDataFrame
buff.final <- SpatialPointsDataFrame(ranbuff.merged, data=data.frame(poly_id=buffpt_id))
#Plot buff.final on buffered circles
plot(settbuff); points(buff.final, col=random.final$poly_id, pch=3, cex=0.5)
buff.final
#Make 'buff.final' a data.frame
buffer.df = as.data.frame(buff.final,coords=coords)
names(buffer.df) = c("ID", "x", "y")
head(buffer.df)
# The first line below can take 5+ minutes
str(random.df)
str(buffer.df)buff_avail = grab.values(layers2, buffer.df$x, buffer.df$y)
buff_avail$x = buffer.df$x
buff_avail$y = buffer.df$y
buff_avail$animal_id = buffpt_id
buff_avail$use = 0
buff_avail[1:10,]data2 = rbind(buff_avail, used)
str(data2)#Save workspace so all analysis are available
save.image("RSF_dataprep.RData")
We are going to focus the remainder of this chapter on Selection Ratios and Resource Selection Functions (RSFs) because Selection Ratios identify a general use of habitat given what is available that can be further explored and studied through use of RSFs. Resource Selection Functions are spatially-explicit models that predict the (relative) probability of use by an animal at a given area/location during a given time, based on the environmental conditions that influence or account for selection. There are numerous types of RSFs that can be performed based on the availability of data collected during the study and there are volumes of literature devoted to the topic of resource selection and sampling designs for radiotelemetry studies (Manly et al. 2002, Cooper and Millspaugh 2001, Erickson et al. 2001, Leban et al. 2001).
