R/GetBioVars.R

In SMBC-NZP/BBSclim: Analysis of BBS data using correlated occupancy models to examine relationship between range dynamics and climate variables

#' GetBioVars
#'
#' Extract scaled annual bioclim variable values for each BBS route
#' @param count Data frame containing the (buffered) count data for the focal species
#' @param index Integer vector containing the bioclim variables of interest
#' @param ind_name Character vector containing abbreviated name of the bioclim variables of interest
#' @return A .csv files containing the following fields:
#' @return   routeID The unique 8 digit route ID for each route
#' @return   Latitude The latitude for the route
#' @return   Longitude The longitude for the route
#' @return   Ind_Year The value for biovar[Ind] in year Year
#' @export

GetBioVars <- function(alpha, index = c(1, 2, 8, 12, 18),
                       ind_name = c("tmp", "dtr", "Twet", "Prec", "Pwarm")){

  counts <- read.csv(paste0("inst/output/", alpha, "/count_buff.csv"))

  rxy <- dplyr::select(counts, routeID, Longitude, Latitude)
  rxy <- rxy[!duplicated(rxy),]

  xy <- dplyr::select(rxy, Longitude, Latitude)

  start_yr <- min(counts$Year)
  end_yr <- max(counts$Year)
  bbs_years <- seq(from = start_yr, to = end_yr)

  if(sum(bbs_years %in% seq(from = 1971, to = 2014)) < length(bbs_years)) stop("Count data contains years with no climate data")

  for (jj in 1:length(index)) {
    for (ii in 1:(end_yr - start_yr + 1)) {
      out	<- raster::extract(NA_biovars[[ii]][[index[jj]]], xy)
      varname <- paste(ind_name[[jj]], bbs_years[[ii]], sep = "_")
      rxy[[varname]] <- out
    }
  }

  ## For locations w/ NA for climate variables, fill in mean of 8 neighboring cells
  problem_routes <- which(is.na(rxy[, ncol(rxy)]))
  first_climate <- grep(ind_name[1], colnames(rxy))[1]

  if(length(problem_routes) > 0){
    ind <- 0
    while(ind < 1){
      problem_xy <- dplyr::select(dplyr::slice(rxy, problem_routes), Longitude, Latitude)

      fix_routes <- problem_out <- problem_nearby <- NULL
      adj <- c(-.5,0,.5)

      for (jj in 1:length(index)) {
        for (ii in 1:(end_yr - start_yr + 1)) {
          for (i_lat in 1:3) {
            for (i_lon in 1:3) {
              mjc	<- raster::extract(NA_biovars[[ii]][[index[jj]]],
                                     problem_xy + matrix(c(adj[i_lon], adj[i_lat]),
                                                         nrow = length(problem_routes), ncol=2), byrow=T)
              problem_nearby	<- cbind(problem_nearby, mjc)
            }
          }
          problem_out <- rowMeans(problem_nearby, na.rm = TRUE)
          fix_routes <- cbind(fix_routes, problem_out)
        }
      }
      rxy[problem_routes, first_climate:ncol(rxy)] <- fix_routes
      problem_routes2 <- which(is.na(rxy[, ncol(rxy)]))
      if(length(problem_routes2) == 0 | length(problem_routes2) == length(problem_routes)){
        ind <- 1
        problem_routes <- problem_routes2
      }else{
        problem_routes <- problem_routes2
      }
    }
    if(length(problem_routes) > 0){
      ind <- 0
      min.adj <- -1
      max.adj <- 1
      while(ind < 1){
        problem_xy <- dplyr::select(dplyr::slice(rxy, problem_routes), Longitude, Latitude)

        fix_routes <- problem_out <- problem_nearby <- NULL
        adj <- c(min.adj,0,max.adj)

        for (jj in 1:length(index)) {
          for (ii in 1:(end_yr - start_yr + 1)) {
            for (i_lat in 1:3) {
              for (i_lon in 1:3) {
                mjc	<- raster::extract(NA_biovars[[ii]][[index[jj]]],
                                       problem_xy + matrix(c(adj[i_lon], adj[i_lat]),
                                                           nrow = length(problem_routes), ncol=2), byrow=T)
                problem_nearby	<- cbind(problem_nearby, mjc)
              }
            }
            problem_out <- rowMeans(problem_nearby, na.rm = TRUE)
            fix_routes <- cbind(fix_routes, problem_out)
          }
        }
        rxy[problem_routes, first_climate:ncol(rxy)] <- fix_routes
        problem_routes2 <- which(is.na(rxy[, ncol(rxy)]))
        if(length(problem_routes2) == 0 | length(problem_routes2) == length(problem_routes)){
          ind <- 1
          problem_routes <- problem_routes2
        }else{
          problem_routes <- problem_routes2
          min.adj <- min.adj - 0.5
          max.adj <- max.adj + 0.5
        }
      }
    }
  }


  ### Center and scale climate variables
  clim_scale <- matrix(99, nrow = length(ind_name), ncol = 2, dimnames = list(ind_name, c("mean", "sd")))

  for(ii in seq_along(ind_name)){
    clim_scale[ii, "mean"] <- mean(as.matrix(dplyr::select(rxy, grep(ind_name[ii], names(rxy)))), na.rm = TRUE)
    clim_scale[ii, "sd"] <- sd(as.matrix(dplyr::select(rxy, grep(ind_name[ii], names(rxy)))), na.rm = TRUE)

    rxy[,grep(ind_name[ii], names(rxy))] <- (rxy[,grep(ind_name[ii], names(rxy))] - clim_scale[ii, "mean"]) / clim_scale[ii, "sd"]
  }

  write.csv(clim_scale,
            paste0("inst/output/", alpha, "/clim_scale.csv"),
            row.names = FALSE)

  ### Add squared climate variables
  sq_rxy <- dplyr::as_data_frame(rxy[, first_climate:ncol(rxy)] ^ 2)
  colnames(sq_rxy) <- paste("sq", colnames(sq_rxy), sep = "_")
  rxy <- dplyr::bind_cols(rxy, sq_rxy)

  write.csv(rxy, paste0("inst/output/", alpha, "/route_clim.csv"),
            row.names = FALSE)
}