IDFG Regional Projects

Headquarters

While most of my projects for IDFG generally have a statewide-focus, I occasionally help with specific asks from bureau staff.

FSVM predictions for Lodgepole pine distributions

Using the fine scale vegetation model I developed, I was able to provide our staff ecologist with predictions for lodgepole pine to help guide their field sampling efforts.

Lodepole pine predicted probability of presence

R Functions (SQL Data Download)

To assist with my regular tasks accessing collar and capture data, I wrote R code to easily access those SQL databases.

## Download Collar Locations from SQL Database (filtered by species)
getSQLData_locs <- function(sqltable = "Collars_Locations_ALL", 
                            species){

  species <- match.arg(species,
                       c("Elk","Mule Deer","Pronghorn Antelope","Moose","Wolf","Black Bear","White-tailed Deer",
                         "Mountain Goat","Rocky Mountain Bighorn Sheep","California Bighorn Sheep","Mountain Lion"),
                       several.ok = T)
  con <- DBI::dbConnect(
    odbc::odbc(),
    driver = "SQL Server",
    database = "IFWIS_WildlifeReporting",
    uid = "CollarManager",
    pwd = #REDACTED,
    server = #REDACTED,
    port = #REDACTED)
  locdata <- DBI::dbReadTable(con,sqltable)
  DBI::dbDisconnect(con)
  sel <- dplyr::filter(locdata,Game %in% species)
  rm(con,sqltable,locdata)
  return(sel)
}

## Download Necropsy data from SQL Database (filtered by species)
getSQLData_nec = function(sqltable = "SAMM_Necropsy", species) {
  species <- match.arg(species,
                       c("Elk","Mule Deer","Pronghorn Antelope","Moose","Wolf","Black Bear","White-tailed Deer",
                         "Mountain Goat","Rocky Mountain Bighorn Sheep","California Bighorn Sheep","Mountain Lion"),
                       several.ok = T)
  con <- DBI::dbConnect(
    odbc::odbc(),
    driver = "SQL Server",
    database = "IFWIS_WildlifeReporting",
    uid =  "ShinyUserInternal",
    pwd = #REDACTED,
    server = #REDACTED,
    port = #REDACTED)
  
  columns = c("CaptureID","Animal_ID","GameID","FateID","FateDate",
              "CensorID","CensorDate","RadFreq","Latitude","Longitude",
              "Bnumber","LastLoc","FateDate_OG","CensorDate_OG","NecID",
              "BGMR","BlueEarTag")
  
  queries = paste0("SELECT ", 
                   columns, 
                   " from ", sqltable)
  
  res = lapply(queries, function(query) {
    res = DBI::dbGetQuery(con, query)
    res
  })
  res = dplyr::bind_cols(res)
  
  lut_game <- DBI::dbGetQuery(con, "SELECT GAMEID, Game from GAME_PIC_GAME") %>%
    dplyr::filter(Game %in% species) %>%
    dplyr::select(GAMEID) %>%
    unique() %>%
    unlist()
  
  lut_fate <- DBI::dbGetQuery(con, "SELECT FateID, FateDesc from SAMM_PIC_FateTypes")

  sel <- dplyr::filter(res,GameID %in% lut_game) %>%
    dplyr::left_join(.,lut_fate,by="FateID")
  
  rm(con,sqltable,res,queries,columns,lut_game,lut_fate)
  
  return(sel)
}

## Download Capture data from SQL Database (filtered by species)
getSQLData_cap = function(sqltable = "SAMM_Capture", species) {
  library(dplyr)
  library(dbplyr)
  species <- match.arg(species,
                       c("Elk","Mule Deer","Pronghorn Antelope","Moose","Wolf","Black Bear","White-tailed Deer",
                         "Mountain Goat","Rocky Mountain Bighorn Sheep","California Bighorn Sheep","Mountain Lion"),
                       several.ok = T)
  con <- DBI::dbConnect(
    odbc::odbc(),
    driver = "SQL Server",
    database = "IFWIS_WildlifeReporting",
    uid =  "ShinyUserInternal",
    pwd = #REDACTED,
    server = #REDACTED,
    port = #REDACTED)
  
  tab = tbl(con,sqltable)
  columns = colnames(tab)
  queries = paste0("SELECT ", 
                   columns, 
                   " from ", sqltable)
  res = lapply(queries, function(query) {
    res = DBI::dbGetQuery(con, query)
    res
  })
  res = dplyr::bind_cols(res)
  
  lut <- DBI::dbGetQuery(con, "SELECT GAMEID, Game from GAME_PIC_GAME") %>%
    dplyr::filter(Game %in% species) %>%
    dplyr::select(GAMEID) %>%
    unique() %>%
    unlist()
  
  sel <- dplyr::filter(res,GameID %in% lut)
  
  rm(con,sqltable,res,queries,columns,lut)
  
  return(sel)
}
#cap <- getSQLData_cap(species = "Mule Deer")

## Download Vegetation field data from SQL Database 
getSQLData_veg <- function(sqltable = "Veg_fsvm_understory_model_data"){
  con <- DBI::dbConnect(
    odbc::odbc(),
    driver = "SQL Server",
    database = "IFWIS_WildlifeReporting",
    uid = "CollarManager",
    pwd = #REDACTED,
    server = #REDACTED,
    port = #REDACTED)
  fielddata <- DBI::dbReadTable(con,sqltable)
  DBI::dbDisconnect(con); rm(con,sqltable)
  return(fielddata)
}

Northern Idaho Ground Squirrel Covariates

My expertise with using R for GIS manipulations has been previously called upon to assist with covariate formatting for northern Idaho ground squirrel habitat modeling efforts.

# Load packages
source('B:/Seasonal Range Analysis/Mule Deer/Smokey-Boise/winter/RScripts/instaload_function.R')
instaload(c('sf','terra','exactextractr'))

# Load Grid Shapefile
fp <- "Q:/RegionMcc/Diane EM/NIDGS/Covariates/2022_GISLayers_for_covariates/grid 100m_entire range_edited 2018 for Tamarack S_2021 S3 Revision.shp"
grid_shp <- sf::st_read(fp)

# List Veg Class Rasters
hveg_list <- dir("B:/Seasonal Range Analysis/Mule Deer/Covariates/hveg",pattern = "_2020.tif$",full.names = T)

# Load, Clamp, and Stack Rasters
## Function
load_clamp <- function(x){
  r <- terra::rast(x) %>% terra::clamp(.,lower=1,values=F)
  return(r)
}
hveg_stack <- do.call(load_clamp,list(hveg_list))

# ReProject grid shape (match raster stack)
grid_proj <- sf::st_transform(grid_shp,terra::crs(hveg_stack))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(hveg_stack,
                                             grid_proj,
                                             coverage_area = T,
                                             fun="sum",
                                             force_df =T)

# Rename columns
colnames(grid_extract) <- gsub("sum.","",colnames(grid_extract))
colnames(grid_extract) <- paste0(colnames(grid_extract),"_m2")

# Add to original
grid_final <- cbind(grid_shp,grid_extract)

# Load Canopy Cover Raster
cancov <- terra::rast("C:/Users/rritson/Documents/Projects/forErin/treecc30.tif")

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_shp,terra::crs(cancov))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(cancov,
                                             grid_proj,
                                             coverage_area = T,
                                             fun="mean",
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "NLCD_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)

#Save New Grid Shape
sf::st_write(grid_final,dsn = "Q:/RegionMcc/Diane EM/NIDGS/Covariates/2022_GISLayers_for_covariates/grid 100m_entire range_edited 2018 for Tamarack S_2021 S3 Revision_covs.shp")

#clean-up 
terra::tmpFiles(remove = T)
rm(cancov,hveg_list,hveg_stack,grid_shp,grid_proj,grid_extract,grid_final,fp,load_clamp,instaload)
gc()

###
# Load Grid
grid_final <- sf::st_read("Q:/RegionMcc/Diane EM/NIDGS/Covariates/2022_GISLayers_for_covariates/grid 100m_entire range_edited 2018 for Tamarack S_2021 S3 Revision_covs.shp")

# Load remaining covariates
## Soil Bulk Density
bulkdens <- terra::rast("Q:/RegionMcc/Diane EM/NIDGS/Covariates/EnvironmentalData/Soils/POLARIS/bulkdens")

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_final,terra::crs(bulkdens))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(bulkdens,
                                             grid_proj,
                                             coverage_area = T,
                                             fun="mean",
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "BulkDensity_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)

## Percent Silt
siltperc <- terra::rast("Q:/RegionMcc/Diane EM/NIDGS/Covariates/EnvironmentalData/Soils/POLARIS/siltperc")

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_final,terra::crs(siltperc))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(siltperc,
                                             grid_proj,
                                             coverage_area = T,
                                             fun="mean",
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "PercSilt_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)

## Soil Depth
resdep <- terra::rast("Q:/RegionMcc/Diane EM/NIDGS/Covariates/EnvironmentalData/Soils/POLARIS/resdep")

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_final,terra::crs(resdep))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(resdep,
                                             grid_proj,
                                             coverage_area = T,
                                             fun="mean",
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "SoilDepth_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)


#install.packages('spatialEco')
require(spatialEco)
dem <- terra::rast("C:/Users/rritson/Documents/Projects/forErin/DEM10M_1.tif")
## CLIP IS TOO SMALL!!!!!

## Aspect
asp <- terra::terrain(dem,"aspect")

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_final,terra::crs(asp))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(asp,
                                             grid_proj,
                                             coverage_area = T,
                                             fun='mean',
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "Aspect_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)

## Heat Load Index
hli <- spatialEco::hli(dem)

# ReProject grid shape (match raster)
grid_proj <- sf::st_transform(grid_final,terra::crs(hli))

# Extract coverage areas
grid_extract <- exactextractr::exact_extract(hli,
                                             grid_proj,
                                             coverage_area = T,
                                             fun='mean',
                                             force_df =T)
# Rename column
colnames(grid_extract) <- "HEAT_mean"

# Add to original
grid_final <- cbind(grid_final,grid_extract)

grid_final$AspectClass <- ifelse(23 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 68,"Northeast",
                                 ifelse(69 <= grid_final$Aspect_mean  & grid_final$Aspect_mean<= 112,"East",
                                        ifelse(113 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 158,"Southeast",
                                               ifelse(159 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 202,"South",
                                                      ifelse(203 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 248,"Southwest",
                                                             ifelse(249 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 292,"West",
                                                                    ifelse(293 <= grid_final$Aspect_mean & grid_final$Aspect_mean <= 338,"Northwest","North")))))))

#Save New Grid Shape
sf::st_write(grid_final,dsn = "Q:/RegionMcc/Diane EM/NIDGS/Covariates/2022_GISLayers_for_covariates/grid 100m_entire range_edited 2018 for Tamarack S_2021 S3 Revision_covs2.shp")

#clean-up 
terra::tmpFiles(remove = T)
rm(grid_proj,grid_extract,grid_final,asp,dem,resdep,siltperc,bulkdens,hli)
gc()

MTBS Fire Covariates

Recently, I developed code for manipulating fire data from MTBS into covariate rasters relevant to the FSVM and ungulate seasonal range modeling efforts.

# Load packages
require(sf)
require(terra)
require(foreign)
require(dplyr)
require(stringr)
require(lubridate)

### Set-up ------------------------------------------
## significant digits (7 digit coordinate + 10 digits after decimal (for nanometers) = 17 significant digits)
if(options()$digits != 17){
  options(digits = 17)
}
## 'terra' options
#'INT2S'
terra::terraOptions(datatype = 'INT2S') # same data type Brendan used, default in 'terra' is 'FLT4S'
terra::terraOptions(tolerance = 1e-10)
terra::terraOptions()
## proj4strings
#WGS84 Longitude and Latitude
lonlat = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#Idaho Transverse Mercator
#idtm = "+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
idtm <- "epsg:8826"

#Load base grid
base <- terra::rast('B:/Covariates/base_grid/idaho_buffer_100km_30m.tif')
terra::res(base)
terra::crs(base)
base <- terra::project(base, idtm)

#Load Idaho state boundary
idaho <- sf::st_read("B:/GIS_ReferenceLayers/ID_StateBoundary/IDTM", "IdahoBoundary_IDTM")  %>% 
  terra::vect(.) %>%
  terra::project(.,base) #ensure projections match
terra::plot(terra::ext(base))
terra::plot(idaho, add = T)

#Create bounding box from base grid
bbox <- base %>% 
  sf::st_bbox(.,crs = sf::st_crs(base)) %>% 
  sf::st_as_sfc(.) %>% 
  sf::st_as_sf(.) %>% 
  dplyr::mutate(PolyID = "NDVI_30m_Idaho_100kmBuff") %>%
  terra::vect(.)

# Load MTBS perimeter file
## Download at: https://mtbs.gov/direct-download; Burned Areas Boundaries Dataset
fire <- sf::st_read("C:/Users/rritson/Documents/Covariates/Fire/MTBS_perimeter_shape/mtbs_perims_DD.shp",quiet=T) %>%
  sf::st_transform(.,idtm) %>%
  dplyr::select(Event_ID,Incid_Name,Incid_Type,BurnBndAc,Ig_Date,Comment) %>%
  dplyr::rowwise(.) %>%
  dplyr::mutate(Year = stringr::str_split_fixed(Ig_Date,"-",3)[[1]],
                Month = stringr::str_split_fixed(Ig_Date,"-",3)[[2]],
                Day = stringr::str_split_fixed(Ig_Date,"-",3)[[3]]) %>%
  dplyr::select(Event_ID,Incid_Name,Incid_Type,BurnBndAc,Year,Month,Day,Comment) %>%
  terra::vect(.) %>%
  terra::crop(.,bbox) %>%
  sf::st_as_sf(.) #Important for filtering!!!

# Create Binary Burned Rasters
## Filter by Year
yrs <- 2000:2021
for(yr in yrs){
  print(paste0("Calculating ",yr,"..."))
  ### since 1984
  print(paste("All Years"))
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="max",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "Burned"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/Burned_allyrs_",yr,".tif"))
  
  ### Last 10 yrs
  print(paste("Last Ten Years"))
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr & Year >= (yr-10)) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="max",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "Burned"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/Burned_last10_",yr,".tif"))
  
  ### Only Year of
  print(paste("Year of Only"))
  fire_yr <- fire %>%
    dplyr::filter(Year == yr) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="max",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "Burned"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/Burned_yrof_",yr,".tif"))
}

## Frequency of Fire (N unique)
yrs <- 2000:2021
for(yr in yrs){
  print(paste0("Calculating ",yr,"..."))
  ### since 1984
  print(paste("All Years"))
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="sum",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "N_fires"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/FireFrequency_allyrs_",yr,".tif"))
  
  ### Last 10 yrs
  print(paste("Last Ten Years"))
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr & Year >= (yr-10)) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="sum",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "N_fires"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/FireFrequency_last10_",yr,".tif"))
  
  ### Only Year of
  print(paste("Year of Only"))
  fire_yr <- fire %>%
    dplyr::filter(Year == yr) %>%
    dplyr::mutate(Burned = 1) %>%
    dplyr::select(Burned) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,fun="sum",background=0,touches=F,na.rm=T,update=F,by=NULL,cover=F)
  names(fire_rast) <- "N_fires"
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/FireFrequency_yrof_",yr,".tif"))
}

## Time since last fire (Year)
yrs <- 2000:2021
for(yr in yrs){
  print(paste0("Calculating ",yr,"..."))
  
  ## since 1984
  print(paste("All Years"))
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr) %>%
    #dplyr::rowwise(.) %>%
    dplyr::mutate(TSF = yr-as.numeric(Year)) %>%
    dplyr::select(TSF,geometry) #%>%
    #dplyr::ungroup(.) %>%
    #terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,"TSF",fun=min,background=-9999,touches=F,na.rm=T,update=F,by=NULL,cover=F) #9999 never burned or not burned within defined period
  names(fire_rast) <- "TSF" #time since last fire (years)
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/YearSinceFire_allyrs_",yr,".tif"))

  ## Last 10 years
  print(paste("Last Ten Years")) 
  fire_yr <- fire %>%
    dplyr::filter(Year <= yr & Year >= (yr-10)) %>%
    dplyr::rowwise(.) %>%
    dplyr::mutate(TSF = yr-as.numeric(Year)) %>%
    dplyr::select(TSF) %>%
    terra::vect(.)
  fire_rast <- terra::rasterize(fire_yr,base,"TSF",fun=min,background=-9999,touches=F,na.rm=T,update=F,by=NULL,cover=F) #9999 never burned or not burned within defined period
  names(fire_rast) <- "TSF" #time since last fire (years)
  terra::writeRaster(fire_rast,paste0("E:/Seasonal_range_covars/fire/YearSinceFire_last10_",yr,".tif"),overwrite=T)
}

Region 5

Being stationed in the regional office, I occasionally assist staff with GIS-related questions and tasks. I have also helped out with field tasks including mule deer drive nets, pronghorn surveys, and fisheries.

Sterling WMA Spatially Balanced Random Sampling

Generalized random tessellation stratified (GRTS) algorithm

Blackfoot Reservoir Lowland Lake Survey Sites

R code: fishnetR

fishnetR <- function(shp,cell_size,n,stratify=TRUE,seed=1){
  require(dplyr)
  
  # Create grid (fishnet)
  grid <- sf::st_make_grid(shp, cellsize = cell_size, what = "polygons", square = T) %>% 
    terra::vect(.) %>%
    terra::crop(.,terra::vect(shp)) %>% 
    sf::st_as_sf(.) %>%
    sf::st_cast(.,"POLYGON") %>%
    dplyr::mutate(id = seq(1,nrow(.), by=1))
  
  #Calculate Stratas
  stratas <- rep(c(1:n),each=floor(nrow(grid)/n))
  
  # randomly assign remainder
  extra <- sample.int(n,nrow(grid)%%n,F)
  
  # combine strata
  stratas <- sort(c(stratas,extra))
  
  # Set stratas
  grid <- dplyr::mutate(grid,strata = as.factor(stratas))
    
  #Stratify Random Sample
  set.seed(seed)
  if(stratify){
    i=sample(1:1000,1)
    repeat{
      i= i+1
      set.seed(i)
      spx <- sf::st_drop_geometry(grid)
      if(!inherits(spx[,"strata"], "factor")) 
        spx[,"strata"] <- factor(spx[,"strata"]) 
      spx$REP <- NA
      reps=1
      nn=1
      results <- list()
      for(j in levels(spx[, "strata"])) {
        d <- spx[spx[,"strata"] == j,]
        d$rowname <- rownames(d)
        if(nrow(d) > n) {   
          for (i in 1:reps) {   
            s <- lapply(1, function(ij) {
              d[sample(1:nrow(d), nn),]})
            s[[1]]$REP <- i
            results[[paste(j,i,sep="_")]] <-s[[1]] 
          }
        } else {
          d$REP <- 1
          results[[paste(j,i,sep="_")]] <- d
        }
      }
      results
      results <- do.call(rbind, results)
      replace=F
      if(!replace){
        if(any(duplicated(results$rowname))){
          results <- results[-which(duplicated(results$rowname)),]
        }
      }
      results <- stats::na.omit(results[,c("rowname","REP")])
      results <- merge(grid, results, by.y="rowname", by.x = 'row.names', 
                       all.x = FALSE, all.y = TRUE)
      strat_rand_samp <- results %>% 
        sf::st_as_sf(.)
      if(any(sf::st_relate(strat_rand_samp, pattern = "F***1****",sparse=F)==T)==F && nrow(strat_rand_samp)==21){
        break
      }
    }
    #return(strat_rand_samp)
    return(list(strat_rand_samp = strat_rand_samp, grid = grid))
    
    #Simple Random Sample
  }else{
    randsamp <- sample.int(nrow(grid), size = n, replace=F)
    locs <- grid[randsamp,]
    #sf::st_write(locs,filepath) #write grids to filepath (forthcoming)
    return(locs)
  }
}