8. Continue with same script loaded after finishing Exercise 9.1
9. We loaded moose sample locations and created our sampling grid above so now it is time to work with covariate data that was provided in various forms. We imported DEM in the last exercise because we needed to determine the projection information early on to prepare our grid. It is easier to project moose data to fit a raster projection than vice versa so now let's continue adding additional covariates from raster data.
#The layer below is a mule deer HSI raster layer without disturbance created based on data from Sawyer et al. 2009 incorporated into a layer because mule deer are considered host for the parasite we are investigating
nodis <-raster("snowynodis")
nodis
plot(nodis)
summary(nodis)
#Need to remove NoData from mule deer HSI layer
nodis[is.na(nodis[])] <- 0
10. Using functions from the raster package, we can calculate slope and aspect from DEM layer imported above
slope = terrain(dem,opt='slope', unit='degrees')
aspect = terrain(dem,opt='aspect', unit='degrees')
dem #Now let's see metadata for each layer
slope
aspect
plot(dem)
plot(grid, add=T)
11. We also want to look at Land Cover data for this region and reclassify it into fewer categories for comparison and manipulation
nlcdall <- raster("nlcd_snowy")
nlcdall #Look at raster values for the habitat layer
#Values range from 11 to 95
#Or plot to visualize categories in legend
plot(nlcdall)
#Reclassify the values into 7 groups
#all values between 0 and 20 equal 1, etc.
m <- c(0, 19, 1, 20, 39, 2, 40, 50, 3, 51,68, 4, 69,79, 5, 80, 88, 6, 89, 99, 7)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rc <- reclassify(nlcdall, rclmat)
plot(rc) #Now only 7 categories
rc #Now only 7 categories
class(rc)
#Check to be sure all raster have same extent for Stack creation
compareRaster(dem,slope,aspect,nodis,rc)
#[1] TRUE
12. Minimize the size of the data for demonstartion purposes
##################################################
##NOTE: Code in this box was simply for demonstration purposes to reduce ##overall time for processing during class. Skip this section of code if using your
##own data and your computer has the appropriate processing capabilities.
###################################################First we will clip out the raster layers by zooming into only a few locations
plot(rc)
plot(grid, add=T)
points(snowy.spdf)
#Code below is used to just zoom in on grid using raster layer
e <- drawExtent()
#click on top left of crop box and bottom right of crop box create zoom
newclip <- crop(rc,e)
plot(newclip)
plot(grid, add=T)
points(snowy.spdf, col="red")
#Clip locations within extent of raster
samp_clip <- crop(snowy.spdf,newclip)
plot(newclip)
plot(samp_clip, add=T)
grid_clip <- crop(grid, newclip)
plot(grid_clip, add=T)
slope2 <- crop(slope,newclip)
aspect2 <- crop(aspect,newclip)
dem2 <- crop(dem,newclip)
HSI <- crop(nodis, newclip)
#Check to be sure all rasters have same extent for Stack creation
compareRaster(dem2,slope2,aspect2,HSI,newclip)
#[1] TRUE
grid <- grid_clip #rename clipped grid as grid to match code below
rc <- newclip #rename clipped Land Cover as "rc" to match code below
snowy.spdf <- samp_clip
#Create a Stack of all Raster layers
#This will take a long long time if rasters have a large extent
r <- stack(list(dem=dem2, slope=slope2, aspect=aspect2, "mule deer HSI"=HSI, nlcd=newclip))
#END Demonstration code
##################################################
13. Now we want to combine all raster layers into a multi-layered raster called a "stack" below before proceeding if using your own data. Otherwise skip Line 13 if using demonstration code above and continue on with Line 14
#Create a Stack of all Rasters
#This will take a long long time if rasters have a large extent
r <- stack(list(dem=dem, slope=slope, aspect=aspect, "HSI"=nodis, nlcd=rc))
nlayers(r) #Show how many layers are in the "r" stack
plot(r) #Visualize "r"
names(r) #Names of each raster in the "r" stack
