Basic Tutorial

# First step is to load the libraries. Not all of these libraries are stricly
# needed; some are used for convenience and visualization for this tutorial.
library("samc")
library("raster")
library("viridisLite")

# "Load" the data. In this case we are using data built into the package.
# In practice, users will likely load raster data using the raster() function
# from the raster package.
res_data <- samc::ex_res_data
abs_data <- samc::ex_abs_data
occ_data <- samc::ex_occ_data


# Plot the data and make sure it looks good. The built-in data is in matrices, 
# so we use the raster() function to help with the plotting. Note that when
# matrices are used by the package, it sets the extents based on the number of
# rows/cols. We do the same thing here when converting to a raster, otherwise
# the default extents will be (0,1) for both x and y, which is not only
# uninformative, but can result in "stretching" when visualizing datasets
# based non-square matrices.
plot(raster(res_data, xmn = 1, xmx = ncol(res_data), ymn = 1, ymx = nrow(res_data)),
     main = "Example Resistance Data", xlab = "x", ylab = "y", col = viridis(256))

plot(raster(abs_data, xmn = 1, xmx = ncol(abs_data), ymn = 1, ymx = nrow(abs_data)),
     main = "Example Absorption Data", xlab = "x", ylab = "y", col = viridis(256))

plot(raster(occ_data, xmn = 1, xmx = ncol(occ_data), ymn = 1, ymx = nrow(occ_data)),
     main = "Example Occupancy Data", xlab = "x", ylab = "y", col = viridis(256))

# Setup the details for our transition function
tr <- list(fun = function(x) 1/mean(x), # Function for calculating transition probabilities
           dir = 8, # Directions of the transitions. Either 4 or 8.
           sym = TRUE) # Is the function symmetric?

# Create a `samc-class` object using the resistance and absorption data. We use the
# recipricol of the arithmetic mean for calculating the transition matrix. Note,
# the input data here are matrices, not RasterLayers.
samc_obj <- samc(res_data, abs_data, tr_args = tr)


# Print out the samc object and make sure everything is filled out. Try to
# double check some of the values, such as the nrows/ncols of the landscape
# data. The dimensions of the matrix (slot p) should be the number of non-NA
# cells in your data +1. In this case, our data has 2624 non-NA cells, so the
# matrix should be 2625 x 2625
str(samc_obj)
#> Formal class 'samc' [package "samc"] with 7 slots
#>   ..@ data    :Formal class 'samc_data' [package "samc"] with 3 slots
#>   .. .. ..@ q    :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. .. .. ..@ i       : int [1:19786] 1 115 116 117 0 1 2 116 117 118 ...
#>   .. .. .. .. ..@ p       : int [1:2625] 0 4 10 16 22 28 34 40 46 52 ...
#>   .. .. .. .. ..@ Dim     : int [1:2] 2624 2624
#>   .. .. .. .. ..@ Dimnames:List of 2
#>   .. .. .. .. .. ..$ : chr [1:2624] "1" "2" "3" "4" ...
#>   .. .. .. .. .. ..$ : chr [1:2624] "1" "2" "3" "4" ...
#>   .. .. .. .. ..@ x       : num [1:19786] 0.226 0.138 0.163 0.104 0.293 ...
#>   .. .. .. .. ..@ factors : list()
#>   .. .. ..@ t_abs: Named num [1:2624] 2e-04 2e-04 2e-04 2e-04 2e-04 ...
#>   .. .. .. ..- attr(*, "names")= chr [1:2624] "1" "2" "3" "4" ...
#>   .. .. ..@ c_abs: num[0 , 0 ] 
#>   ..@ source  : chr "map"
#>   ..@ map     :Formal class 'RasterLayer' [package "raster"] with 12 slots
#>   .. .. ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
#>   .. .. .. .. ..@ name        : chr ""
#>   .. .. .. .. ..@ datanotation: chr "FLT4S"
#>   .. .. .. .. ..@ byteorder   : chr "little"
#>   .. .. .. .. ..@ nodatavalue : num -Inf
#>   .. .. .. .. ..@ NAchanged   : logi FALSE
#>   .. .. .. .. ..@ nbands      : int 1
#>   .. .. .. .. ..@ bandorder   : chr "BIL"
#>   .. .. .. .. ..@ offset      : int 0
#>   .. .. .. .. ..@ toptobottom : logi TRUE
#>   .. .. .. .. ..@ blockrows   : int 0
#>   .. .. .. .. ..@ blockcols   : int 0
#>   .. .. .. .. ..@ driver      : chr ""
#>   .. .. .. .. ..@ open        : logi FALSE
#>   .. .. ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
#>   .. .. .. .. ..@ values    : logi [1:6868] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>   .. .. .. .. ..@ offset    : num 0
#>   .. .. .. .. ..@ gain      : num 1
#>   .. .. .. .. ..@ inmemory  : logi TRUE
#>   .. .. .. .. ..@ fromdisk  : logi FALSE
#>   .. .. .. .. ..@ isfactor  : logi FALSE
#>   .. .. .. .. ..@ attributes: list()
#>   .. .. .. .. ..@ haveminmax: logi TRUE
#>   .. .. .. .. ..@ min       : int 0
#>   .. .. .. .. ..@ max       : int 1
#>   .. .. .. .. ..@ band      : int 1
#>   .. .. .. .. ..@ unit      : chr ""
#>   .. .. .. .. ..@ names     : chr ""
#>   .. .. ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
#>   .. .. .. .. ..@ type      : chr(0) 
#>   .. .. .. .. ..@ values    : logi(0) 
#>   .. .. .. .. ..@ color     : logi(0) 
#>   .. .. .. .. ..@ names     : logi(0) 
#>   .. .. .. .. ..@ colortable: logi(0) 
#>   .. .. ..@ title   : chr(0) 
#>   .. .. ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
#>   .. .. .. .. ..@ xmin: num 0.5
#>   .. .. .. .. ..@ xmax: num 202
#>   .. .. .. .. ..@ ymin: num 0.5
#>   .. .. .. .. ..@ ymax: num 34.5
#>   .. .. ..@ rotated : logi FALSE
#>   .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
#>   .. .. .. .. ..@ geotrans: num(0) 
#>   .. .. .. .. ..@ transfun:function ()  
#>   .. .. ..@ ncols   : int 202
#>   .. .. ..@ nrows   : int 34
#>   .. .. ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
#>   .. .. .. .. ..@ projargs: chr NA
#>   .. .. ..@ history : list()
#>   .. .. ..@ z       : list()
#>   ..@ clumps  : int 1
#>   ..@ override: logi FALSE
#>   ..@ threads : num 1
#>   ..@ .cache  :<environment: 0x55bb1d792ee0>

# Convert the occupancy data to probability of occurrence
occ_prob_data <- occ_data / sum(occ_data, na.rm = TRUE)


# Calculate short- and long-term mortality metrics and long-term dispersal
short_mort <- mortality(samc_obj, occ_prob_data, time = 4800)
long_mort <- mortality(samc_obj, occ_prob_data)
long_disp <- dispersal(samc_obj, occ_prob_data)
#> 
#> Cached diagonal not found.
#> Performing setup. This can take several minutes... Complete.
#> Calculating matrix inverse diagonal...
#> 
Computing: 100% (done)                         
#> 
Complete                                                      
#> Diagonal has been cached. Continuing with metric calculation...

# Create rasters using the vector result data for plotting.
short_mort_map <- map(samc_obj, short_mort)
long_mort_map <- map(samc_obj, long_mort)
long_disp_map <- map(samc_obj, long_disp)


# Plot the mortality and dispersal results
plot(short_mort_map, main = "Short-term Mortality", xlab = "x", ylab = "y", col = viridis(256))
plot(long_mort_map, main = "Long-term Mortality", xlab = "x", ylab = "y", col = viridis(256))
plot(long_disp_map, main = "Long-term Dispersal", xlab = "x", ylab = "y", col = viridis(256))