R/vg_calc.R

In jnghiem/bfasttools: Tools for BFAST spatial analysis

#' Perform variogram analysis on \code{bfastSpatial()} shift and slope
#' statistics
#'
#' This function uses the \code{gstat} package to compute variogram statistic
#' for the shift and slope output from \code{bfastSpatial} output. This function
#' uses the exponential variogram model.
#'
#' This function returns a data frame with the following fields: \enumerate{
#' \item \code{year}: The year of the observation. \item \code{month}: The month
#' of the observations. \item \code{stat}: The BFAST statistics (\code{"slope"}
#' or \code{"shift"}) \item \code{model}: The variogram model (\code{"Nug"} or
#' \code{"Exp"}) \item \code{psill}: The partial sill. \item \code{range}: The
#' range parameter of the model. Only valid for Exp model. \item
#' \code{effective_range}: The effective range of the model (within 5% of the
#' sill). Only valid for Exp model. \item \code{SSerr}: The sum of squared
#' errors of the model from the sample variogram. \item \code{n}: The number of
#' points in the sample variogram. }
#'
#' If the variogram model is singular or does not converge after 200 iterations,
#' no model is fit and \code{NA}s are returned for the entire row. Since this
#' function tries to fit a large number of models, it may fail to fit a
#' substantial proportion of models if the sample variograms do not match the
#' model well (or require transformation).
#'
#' @param bfast_in A data frame generated by \code{bfastSpatial()}.
#' @param boundary A SpatialPolygons object defining the boundary of the
#'   variogram analysis. If this argument is not supplied, all records in
#'   \code{bfast_in} will be used.
#' @param template A Raster object with the same resolution, extent, and
#'   projection as the raster dataset used to generate \code{bfast_in}.
#'   Alternatively, an XML file generated by \code{create_raster_metadata()}.
#' @param cutoff_dist An optional data frame containing the output of
#'   \code{pcf_calc()} for the same study area.
#' @param mc.cores A numeric indicating the number of cores to be used in
#'   parallel computing.
#' @return A data frame with variogram results.
#' @examples
#' \dontrun{
#' vg_calc(bf_df, boundary=rgdal::readOGR(dsn="C:/Desktop/shapefiles", layer="Mojave"), template="C:/Desktop/raster_metadata.xml", cutoff_dist=pcf_results, mc.cores=6)
#' }
#' @importFrom lubridate year month
#' @importFrom foreach %dopar% foreach
#' @importFrom evaluate evaluate
#' @importFrom dplyr arrange filter
#' @export
vg_calc <- function(bfast_in, boundary=NULL, template, cutoff_dist=NULL, mc.cores=6) {
  template <- create_template_raster(template)
  if (!is.null(boundary)) {
    cells <- extract(template, boundary, cellnumbers=TRUE, df=TRUE)[,"cell"]
  } else {
    cells <- 1:ncell(template)
  }
  brks_only <- bfast_in %>%
    clean_bfast() %>%
    filter(no_cell %in% cells, brk==1)
  years <- brks_only[,"year"]
  yr.range <- min(years):max(years)
  stat <- c("shift", "slope")
  if (mc.cores>1) {
    cl <- snow::makeCluster(mc.cores, type="SOCK")
    doSNOW::registerDoSNOW(cl)
    variogram_stat <- foreach(i=yr.range, .combine=rbind, .inorder=FALSE, .packages=c("evaluate", "gstat", "raster", "dplyr"), .export="create_raster") %dopar% {
      add.df <- data.frame()
      for (j in 1:12) {
        brk_stat <- filter(brks_only, year==i, month==j)
        if (nrow(brk_stat)==0) next
        for (k in 1:length(stat)) {
          stat_pts <- rasterToPoints(create_raster("no_cell", stat[k], brk_stat,  template), spatial=TRUE)
          cut <- filter(cutoff_dist, year==i, month==j, type=="first peak")[,"dist"]*0.8
          if (length(cut)==0 | is.na(cut)) {
            vg.samp <- variogram(layer~x+y, data=stat_pts)
          } else {
            vg.samp <- variogram(layer~x+y, data=stat_pts, cutoff=cut)
            if (is.null(vg.samp)) vg.samp <- variogram(layer~x+y, data=stat_pts)
          }
          if (!is.null(vg.samp)) {
            test.mod <- evaluate(fit.variogram(vg.samp, vgm("Exp")))
            vg.mod <- fit.variogram(vg.samp, vgm("Exp"))
          }
          if (any(grepl("No convergence", sapply(test.mod, as.character))) | attr(vg.mod, "singular") | is.null(vg.samp)) {
            new_rows <- data.frame(year=i, month=j, stat=NA, model=NA, psill=NA, range=NA, effective_range=NA, SSerr=NA, n=NA)
          } else {
            new_rows <- data.frame(year=i, month=j, stat=stat[k], dplyr::select(as.data.frame(vg.mod), model, psill, range), effective_range=c(0, effective_range(vg.mod)), SSerr=c(0, attr(vg.mod, "SSErr")), n=nrow(vg.samp))
          }
          add.df <- rbind(add.df, new_rows)
        }
      }
      return(add.df)
    }
    snow::stopCluster(cl)
  } else {
    variogram_stat <- data.frame()
    for (i in yr.range) {
      for (j in 1:12) {
        brk_stat <- filter(brks_only, year==i, month==j)
        if (nrow(brk_stat)==0) next
        for (k in 1:length(stat)) {
          stat_pts <- rasterToPoints(create_raster("no_cell", stat[k], brk_stat, template), spatial=TRUE)
          cut <- filter(cutoff_dist, year==i, month==j, type=="first peak")[,"dist"]*0.8
          if (length(cut)==0 | is.na(cut)) {
            vg.samp <- variogram(layer~x+y, data=stat_pts)
          } else {
            vg.samp <- variogram(layer~x+y, data=stat_pts, cutoff=cut)
            if (is.null(vg.samp)) vg.samp <- variogram(layer~x+y, data=stat_pts)
          }
          if (!is.null(vg.samp)) {
            test.mod <- evaluate(fit.variogram(vg.samp, vgm("Exp")))
            vg.mod <- fit.variogram(vg.samp, vgm("Exp"))
          }
          if (any(grepl("No convergence", sapply(test.mod, as.character))) | attr(vg.mod, "singular") | is.null(vg.samp)) {
            new_rows <- data.frame(year=i, month=j, stat=NA, model=NA, psill=NA, range=NA, effective_range=NA, SSerr=NA, n=NA)
          } else {
            new_rows <- data.frame(year=i, month=j, stat=stat[k], dplyr::select(as.data.frame(vg.mod), model, psill, range), effective_range=c(0, effective_range(vg.mod)), SSerr=c(0, attr(vg.mod, "SSErr")), n=nrow(vg.samp))
          }
          variogram_stat <- rbind(variogram_stat, new_rows)
        }
      }
    }
  }
  return(arrange(variogram_stat, year, month))
}