R/ceip_regression.r

In khufkens/ceipr: Centre on Emission Inventories and Projections processing code

#' Run a simple linear regression on the
#' the Centre on Emission Inventories and Projections (CEIP)
#' emission data as formatted using ceip_download()
#'
#' @param file path to a file generated by ceip_download()
#' @param out_dir output directory
#' @param internal logical TRUE / FALSE
#' @return geotiff written to disk
#' @keywords emission, voc
#' @export
#' @examples
#'
#' \dontrun{
#' ceip_download()
#'}

ceip_regression <- function(file = "~/SOx_A.tif",
                            out_dir = "~",
                            internal = TRUE){
  # read in the data
  s <- raster::brick(file)

  # for statistical consistency set seed
  set.seed(0)

  # regression of values in one brick (or stack) with 'time'
  time <- 1:raster::nlayers(s)

  # extract regression parameters
  slope_fun <- function(x) {
    if(all(is.na(x))){
      rep(NA,3)
    } else {
      fit <- stats::lm(x ~ time)
      slope <- fit$coefficients[2]

      p <- try(summary(fit)$coefficients[2,4], silent = TRUE)
      if(inherits(p, "try-error")){
        p <- NA
      }

      r2 <- try(summary(fit)$r.squared, silent = TRUE)
      if(inherits(r2, "try-error")){
        r2 <- NA
      }

      # return data
      return(as.numeric(c(slope, r2, p)))
    }
  }

  # calculate statistics
  fit_data <- raster::calc(s, slope_fun)

  # combine data, assign layer names
  names(fit_data) <- c("slope","r_squared", "p_value")

  # return a raster stack with regression
  # parameters
  if(internal){
    return(fit_data)
  } else {
    raster::writeRaster(fit_data,
                        paste0(tools::file_path_sans_ext(file),"_regression.tif"),
                        overwrite = TRUE,
                        options = "COMPRESS=DEFLATE")
    }
}