1. Exercise 9.1 - Download and extract zip folder into your preferred location
  2. Set working directory to the extracted folder in R under File - Change dir...
  3. 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)

  4. 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))

  5. 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)

  6. 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)

  7. 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

  8. 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

  9. 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)

  10. 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

  11. 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
    #############################################################
    #############################################################