- Exercise 9.1 - 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 use "raster" function
library(adehabitatHR)
library(maptools)#readAsciiGrid
library(zoo) - Now open the script "Elaeo_dataprep.R" and run code directly from the script
snowy <-read.csv("SnowySamples.csv", header=T)
str(snowy)
#Clean up by deleting extraneous columns if needed
snowy <- snowy[c(-20:-38, -40:-64)]
snowy$Status <- snowy$E_schneide
#Make a spatial data frame of locations after removing outliers
coords<-data.frame(x = snowy$X_Coordina, y = snowy$Y_Coordina)
utm.crs<-"+proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
utm.spdf <- SpatialPointsDataFrame(coords= coords, data = snowy, proj4string = CRS(utm.crs)) - We now need to load some raster layers of covariates that may be related to disease occurrence
#Load DEM raster layer
dem <-raster("snowydem")
image(dem)
class(dem)
proj4string(dem)
#Now transform projections all to match DEM (i.e., Albers)
Albers.crs <-CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
snowy.spdf <-spTransform(utm.spdf, CRS=Albers.crs) - Now we can create a sampling grid that overlaps our disease locations by getting boundary box information from our locations. We added 3 rows of cells (3610 x 3 = 10830) around our outer most samples to encompass all disease samples and neighboring cells until we can figure out how to expand grid polygons in a simpler way. Alternatively, simply use the coordinates from the boundary box (bbox code) of your locations to create your sampling grid.
sublette.df <- as.data.frame(sublette.spdf)
str(sublette.df)
minx <- (min(sublette.df$x)-10830)
maxx <- (max(sublette.df$x)+10830)
miny <- (min(sublette.df$y)-10830)
maxy <- (max(sublette.df$y)+10830)
## create vectors of the x and y points
x <- seq(from = minx, to = maxx, by = 3610)
y <- seq(from = miny, to = maxy, by = 3610)
#Alternate bbox code for spatial points
# min max
#x -854784.4 -724665.0
#y 156859.0 247343.2
#Create vectors of the x and y points
#x <- seq(from = -854784.4, to = -724665.0, by = 3610)
#y <- seq(from = 156859.0, to = 247343.2, by = 3610)
#Create a grid of all pairs of coordinates (as a data.frame)
xy <- expand.grid(x = x, y = y)
class(xy)
str(xy)
#Identifiy projection before creating Spatial Points Data Frame
Albers.crs2 <-"+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
#NOTE: Albers.crs2 is needed because SPDF needs different projection command than spTransform above
grid.pts<-SpatialPointsDataFrame(coords= xy, data=xy, proj4string = CRS(Albers.crs2))
proj4string(grid.pts)
plot(grid.pts)
gridded(grid.pts)
class(grid.pts)
#Need to define points as a grid to convert to Spatial Polygons below
gridded(grid.pts) <- TRUE
gridded(grid.pts)
str(grid.pts)
plot(grid.pts)
#Convert grid points to Spatial Polygons in essence converting to a shapefile
gridsp <- as(grid.pts, "SpatialPolygons")
str(gridsp)
plot(gridsp)
class(gridsp)
summary(gridsp) - Now convert gridpts to Spatial Polygons Data Frame for added flexibility in manipulating layer
grid <- SpatialPolygonsDataFrame(gridsp, data=data.frame(id=row.names(gridsp), row.names=row.names(gridsp)))
class(grid)
plot(grid)
names.grd<-sapply(grid@polygons, function(x) slot(x,"ID"))
text(coordinates(grid), labels=sapply(slot(grid, "polygons"), function(i) slot(i, "ID")), cex=0.3)
#Let's check to see if all grid cells are the same size?
summary(grid)
getSlots(class(grid))
class(slot(grid, "polygons")[[1]])
getSlots(class(slot(grid, "polygons")[[1]]))
#Check area of each cell in the sampling grid in square meters
sapply(slot(grid, "polygons"), function(x) slot(x,"area"))
#[1] 13032100
#Grid cell size converted strto square kilometers
13032100/1000000
#[1] 13.0321 is grid cell size in square kilometers - After loading moose sample locations and creating our sampling grid it is time to work with more covariate data along with the DEM raster previously imported. We imported DEM first earlier because we need to determine the projection information early on to prepare our grid. It is easier to project moose data to fit a raster projection that vice versa.
#The layer below is a mule deer HSI raster layer without disturbance 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 - 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) - We also want to look at Land Cover data for this region and reclassify it into fewer categories for each of 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 (i.e., ll 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 - Minimize the size of the data for demonstration 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 raster 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
#############################################################
#############################################################
