# SPATIAL DATA

# Loading packages
pcks <- list("sp", "sf", "raster", "rgdal", "rgeos", "ggmap")
sapply(pcks, require, char = TRUE)

# Loading and visualizing fire data (from a shapefile)
fires.sp <- readOGR("./_data/fmch_fire_tws.shp")
plot(fires.sp, main = "Fires around the Steese National Conservation Area", axes=TRUE, col=alpha(("red"),.5))

# Looking at different "slots" or attributes of the fire data as an `sp` object (SpatialPolygonsDataFrame)
head(fires.sp@data)
fires.sp@proj4string
fires.sp@bbox

# Now loading the fires shapefile with the new, fancy, 'sf' package
fires <- st_read("./_data/fmch_fire_tws.shp")
class(fires)

# making a quick map of the fires
mymap <- get_map(location = unname(st_bbox(fires)), maptype="hybrid")
ggmap(mymap) + 
  geom_sf(data = fires, inherit.aes = FALSE, fill = "red", alpha = 0.5)

# loading the roads data (polyline) in the 'sf' package
roads <- st_read("./_data/fmch_roads_tws.shp")
str(roads)

# a map with fires, roads and caribou
load("./_data/caribou.df.rda")
ggmap(mymap) + 
  geom_sf(data = fires, inherit.aes = FALSE, fill = "red", alpha = 0.5) +
  geom_point(data = caribou.df, aes(x = lon, y = lat, col = id)) +
  geom_sf(data = roads, inherit.aes = FALSE, col = "white")

# bringing in the elevation data (a raster) and plotting it
elevation <- raster("./_data/fmch_elev_tws.tif")
plot(elevation)

# loading the projected caribou data (which is a 'MoveStack' object from the 'move' package), 
# checking the projection of all the layers
load("./_data/caribou.utm.rda")
projection(caribou.utm)
projection(elevation)
st_crs(fires)
st_crs(roads)

# projecting the fires and roads layers to match the caribou layer
fires.utm <- st_transform(fires, projection(caribou.utm))
roads.utm <- st_transform(roads, projection(caribou.utm))

# doing the same for the raster layer
elevation.utm <- projectRaster(elevation, crs = "+proj=utm +zone=6 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", res=100)


# creating a mega figure with all the habitat data, plus the caribou data
par(mgp = c(2,.5,0), cex.axis = 0.8)
plot(elevation.utm, col = grey.colors(100, start = 0))
plot(fires.utm, col = alpha(("red"), 0.5), add = TRUE)
plot(roads.utm, col = "black", add = TRUE, lwd=2)
points(caribou.utm, pch = 19, cex = 0.5, add= TRUE, col=c("red","yellow","darkgreen","blue"))

# creating a time since fire variable
fires.utm$startDate <- lubridate::ymd(fires.utm$DiscDate)
fires.utm$endDate <- lubridate::ymd(fires.utm$OutDate)
sum(is.na(fires.utm$endDate))
fires.utm$fireDate <- pmax(fires.utm$startDate, fires.utm$endDate, na.rm = TRUE)
fires.utm$tsf <- as.numeric(difftime(lubridate::ymd("2017-08-10"), fires.utm$fireDate)/365.25)

# converting the 'sf' objects back to sp objects because 'sf' package doesn't work yet with rasters
roads.utm <- st_zm(roads.utm)
fires.utm <- st_zm(fires.utm)
roads.sp <- as(roads.utm, "Spatial")
fires.sp <- as(fires.utm, "Spatial")

# creating raster from the fire data (rasterizing), 
# by first creating an empty raster and then filling it
empty.fire.raster <- raster()
extent(empty.fire.raster) <- extent(elevation.utm) 	
projection(empty.fire.raster) <- projection(elevation.utm)
res(empty.fire.raster) <- 100
fire.raster <- rasterize(fires.sp, empty.fire.raster, field="tsf", fun=min)

# converting all "unburned" pixels to a time since fire of 100 years
fire.raster[is.na(fire.raster[])] <- 100 
plot(fire.raster)

# Don't bother doing this step in the workshop because it will take forever
# These steps create an empty raster, then fill it with values for distance to roads

# empty.roads.raster <- raster()
# extent(empty.roads.raster) <- extent(elevation.utm) 	
# projection(empty.roads.raster) <- projection(elevation.utm)
# res(empty.roads.raster) <- 100
# empty.roads.raster[] <- NA
# roads.raster <- rasterize(roads.sp, empty.roads.raster, 1)
# dist.roads <- distance(roads.raster)


# Instead, just load this file we made for you
load("./_data/dist.roads.rda")
plot(dist.roads)

# making a raster stack of elevation, time since fire, and distance to roads, and plotting it
fmch.habitat <- stack(elevation.utm, fire.raster, dist.roads)
names(fmch.habitat) <- c("elevation", "time.since.fire", "dist.roads")
plot(fmch.habitat, main=c("Elevation (m)", "Time since fire (years)", "Distance to road (m)"))