#install.packages("devtools")
#devtools::install_github("ninoxconsulting/skmoose")
library(devtools)
library(usethis)
library(skmoose)
library(sf)
library(dplyr)
library(terra)
library(foreach)Introduction
This package assists in determining moose quality and extent of usable habitat for each block given study area. This vignette provides a detailed version of the functions and how to step through each component in the process.
Preparation
Firstly we will read in the survey block layout. Ideally this geospatial file will be in .gpkg format, however a helper function is provide to convert a .shp file to .gpkg. An example data set is provided within the package. The workflow will use this example dataset, however notes are provided to assist in using your own data set.
# read in the test data
data_location <- fs::path_package("extdata", package = "skmoose")
data_file <- skmoose::skmoose_example()
bks_name <- file.path(data_location, data_file)
# Alternatively you can read in an entire data set by pointing to the file location as shown below:
#bks_name <- file.path("/home/user/Documents/r_repo/2023_moose_block_surveyR/02.Data/Tweedsmuir-SRBblocks.gpkg")
bks_name <- file.path("C:\\Users\\genev\\OneDrive\\Documents\\02.Contracts\\2023_oddjobs\\2023_moose_block_surveyR\\02.Data\\Tweedsmuir-SRBblocks.gpkg")Read in the file and generate an individual block id number
# read in and format the data
aoi <- sf::st_read(file.path(bks_name))
aoi <- aoi %>%
dplyr::mutate(bkname = seq(1,length(aoi$Name), 1)) %>%
dplyr::select(bkname)%>%
st_zm()
# generate a list with all unique block numbers to be used to iterate through the spatial extraction
blockno <- as.list(unique(aoi$bkname))Create an output folder where you want the raw data geo spatial data to be saved. Note a folder will be generated for each block number (labelled numerically).
# create an output folder:
out_dir <- "C:\\Users\\genev\\OneDrive\\Documents\\02.Contracts\\2023_oddjobs\\2023_moose_block_surveyR\\02.Data\\data"
#out_dir <- "/home//user/Documents/r_repo/2023_moose_block_surveyR/02.Data/data"
# check if the folder already exists and if it doesn't create a folder
if (!dir.exists(out_dir )){
dir.create(out_dir )
}else{
print("dir exists")
}
# write out a copy of the full aoi
sf::st_write(aoi, file.path(out_dir, "aoi.gpkg"), append = FALSE)We can now begin to extract the base data for each block. This uses the bcdata and bcmaps package to pull data directly from the bcdata catalogue. The process is set up to iterate over each block. While this is not very fast, it provides a more stable method to reduce the incidence of timing out during large data downloads, specifically with the VRI dataset.
WARNING: this is the most time consuming part of the process and will take some time to run if extracting many blocks. Best to run overnight. If there is a time out error or
# set up parrallel processing parameters
registerDoSEQ()
# loop through each block and extract the spatial basedata. These will be filtered based on criteria provided
# and a copy saved under the respective file folder.
# 247 error
basedata <- foreach(x = blockno) %dopar% {
# test line - if you want to run a single block
#x = blockno[92] #42, #50, #38
print(x)
tmp_aoi <- aoi %>% dplyr::filter(bkname == x)
#st_bbox(tmp_aoi)
temp_out_dir <- file.path(out_dir, x)
if (!dir.exists(temp_out_dir)){
dir.create(temp_out_dir )
}else{
print("dir exists")
}
get_basedata(in_aoi = tmp_aoi, out_path = temp_out_dir, overwrite = FALSE)
#print(x)
}Extract uninhabitable area
Once all the base data is extracted we can filter uninhabitable areas for moose. This includes:
- Rock and Ice
- Waterbodies >1 km2
- Elevations >1300 m
- Slopes >55 degrees
If there are sections within the study area that meet these criteria they will be saved as a new .gpkg within the folder of the particular block. If no additional files are written out, then the criteria are not met.
# set thresholds
eval_threshold = 1300
slope_threshold = 55
uninhab <- foreach(x = blockno) %dopar% {
#x = blockno[26]
tmp_aoi <- aoi %>% dplyr::filter(bkname == x)
temp_out_dir <- file.path(out_dir, x)
# 1. extract rock and ice areas
rockice <- st_read(file.path(temp_out_dir, "vri.gpkg")) %>%
dplyr::filter(BCLCS_LEVEL_3 == "A") %>%
dplyr::select(id) %>%
st_union()%>%
st_cast("MULTIPOLYGON")
if(length(st_is_empty(rockice)< 1)){
sf::st_write(rockice, file.path(temp_out_dir, "rockice.gpkg"), append = FALSE)}
# 2. filter lakes with >1 km2
if(file.exists(file.path(temp_out_dir,"lakes.gpkg"))){
largelakes <- st_read(file.path(temp_out_dir, "lakes.gpkg")) %>%
dplyr::filter(AREA_HA > 100) %>%
st_union()
if(length(st_is_empty(largelakes)< 1)){
st_write(largelakes, file.path(temp_out_dir, "largelakes.gpkg"), append = FALSE)}
}
# 3. filter elevations above threshold value
trim <- terra::rast(file.path(temp_out_dir, "dem.tif"))
high_elev_sf <- high_elev(trim, eval_threshold )
if(length(st_is_empty(high_elev_sf )< 1)){
st_write(high_elev_sf , file.path(temp_out_dir, "high_elevation.gpkg"), append = FALSE)
}
# 4. filter high slope
steep_sf <- steep_slope(trim, slope_threshold)
if(length(st_is_empty(steep_sf )< 1)){
st_write(steep_sf, file.path(temp_out_dir, "steep.gpkg"), append = FALSE)
}
}Combine uninhabitable area.
For each block we then combine the uninhabitable area into a single polygon and clip to the block boundary. The output file will be called “uninhabitable.gpkg”
Extract Moose Habitat
We now compiled all the moose habitat areas. This involves reading in the base data and filtering for
Deciduous tree species, queried in VRI: [SPECIES_CD LIKE ‘A%’ OR SPECIES_CD = ‘EP’ OR SPECIES_CD = ‘SB’ OR SPECIES__1 LIKE ‘A%’ OR SPECIES__1 = ‘EP’ OR SPECIES__1 = ‘SB’ OR SPECIES__3 LIKE ‘A%’ OR SPECIES__3 = ‘EP’ OR SPECIES__3 =‘SB’ OR SPECIES__5 LIKE ‘A%’ OR SPECIES__5 LIKE ‘E%’]
Early seral/shrub dynamic habitats: Fires that are >10 and <25 years old
Cutblocks that are >5 and < 40 years old
Buffered 3-8 order streams from FWA stream dataset by 150m (dissolved)
9th order streams buffered by 500 meters.
Wetlands from FWA ≤ 1 km2
Skeena Wildlife Ecological Resource Model- Winter forage output
# set thresholds
cutblock_min_yr = 5
cutblock_max_yr = 25
fire_min_yr = 10
fire_max_yr = 25
habitable <- foreach(x = blockno) %dopar% {
#x = blockno[20]
print(x)
tmp_aoi <- aoi %>% dplyr::filter(bkname == x)
temp_out_dir <- file.path(out_dir, x)
# get deciduous leading species
vri <- st_read(file.path(temp_out_dir, "vri.gpkg"),quiet = TRUE)
species_codes = c("AT", "AC","EP","SB")
decid <- vri_browse(vri, species_codes)
if(length(st_is_empty(decid)< 1)){
print("contains deciduous")
st_write(decid , file.path(temp_out_dir, "vri_decid.gpkg") , append = FALSE,quiet = TRUE)
}
# get recent harvest
if("cutblocks.gpkg" %in% list.files(temp_out_dir)){
cutblocks <- st_read(file.path(temp_out_dir, "cutblocks.gpkg"),quiet = TRUE)
cutblocks_yrs <- cutblocks_recent(cutblocks, cutblock_min_yr, cutblock_max_yr)
if(length(st_is_empty(cutblocks_yrs)) > 0 ){
print("contains cutblocks")
st_write(cutblocks_yrs, file.path(temp_out_dir, "cutblocks_filtered.gpkg"), quiet = TRUE, append = FALSE)
}
}
# get recent fire years
if("fire.gpkg" %in% list.files(temp_out_dir)){
fires <- st_read(file.path(temp_out_dir, "fire.gpkg"),quiet = TRUE)
fires_filtered <- fires_recent(fires, fire_min_yr, fire_max_yr)
if(length(st_is_empty(fires_filtered)) > 0 ){
print("contains fire")
st_write( fires_filtered, file.path(temp_out_dir, "fires_filtered.gpkg"), append = FALSE,quiet = TRUE)
}
}
# get stream order and buffer by stream order, stream order 3 - 8 is 150m buffer, stream order 9 bufferd by 500m
if("streams.gpkg" %in% list.files(temp_out_dir)){
streams <- st_read(file.path(temp_out_dir, "streams.gpkg"), quiet = TRUE)
stream38 <- buffer_streams(streams)
stream9 <- buffer_streams(streams, 9, 500)
if(length(st_is_empty(stream38 )) > 0 ){
print("contains streams")
st_write(stream38 , file.path(temp_out_dir, "streams3_8.gpkg"),quiet = TRUE, append = FALSE)
}
if(length(st_is_empty(stream9 )) > 0 ){
st_write(stream9 , file.path(temp_out_dir, "streams9.gpkg"), quiet = TRUE, append = FALSE)
}
}
# select small lakes and all wetlands (previously filtered to < 100 Ha Area.
if("lakes.gpkg" %in% list.files(temp_out_dir)){
smalllakes <- st_read(file.path(temp_out_dir, "lakes.gpkg"),quiet = TRUE)%>%
dplyr::filter(AREA_HA < 100) %>%
st_union()
if(length(st_is_empty(smalllakes)) > 0 ){
print("contains lakes")
st_write(smalllakes, file.path(temp_out_dir, "smalllakes.gpkg"), quiet = TRUE, append = FALSE)
}
}
}Merge habitable layers into single file per block
For each block we then combine all the outputs into a single moose habitat layer called “habitat.gpkg”.
# merge uninhabitable area
habitable <- foreach(x = blockno) %dopar% {
#x = blockno[20]
tmp_aoi <- aoi %>% dplyr::filter(bkname == x)
temp_out_dir <- file.path(out_dir, x)
merge_habit(tmp_aoi, temp_out_dir)
# remove an uninhabitable areas from within habitable
if(file.exists(file.path(temp_out_dir,"uninhabitable.gpkg"))){
uninh <- st_read(file.path(temp_out_dir, "uninhabitable.gpkg"),quiet = TRUE)
hab <- st_read(file.path(temp_out_dir, "habitable.gpkg"),quiet = TRUE)
habit <- st_intersection(uninh, hab)
final_hab <- st_difference(hab, uninh)
st_write(final_hab, file.path(temp_out_dir, "habitable.gpkg"), quiet = TRUE, append = FALSE)
}
}
Check for Burn Severity extents
Calculate the burn severity for each of the plots and out put a gpkg called non_wet_fire.gpkg
burn_strat <- foreach(x = blockno) %dopar% {
#x = blockno[42]
tmp_aoi <- aoi %>% dplyr::filter(bkname == x)
temp_out_dir <- file.path(out_dir, x)
# merge_habit(tmp_aoi, temp_out_dir)
# remove an uninhabitable areas from within habitable
if(file.exists(file.path(temp_out_dir,"fire_int.gpkg"))){
fire_int <- st_read(file.path(temp_out_dir, "fire_int.gpkg"),quiet = TRUE)%>%
sf::st_union()
if(file.exists(file.path(temp_out_dir,"streams3_8.gpkg"))){
stream <- st_read(file.path(temp_out_dir, "streams3_8.gpkg"),quiet = TRUE)
non_wet_fire<- rmapshaper::ms_erase(fire_int, stream)
} else {
non_wet_fire <- fire_int
}
#final_hab <- st_difference(hab, uninh)
st_write(non_wet_fire, file.path(temp_out_dir, "non_wet_fire.gpkg"), quiet = TRUE, append = FALSE)
}
}Calculate areas for habitat and non-habitable
For each block we will read in the habitable and inhabitable areas and calculate the areas and proportions for each block. A table is
out_table <- calculate_areas(blockno, aoi, out_dir)
write.csv(out_table, file.path(out_dir, "moose_stratification_outputs_raw.csv"), row.names = FALSE)
out_table <- read.csv(file.path(out_dir, "moose_stratification_outputs_raw.csv"))Optional: Moose Habitat Strata Classification
Assign Moose blocks into stratification of High, Medium, Low for further review. This will be based on the Proportion of Moose Habitat/Block Area (prop_habit_block_km2)
Two proposed methods include 1) Assign values based on set thresholds, i.e: 0-40% =L, #41-70%= Moderate, 71-100%, 71-100% = High 2) Assign values based on distribution of all block values.
# Option 1: set thresholds based on the top value in each catergorty.
# Note you can adjust low, med and high values within the function
out <- assign_categories(out_table, low = 0.4, med = 0.7, high = 1, burn_strat = 0.5)
# OR
# option 2: estimate thresholds based on the quantiles ie:
outq <- assign_categories(out_table, quartile = TRUE)
write.csv(out, file.path(out_dir, "moose_stratification_outputs_classed.csv"), row.names = FALSE)
#out <- read.csv(file.path(out_dir, "moose_stratification_outputs_classed.csv"))
# write out a spatial output with blocks and values;
sp_out <- merge(aoi, out)
st_write(sp_out, file.path(out_dir, "moose_stratification_output.gpkg"),quiet = TRUE,append = FALSE)
st_write(sp_out, file.path(out_dir, "moose_stratification_output.shp"),append = FALSE)