A. EOSD raster - cropped and masked

This is (in my experience) the fussiest step. We’ll start with the fire modified EOSD:

eosd <- raster("_data/rasters/EOSD_withfire")
plot(eosd)
plot(hayriver.lines %>% st_transform(projection(eosd)), add= TRUE)  

A lot of this data is outside the range of our eventual RSF, so we crop it and mask it to the extent of the movement data.

eosd.cropped <- crop(eosd, as(hayriver.buffer, "Spatial"))
eosd.masked <- mask(eosd.cropped, hayriver.buffer)
plot(eosd.masked)
plot(hayriver.lines %>% st_transform(projection(eosd)), add= TRUE)

Recall, the details of the raster as as follows:

eosd.masked

## class       : RasterLayer 
## dimensions  : 5731, 7297, 41819107  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -423996.4, -205086.4, 2214172, 2386102  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-112 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : in memory
## names       : EOSD_NWT_HRL_clip 
## values      : 0, 250  (min, max)
## attributes  :
##        ID OBJECTID Value Description   Color   Count
##  from:  0        1    12      Shadow #333333     365
##  to  : 19       20   250 Recent Burn #FF5555 6731695

i.e.~ it is at a 30m resolution spanning:

(bbox(eosd.masked) %>% apply(1, diff))/1e3

##     s1     s2 
## 218.91 171.93

about 220 x 171 km, with about 42 million cells. In other words, it’s LARGE. Eventually - to create a predicitive map - ALL the other raster layers will have to exactly match the resolution and extent of this raster, which means they will have to be LARGE as well.

B. BEAD line density

To get the line densities, we adapt the code from Step 5b to obtain a raster of linear densities. Specifically, we load the linear elements data and compute a “density” of linear elements.

This is, if you recall, a computationally intensive step, but some of the functions we wrote earlier will be useful for this operation. These two functions are getLineDensity and makeLineDensityRaster. The code is included below (click on the button to reveal), and is only slightly modified from the versions here.

getLineDensity():

getLineDensity <- function(pt, lines, r, plotme = FALSE){
  
  pt.circle <- st_buffer(pt, r) 
  intersect <- st_intersection(pt.circle, lines)
  length <- sum(st_length(intersect))/1e3 # in km
  
  density <- (as.numeric(length)/(pi*(r/1e3)^2)) # in km/km^2
  
  if(plotme){
    plot(st_geometry(lines))
    plot(pt.circle, add = TRUE, border = 2)
    plot(intersect, add = TRUE, col =3)
  }
  
  return(density)
}

makeLineDensityRaster():

makeLineDensityRaster <- function(lines, resolution, radius){
  
  require(raster)
  require(doParallel)
  require(magrittr)
  require(sf)
  require(spex)
  
  # data prep
  density.raster <- raster(extent(lines), 
                           resolution = resolution, 
                           crs = st_crs(lines))
  
  print("Getting polygons")
  
  density.poly <- polygonize(density.raster) 
  density.poly.sf <- st_as_sf(density.poly)
  st_crs(density.poly.sf) <- st_crs(lines)
  
  # buffering lines
  
  print("Tidying outside of buffer")
   
  lines.buffered <- st_buffer(lines, radius) %>% st_union
  which.in.buffer <- st_intersection(density.poly.sf, lines.buffered) %>% 
    row.names %>% as.numeric
  
  # cluster prep
  
  density.poly.list <- split(density.poly.sf[which.in.buffer,], 1:length(which.in.buffer))

  print("Making cluster")
  
  cl <- makeCluster(detectCores()-1)
  clusterEvalQ(cl, library(sf))
  clusterExport(cl, "getLineDensity")
  
  print("Computing densities")
  densities.in.buffer <- parSapply(cl, density.poly.list, 
                                   getLineDensity, lines = lines, r = radius)
  values(density.raster)[which.in.buffer] <- densities.in.buffer
  density.raster[is.na(density.raster)] <- 0
  stopCluster(cl)   
  
  # check that projection is correct
  projection(density.raster) <- st_crs(lines)
  return(density.raster)
}

A few notes:

1. makeLineDensityRaster is made significantly faster because it uses paralleliztion … which might not always work within the function as its written (e.g. it may be a bit operating system dependent), in which case the function might need to be unpacked.
2. makeLineDensityRaster also requires the nes spex package, which has a function called polygonize() which is a much faster version of the rasterToPolygons() function I used in the earlier attempts.

Here, I run this code on the the BEAD.lines data at a 2km distance (which seems, typically, to be relevant for ungulates). Three steps here:

1. Load / project / and crop the data to the masked EOSD layer

BEAD.lines <- st_read("_data/BEAD/BEADlines2015_NT1_Bistcho_Calendar_MaxH_Yates_HRL_clip.shp")
BEAD.lines.cropped <- st_crop(BEAD.lines %>% st_transform(projection(eosd.masked)), extent(eosd.masked))

2. Compute the densities

roadDensity <- makeLineDensityRaster(lines = BEAD.lines.cropped, resolution = 300, 2e3)

This is the SLOW step! Make sure you save the resulting raster!

roadDensity <- writeRaster(roadDensity, "_data/RoadDensity")

3. Match the rasters

Note, we do this to a 300 m resolution (30 m is extremely high for this application!) But, to match the rasters, it is easy enough to increase the resolution via the resample command:

roadDensity.matched <- resample(roadDensity, eosd.masked)

This is the KEY function for building layers for an RSF since it makes the layers have the exact same properties (extent, resolution and projection).

To guarantee that this worked, we combine the rasters into an object called a rasterBrick:

## RSF.brick <- brick(eosd.masked, roadDensity.matched)

Note, there are two very similar objects “stacks” and “bricks” - craeated and accessed viw the stack() and brick() functions respectively. I am not very clear on the difference between the two, except that bricks are apparently faster - presumably at the cost of retaining some information.

In any case, the brick is a sort of list. Here are all of its properties.

RSF.brick

## class       : RasterBrick 
## dimensions  : 5731, 7297, 41819107, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : -423996.4, -205086.4, 2214172, 2386102  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-112 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : c:/eli/caribou/NWT/data/RSF_stack.grd 
## names       : EOSD_NWT_HRL_clip,      layer 
## min values  :         0.0000000, -0.0704177 
## max values  :        250.000000,   4.586764

names(RSF.brick)

## [1] "EOSD_NWT_HRL_clip" "layer"

And if you want to access a single layer, you can do it my its location in the list:

RSF.brick[[1]]

## class       : RasterLayer 
## band        : 1  (of  2  bands)
## dimensions  : 5731, 7297, 41819107  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -423996.4, -205086.4, 2214172, 2386102  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-112 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : c:/eli/caribou/NWT/data/RSF_stack.grd 
## names       : EOSD_NWT_HRL_clip 
## values      : 0, 250  (min, max)

or by its name:

RSF.brick$EOSD_NWT_HRL_clip

## class       : RasterLayer 
## band        : 1  (of  2  bands)
## dimensions  : 5731, 7297, 41819107  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -423996.4, -205086.4, 2214172, 2386102  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-112 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : c:/eli/caribou/NWT/data/RSF_stack.grd 
## names       : EOSD_NWT_HRL_clip 
## values      : 0, 250  (min, max)

The names are awkward … let’s simplify:

names(RSF.brick) <- c("EOSD", "RoadDensity")

And now you can save the brick using the writeRaster command:

writeRaster(RSF.brick, filename = "c:/eli/caribou/NWT/data/RSF_brick_v1.tif")

Note: The .tif format is relatively compact (e.g. compared to R’s native .grd format), but this object is still quite large, about 95 MB. So be warned!

As needed, you can then load the raster brick using the brick command:

RSF.brick <- brick("c:/eli/caribou/NWT/data/RSF_brick_v1.tif")

C. Digital Elevation Model

Just to make the analysis a little bit more interesting, I also include a digital elevation model. These particular data are from a 30m global elevelation model called “ALOS World 3D (AW3D30)” documentation here. I extracted these data in a relatively tight crop around the Hay River Lowlands at 100 m resolution (for tractability).

Note: I tried obtaining 30m resolution data from the new and much vaunted Arctic DEM - but there was really too much missing information, possibly because this study area is really a the edge of the Arctic DEM’s coverage.

Here are the raw DEM data:

DEM <- raster("_data/HayRiverDEM_100m.tif")
plot(DEM)
contour(DEM, add = TRUE, nlevels = 20)

Always fun (and easy) to layer these on top of a mapview:

mapview(DEM, map.types = "Esri.WorldImagery")

There isn’t too much going on here. Aside from Great Slave Lake to the northeast, there is a relatively high (600-800 m) plateau jutting in from the south.

Matching this raster to our brick is straightforward. But note: it is faster to first reproject this raster at a coarse resolution (e.g. the approximate native 100 m) and THEN resample to match the brick, rather than directly project this to. In other words, this:

DEM.matched <- projectRaster(DEM, RSF.brick)

is VERY slow. But this:

DEM.projected <- projectRaster(DEM, crs = projection(RSF.brick), res = 100)
DEM.matched <- resample(DEM.projected, RSF.brick)

is reasonable. We can add this layer to our final brick:

RSF.brick[[3]] <- DEM.matched
names(RSF.brick) <- c("EOSD", "RoadDensity", "DEM")
writeRaster(RSF.brick, 
            filename = "c:/eli/caribou/NWT/data/RSF_brick_v2.tif")

Final tidying

Just a few more items for final data prep.

We probably don’t need the geometry column (that’s a vestige of its sf origins, and the coordinates themselves are in the X and Y columns), but it might be nice to have the lat-long coordinates. The following code is a bit fussy (due to some simple feature weirdness) but ultimately returns a column of lat/long:

ll <- hayriver.RSF$geometry %>% st_sfc %>% st_transform(4326) %>% st_coordinates %>% data.frame %>% plyr::rename(c(X = "longitude", Y = "latitude"))
hayriver.RSF <- cbind(hayriver.RSF, ll) %>% mutate(geometry = NULL)

It would be very convenient to have a column of EOSD land layers by name. Recall that they are stored in the EOSD raster levels:

levels(eosd)

## [[1]]
##    ID OBJECTID Value         Description   Color    Count
## 1   0        1    12              Shadow #333333      365
## 2   1        2    20               Water #3333FF  5559064
## 3   2        3    32         Rock/Rubble #CCCCCC       67
## 4   3        4    33 Exposed/Barren Land #996633    13635
## 5   4        5    34           Developed #CC6699    64239
## 6   5        6    40             Bryoids #FFCCFF    96789
## 7   6        7    51          Shrub Tall #FFFF33    72952
## 8   7        8    52           Shrub Low #FFFF99   668629
## 9   8        9    81       Wetland-treed #9933CC  2358766
## 10  9       10    82       Wetland-shrub #CC66FF  4242600
## 11 10       11    83        Wetland-herb #CC99FF  1160373
## 12 11       12   100               Herbs #CCFF33   514846
## 13 12       13   211    Coniferous-dense #006600  6715586
## 14 13       14   212     Coniferous-open #009900 11162064
## 15 14       15   213   Coniferous-sparse #00CC66  3236051
## 16 15       16   221     Broadleaf-dense #00CC00   917321
## 17 16       17   222      Broadleaf-open #99FF99   985631
## 18 17       18   231     Mixedwood-dense #CC9900  3088338
## 19 18       19   232      Mixedwood-open #CCCC66   600763
## 20 19       20   250         Recent Burn #FF5555  6731695

We need to match our EOSD column to the Value column of this object, and extract the Description. We do that (in one line) by merging and renaming:

hayriver.RSF <- merge(hayriver.RSF, levels(eosd)[[1]][,c("Value","Description")], by.x = "EOSD", by.y = "Value") %>% plyr::rename(c(Description = "EOSD.class")) 

Here’s a useful quick survey of the landscape data:

with(hayriver.RSF, table(EOSD.class, Used))

##                      Used
## EOSD.class            FALSE TRUE
##   Broadleaf-dense         7    3
##   Broadleaf-open          6    6
##   Bryoids                 3    3
##   Coniferous-dense      161 1373
##   Coniferous-open       217 3174
##   Coniferous-sparse      40  460
##   Developed               3    0
##   Exposed/Barren Land     0    0
##   Herbs                   0    0
##   Mixedwood-dense        55  164
##   Mixedwood-open         12   72
##   Recent Burn            72  650
##   Rock/Rubble             0    0
##   Shadow                  0    0
##   Shrub Low               4   32
##   Shrub Tall              1    0
##   Water                  91   91
##   Wetland-herb           25  682
##   Wetland-shrub         105 1404
##   Wetland-treed          53  724

Also - the nick_name is still retaining all 201 levels of the original data frame. Lets drop those unused levels:

hayriver.RSF$nick_name <- droplevels(hayriver.RSF$nick_name)

Finally, often this processing shuffles the order of the data. It would be good to sort the data, by ID and by time, so that they are ordered in a reasonable way:

hayriver.RSF <- arrange(hayriver.RSF, nick_name, timestamp)

Finally, be sure to save this data frame:

save(hayriver.RSF, file = "_data/hayriver_RSF.rda")