diff --git a/.Rbuildignore b/.Rbuildignore index e24c06a..2e6b4a4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,4 @@ ^\README.Rmd$ ^dev_history\.R$ ^\.github$ +^data-raw$ \ No newline at end of file diff --git a/.gitignore b/.gitignore index e92156b..8c466ba 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ src/*.so src/*.dll README.html README.Rmd +data-raw/*.RData \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 2ad2f5d..698faaf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: remote Title: Empirical Orthogonal Teleconnections in R -Version: 1.2.3 -Authors@R: +Version: 1.2.3.9009 +Authors@R: c(person(given = "Tim", family = "Appelhans", role = c("cre", "aut"), @@ -25,24 +25,29 @@ Description: Empirical orthogonal teleconnections in R. space and time, EOT analysis produces patterns that are orthogonal in either space or time. License: GPL (>= 3) | file LICENSE -Depends: +Depends: Rcpp (>= 0.10.3), raster, - methods -Imports: + methods, + R (>= 3.5) +Imports: grDevices, gridExtra, latticeExtra, mapdata, scales, stats, + terra, utils -Suggests: - maps, - lattice, +Suggests: + checkmate, grid, - sp + lattice, + maps, + sp, + tinytest LinkingTo: Rcpp -RoxygenNote: 7.3.2 Roxygen: list(markdown = TRUE) Encoding: UTF-8 +Config/roxygen2/version: 8.0.0 +LazyData: true diff --git a/NAMESPACE b/NAMESPACE index 9b2be82..3bdb4fb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,17 +10,17 @@ export(denoise) export(deseason) export(geoWeight) export(getWeights) -export(lagalize) -export(longtermMeans) +export(nmodes) export(readEot) export(subset) exportClasses(EotMode) exportClasses(EotStack) exportMethods("names<-") exportMethods(eot) +exportMethods(lagalize) +exportMethods(longtermMeans) exportMethods(nXplain) exportMethods(names) -exportMethods(nmodes) exportMethods(plot) exportMethods(predict) exportMethods(writeEot) @@ -32,6 +32,7 @@ import(methods) import(parallel) import(raster) import(scales) +import(terra) importFrom(grDevices,colorRampPalette) importFrom(grDevices,hcl) importFrom(stats,cov.wt) diff --git a/NEWS.md b/NEWS.md index 90bc6d0..65b6477 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,18 @@ +# remote 1.2.3.9009 (2026-07-22) + +#### ✨ features and improvements + + * `anomalize()` fails gracefully if 'reference' is not a raster object, and + uses the 1st layer (with a warning) if a multi-layer object is supplied. + * `denoise()` fails gracefully if neither 'expl.var' nor 'k' are supplied. + +#### 🐛 bug fixes + +#### 💬 documentation etc + +#### 🍬 miscellaneous + + # remote 1.2.3 (2025-04-12) #### 💬 documentation etc diff --git a/R/AAAremoteClasses.R b/R/AAAremoteClasses.R index d9e108c..ad4cce5 100644 --- a/R/AAAremoteClasses.R +++ b/R/AAAremoteClasses.R @@ -6,54 +6,58 @@ #' @slot coords_bp the coordinates of the identified base point #' @slot cell_bp the cell number of the indeified base point #' @slot cum_exp_var the cumulative explained variance of the considered EOT mode -#' @slot r_predictor the RasterLayer of the correlation coefficients between the base point and each pixel of the predictor domain +#' @slot r_predictor `SpatRaster` of the correlation coefficients between the base point and each pixel of the predictor domain #' @slot rsq_predictor as above but for the coefficient of determination of the predictor domain #' @slot rsq_sums_predictor as above but for the sums of coefficient of determination of the predictor domain -#' @slot int_predictor the RasterLayer of the intercept of the regression equation for each pixel of the predictor domain -#' @slot slp_predictor same as above but for the slope of the regression equation for each pixel of the predictor domain -#' @slot p_predictor the RasterLayer of the significance (p-value) of the the regression equation for each pixel of the predictor domain -#' @slot resid_predictor the RasterBrick of the reduced data for the predictor domain -#' @slot r_response the RasterLayer of the correlation coefficients between the base point and each pixel of the response domain +#' @slot int_predictor `SpatRaster` of the intercept of the regression equation for each pixel of the predictor domain +#' @slot slp_predictor as above but for the slope of the regression equation for each pixel of the predictor domain +#' @slot p_predictor `SpatRaster` of the significance (p-value) of the regression equation for each pixel of the predictor domain +#' @slot resid_predictor `SpatRaster` of the reduced data for the predictor domain +#' @slot r_response `SpatRaster` of the correlation coefficients between the base point and each pixel of the response domain #' @slot rsq_response as above but for the coefficient of determination of the response domain -#' @slot int_response the RasterLayer of the intercept of the regression equation for each pixel of the response domain +#' @slot int_response `SpatRaster` of the intercept of the regression equation for each pixel of the response domain #' @slot slp_response as above but for the slope of the regression equation for each pixel of the response domain -#' @slot p_response same the RasterLayer of the significance (p-value) of the the regression equation for each pixel of the response domain -#' @slot resid_response the RasterBrick of the reduced data for the response domain +#' @slot p_response `SpatRaster` of the significance (p-value) of the regression equation for each pixel of the response domain +#' @slot resid_response `SpatRaster` of the reduced data for the response domain #' #' @exportClass EotMode #' @rdname EotMode-class -setClass('EotMode', - slots = c(mode = 'integer', - name = 'character', - eot = 'numeric', - coords_bp = 'matrix', - cell_bp = 'integer', - cum_exp_var = 'numeric', - r_predictor = 'RasterLayer', - rsq_predictor = 'RasterLayer', - rsq_sums_predictor = 'RasterLayer', - int_predictor = 'RasterLayer', - slp_predictor = 'RasterLayer', - p_predictor = 'RasterLayer', - resid_predictor = 'RasterBrick', - r_response = 'RasterLayer', - rsq_response = 'RasterLayer', - int_response = 'RasterLayer', - slp_response = 'RasterLayer', - p_response = 'RasterLayer', - resid_response = 'RasterBrick')) +methods::setClass( + 'EotMode' + , slots = c(mode = 'integer' + , name = 'character' + , eot = 'numeric' + , coords_bp = 'matrix' + , cell_bp = 'integer' + , cum_exp_var = 'numeric' + , r_predictor = 'SpatRaster' + , rsq_predictor = 'SpatRaster' + , rsq_sums_predictor = 'SpatRaster' + , int_predictor = 'SpatRaster' + , slp_predictor = 'SpatRaster' + , p_predictor = 'SpatRaster' + , resid_predictor = 'SpatRaster' + , r_response = 'SpatRaster' + , rsq_response = 'SpatRaster' + , int_response = 'SpatRaster' + , slp_response = 'SpatRaster' + , p_response = 'SpatRaster' + , resid_response = 'SpatRaster' + ) +) NULL #' Class EotStack +#' #' @slot modes a list containing the individual 'EotMode's of the 'EotStack' #' @slot names the names of the modes #' #' @exportClass EotStack #' @rdname EotStack-class -setClass('EotStack', +methods::setClass('EotStack', slots = c(modes = 'list', names = 'character')) diff --git a/R/EotCycle.R b/R/EotCycle.R index 882b04a..8acf2df 100644 --- a/R/EotCycle.R +++ b/R/EotCycle.R @@ -1,139 +1,157 @@ +methods::setGeneric( + "EotCycle" + , function(x, ...) { + standardGeneric("EotCycle") + } +) + #' Calculate a single EOT #' #' @description -#' EotCycle() calculates a single EOT and is controlled by the main eot() function +#' Calculates a single EOT and is controlled by the main [eot()] function. #' -#' @param x a ratser stack used as predictor -#' @param y a RasterStack used as response. If \code{y} is `NULL`, -#' \code{x} is used as \code{y} -#' @param n the number of EOT modes to calculate -#' @param standardised logical. If `FALSE` the calculated r-squared values -#' will be multiplied by the variance +#' @param x,y,n,standardised,write.out,path.out,prefix,type,verbose See [eot()]. #' @param orig.var original variance of the response domain -#' @param write.out logical. If `TRUE` results will be written to disk -#' using \code{path.out} -#' @param path.out the file path for writing results if \code{write.out} is `TRUE`. -#' Defaults to current working directory -#' @param prefix optional prefix to be used for naming of results if -#' \code{write.out} is `TRUE` -#' @param type the type of the link function. Defaults to \code{'rsq'} as in original -#' proposed method from \cite{Dool2000}. If set to \code{'ioa'} index of agreement is -#' used instead -#' @param verbose logical. If `TRUE` some details about the -#' calculation process will be output to the console -#' @param ... If \code{write.out = TRUE}, further arguments passed to -#' [writeEot()]. +#' @param ... Further arguments passed to [writeEot()] if `write.out = TRUE`, +#' or to the underlying `SpatRaster` method in general for `Raster*` input. #' #' @export EotCycle -EotCycle <- function(x, - y, - n = 1, - standardised, - orig.var, - write.out, - path.out, - prefix, - type, - verbose, - ...) { - - ### Identification of the most explanatory pred pixel - - # Extract pixel entries from RasterStack objects - x.vals <- raster::getValues(x) - y.vals <- raster::getValues(y) - type <- type[1] - - # Calculate and summarize R-squared per pred pixel - if (verbose) { - cat("\nCalculating linear model ...", "\n") - } - - type <- type[1] - if (type == "rsq") { - a <- predRsquaredSum(pred_vals = x.vals, resp_vals = y.vals, - standardised = standardised) - } else { - a <- iodaSumC(pred_vals = x.vals, resp_vals = y.vals) - } - - # Identify pred pixel with highest sum of r.squared - if (verbose) { - cat("Locating ", n, ". EOT ...", "\n", sep = "") - } - - maxxy.all <- which(a == max(a, na.rm = TRUE)) - maxxy <- maxxy.all[1] +#' @name EotCycle - if (length(maxxy.all) != 1) { - if (verbose) { - cat("WARNING:", "\n", - "LOCATION OF EOT AMBIGUOUS!", "\n", - "MULTIPLE POSSIBLE LOCATIONS DETECTED, USING ONLY THE FIRST!\n\n") - } - } - if (verbose) { - cat("Location:", raster::xyFromCell(x, maxxy), "\n", sep = " ") +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases EotCycle,RasterStackBrick-method +#' @rdname EotCycle +methods::setMethod( + "EotCycle" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + EotCycle( + terra::rast(x) + , ... + ) } - - ### Regression of most explanatory pred pixel with resp pixels - - ## Fit lm - - # lm(y.vals[i, ] ~ x.vals[maxxy, ]) with T statistics - y.lm.param.t <- respLmParam(x.vals, y.vals, maxxy - 1) # C++ starts at 0! - # Calculate p value from T statistics - y.lm.param.p <- lapply(y.lm.param.t, function(i) { - tmp <- i - tmp[[5]] <- 2 * pt(-abs(tmp[[5]]), df = tmp[[6]]) - - return(tmp) - }) - - - ## Rasterize lm parameters - - # RasterLayer template for R-squared, slope and p value - rst.y.template <- raster::raster(nrows = raster::nrow(y), - ncols = raster::ncol(y), - xmn = raster::xmin(y), - xmx = raster::xmax(y), - ymn = raster::ymin(y), - ymx = raster::ymax(y)) - - rst.y.r <- rst.y.rsq <- rst.y.intercept <- - rst.y.slp <- rst.y.p <- rst.y.template +) - # RasterBrick template for residuals - brck.y.resids <- raster::brick(nrows = raster::nrow(y), - ncols = raster::ncol(y), - xmn = raster::xmin(y), - xmx = raster::xmax(y), - ymn = raster::ymin(y), - ymx = raster::ymax(y), - nl = raster::nlayers(y)) - - # R - rst.y.r[] <- sapply(y.lm.param.p, "[[", 1) - # R-squared - rst.y.rsq[] <- sapply(y.lm.param.p, "[[", 1) ^ 2 - # Intercept - rst.y.intercept[] <- sapply(y.lm.param.p, "[[", 2) - # Slope - rst.y.slp[] <- sapply(y.lm.param.p, "[[", 3) - # P value - rst.y.p[] <- sapply(y.lm.param.p, "[[", 5) - # Residuals - brck.y.resids[] <- matrix(sapply(y.lm.param.p, "[[", 4), - ncol = raster::nlayers(x), byrow = TRUE) - # EOT over time - eot.ts <- as.numeric(raster::extract(x, maxxy)[1, ]) - - ### Regression of most explanatory pred pixel with pred pixels - - # Following code is only executed when pred and resp are not equal +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases EotCycle,SpatRaster-method +#' @rdname EotCycle +methods::setMethod( + "EotCycle" + , signature(x = "SpatRaster") + , function(x + , y + , n = 1 + , standardised + , orig.var + , write.out + , path.out + , prefix + , type = c("rsq", "ioa") + , verbose + , ... + ) { + + ### Identification of the most explanatory pred pixel + + # Extract pixel entries from raster objects + x.vals <- terra::values(x) + y.vals <- terra::values(y) + type <- match.arg(type) + + # Calculate and summarize R-squared per pred pixel + if (verbose) { + cat("\nCalculating linear model ...", "\n") + } + + if (type == "rsq") { + a <- predRsquaredSum( + pred_vals = x.vals + , resp_vals = y.vals + , standardised = standardised + ) + } else { + a <- iodaSumC( + pred_vals = x.vals + , resp_vals = y.vals + ) + } + + # Identify pred pixel with highest sum of r.squared + if (verbose) { + cat("Locating ", n, ". EOT ...", "\n", sep = "") + } + + maxxy.all <- which(a == max(a, na.rm = TRUE)) + maxxy <- maxxy.all[1] + + if (length(maxxy.all) != 1) { + if (verbose) { + cat("WARNING:", "\n", + "LOCATION OF EOT AMBIGUOUS!", "\n", + "MULTIPLE POSSIBLE LOCATIONS DETECTED, USING ONLY THE FIRST!\n\n") + } + } + + if (verbose) { + cat("Location:", terra::xyFromCell(x, maxxy), "\n", sep = " ") + } + + ### Regression of most explanatory pred pixel with resp pixels + + ## Fit lm + + # lm(y.vals[i, ] ~ x.vals[maxxy, ]) with T statistics + y.lm.param.t <- respLmParam(x.vals, y.vals, maxxy - 1) # C++ starts at 0! + # Calculate p value from T statistics + y.lm.param.p <- lapply( + y.lm.param.t + , function(i) { + tmp <- i + tmp[[5]] <- 2 * pt(-abs(tmp[[5]]), df = tmp[[6]]) + + return(tmp) + } + ) + + + ## Rasterize lm parameters + + # multi-layer raster template for residuals + brck.y.resids <- terra::setValues(y, values = NA_real_) + + # single-layer raster template for R-squared, slope and p value + rst.y.template <- brck.y.resids[[1L]] + + rst.y.r <- rst.y.rsq <- rst.y.intercept <- + rst.y.slp <- rst.y.p <- rst.y.template + + # R + rst.y.r[] <- sapply(y.lm.param.p, "[[", 1) + # R-squared + rst.y.rsq[] <- sapply(y.lm.param.p, "[[", 1) ^ 2 + # Intercept + rst.y.intercept[] <- sapply(y.lm.param.p, "[[", 2) + # Slope + rst.y.slp[] <- sapply(y.lm.param.p, "[[", 3) + # P value + rst.y.p[] <- sapply(y.lm.param.p, "[[", 5) + # Residuals + brck.y.resids[] <- matrix(sapply(y.lm.param.p, "[[", 4), + ncol = terra::nlyr(x), byrow = TRUE) + # EOT over time + eot.ts <- as.numeric(terra::extract(x, maxxy)[1, ]) + + ### Regression of most explanatory pred pixel with pred pixels + + # Following code is only executed when pred and resp are not equal + ## Fit lm # lm(x.vals[i, ] ~ x.vals[maxxy, ]) with T statistics @@ -149,26 +167,15 @@ EotCycle <- function(x, ## Rasterize lm parameters - # RasterLayer template for R-squared, slope and p value - rst.x.template <- raster::raster(nrows = raster::nrow(x), - ncols = raster::ncol(x), - xmn = raster::xmin(x), - xmx = raster::xmax(x), - ymn = raster::ymin(x), - ymx = raster::ymax(x)) + # multi-layer raster template for residuals + brck.x.resids <- terra::setValues(x, values = NA_real_) + + # single-layer raster template for R-squared, slope and p value + rst.x.template <- brck.x.resids[[1L]] rst.x.r <- rst.x.rsq <- rst.x.rsq.sums <- rst.x.intercept <- - rst.x.slp <- rst.x.p <- rst.x.template - - # RasterBrick template for residuals - brck.x.resids <- raster::brick(nrows = raster::nrow(x), - ncols = raster::ncol(x), - xmn = raster::xmin(x), - xmx = raster::xmax(x), - ymn = raster::ymin(x), - ymx = raster::ymax(x), - nl = raster::nlayers(x)) - + rst.x.slp <- rst.x.p <- rst.x.template + # R rst.x.r[] <- sapply(x.lm.param.p, "[[", 1) # R-squared @@ -183,103 +190,114 @@ EotCycle <- function(x, rst.x.p[] <- sapply(x.lm.param.p, "[[", 5) # Residuals brck.x.resids[] <- matrix(sapply(x.lm.param.p, "[[", 4), - ncol = raster::nlayers(x), byrow = TRUE) - -# #expl.var <- x[maxxy] / orig.var -# if (!standardised) { -# t <- mean(apply(getValues(brck.y.resids), 1, var, na.rm = TRUE), -# na.rm = TRUE) -# s <- mean(apply(getValues(brck.y.resids), 2, var, na.rm = TRUE), -# na.rm = TRUE) -# resid.var <- t + s -# } else { -# resid.var <- var(as.vector(getValues(brck.y.resids)), na.rm = TRUE) -# } - - resid.var <- calcVar(brck.y.resids, standardised = standardised) - - cum.expl.var <- (orig.var - resid.var) / orig.var - - if (verbose) { - cat("Cum. expl. variance (%):", cum.expl.var * 100, "\n", sep = " ") - } - - xy <- raster::xyFromCell(x, maxxy) - location.df <- as.data.frame(cbind(xy, paste("mode", - sprintf("%02.f", n), - sep = "_"), - cum.expl.var, - if (length(maxxy.all) != 1) - "ambiguous" else "ok", - maxxy), - stringsAsFactors = FALSE) - names(location.df) <- c("x", "y", "mode", "cum_expl_var", "comment", "cell_bp") - mode(location.df$x) <- "numeric" - mode(location.df$y) <- "numeric" - mode(location.df$cum_expl_var) <- "numeric" - mode(location.df$cell_bp) <- "integer" - + ncol = terra::nlyr(x), byrow = TRUE) + + # #expl.var <- x[maxxy] / orig.var + # if (!standardised) { + # t <- mean(apply(getValues(brck.y.resids), 1, var, na.rm = TRUE), + # na.rm = TRUE) + # s <- mean(apply(getValues(brck.y.resids), 2, var, na.rm = TRUE), + # na.rm = TRUE) + # resid.var <- t + s + # } else { + # resid.var <- var(as.vector(getValues(brck.y.resids)), na.rm = TRUE) + # } + + resid.var <- calcVar(brck.y.resids, standardised = standardised) + + cum.expl.var <- (orig.var - resid.var) / orig.var + + if (verbose) { + cat("Cum. expl. variance (%):", cum.expl.var * 100, "\n", sep = " ") + } + + xy <- terra::xyFromCell(x, maxxy) + location.df <- as.data.frame(cbind(xy, paste("mode", + sprintf("%02.f", n), + sep = "_"), + cum.expl.var, + if (length(maxxy.all) != 1) + "ambiguous" else "ok", + maxxy), + stringsAsFactors = FALSE) + names(location.df) <- c("x", "y", "mode", "cum_expl_var", "comment", "cell_bp") + mode(location.df$x) <- "numeric" + mode(location.df$y) <- "numeric" + mode(location.df$cum_expl_var) <- "numeric" + mode(location.df$cell_bp) <- "integer" + ### Output # Output returned by function - out <- new('EotMode', - mode = n, - name = paste("mode", sprintf("%02.f", n), sep = "_"), - eot = eot.ts, - coords_bp = xy, - cell_bp = maxxy, - cum_exp_var = cum.expl.var, - r_predictor = rst.x.r, - rsq_predictor = rst.x.rsq, - rsq_sums_predictor = rst.x.rsq.sums, - int_predictor = rst.x.intercept, - slp_predictor = rst.x.slp, - p_predictor = rst.x.p, - resid_predictor = brck.x.resids, - r_response = rst.y.r, - rsq_response = rst.y.rsq, - int_response = rst.y.intercept, - slp_response = rst.y.slp, - p_response = rst.y.p, - resid_response = brck.y.resids) - + out <- new('EotMode', + mode = n, + name = paste("mode", sprintf("%02.f", n), sep = "_"), + eot = eot.ts, + coords_bp = xy, + cell_bp = maxxy, + cum_exp_var = cum.expl.var, + r_predictor = rst.x.r, + rsq_predictor = rst.x.rsq, + rsq_sums_predictor = rst.x.rsq.sums, + int_predictor = rst.x.intercept, + slp_predictor = rst.x.slp, + p_predictor = rst.x.p, + resid_predictor = brck.x.resids, + r_response = rst.y.r, + rsq_response = rst.y.rsq, + int_response = rst.y.intercept, + slp_response = rst.y.slp, + p_response = rst.y.p, + resid_response = brck.y.resids) + # Output storage (optional) if (write.out) { writeEot(out, path.out = path.out, prefix = prefix, ...) df.name <- paste(prefix, "eot_locations.csv", sep = "_") - if (n == 1) { - write.table(location.df, col.names = TRUE, - paste(path.out, df.name, sep = "/"), - row.names = FALSE, append = FALSE, sep = ",") + if (n == 1) { + write.table( + location.df + , col.names = TRUE + , paste(path.out, df.name, sep = "/") + , row.names = FALSE + , append = FALSE + , sep = "," + ) } else { - write.table(location.df, col.names = FALSE, - paste(path.out, df.name, sep = "/"), - row.names = FALSE, append = TRUE, sep = ",") + write.table( + location.df + , col.names = FALSE + , paste(path.out, df.name, sep = "/") + , row.names = FALSE + , append = TRUE + , sep = "," + ) } rm(list = c("eot.ts", - "maxxy", - "location.df", - # "expl.var", - "rst.x.r", - "rst.x.rsq", - "rst.x.rsq.sums", - "rst.x.intercept", - "rst.x.slp", - "rst.x.p", - "brck.x.resids", - "rst.y.r", - "rst.y.rsq", - "rst.y.intercept", - "rst.y.slp", - "rst.y.p", - "brck.y.resids")) + "maxxy", + "location.df", + # "expl.var", + "rst.x.r", + "rst.x.rsq", + "rst.x.rsq.sums", + "rst.x.intercept", + "rst.x.slp", + "rst.x.p", + "brck.x.resids", + "rst.y.r", + "rst.y.rsq", + "rst.y.intercept", + "rst.y.slp", + "rst.y.p", + "brck.y.resids")) gc() } - - return(out) - -} + + return(out) + + } +) \ No newline at end of file diff --git a/R/anomalize.R b/R/anomalize.R index 2268956..7902f4b 100644 --- a/R/anomalize.R +++ b/R/anomalize.R @@ -1,41 +1,114 @@ -#' Create an anomaly RasterStack +methods::setGeneric( + "anomalize" + , function(x, ...) { + standardGeneric("anomalize") + } +) + +#' Create an anomaly raster series #' -#' @description -#' The function creates an anomaly RasterStack either based on the -#' overall mean of the original stack, or a supplied reference RasterLayer. -#' For the creation of seasonal anomalies use [deseason()]. +#' @description The function creates an anomaly raster series either based on +#' the overall mean of the original series, or a supplied reference raster. +#' For the creation of seasonal anomalies use [deseason()]. #' -#' @param x a RasterStack -#' @param reference an optional RasterLayer to be used as the reference -#' @param ... additional arguments passed to [raster::calc()] (and, in turn, -#' [raster::writeRaster()]) which is used under the hood +#' @param x A `SpatRaster` (or `Raster*`) series. +#' @param reference An optional single-layer `SpatRaster` (or `RasterLayer`) to +#' be used as the reference. Uses the overall mean of the original series if +#' `NULL` (default). +#' @param ... Additional arguments passed to [terra::app()] (e.g. 'cores', +#' 'filename') to calculate the overall mean if 'reference' is `NULL`, or to +#' the underlying `SpatRaster` method for `Raster*` input. #' -#' @return an anomaly RasterStack +#' @return An anomaly `SpatRaster` series. #' #' @seealso -#' [deseason()], [denoise()], [raster::calc()] +#' [deseason()], [denoise()] #' #' @export anomalize +#' @name anomalize #' #' @examples -#' data(australiaGPCP) -#' -#' aus_anom <- anomalize(australiaGPCP) +#' pcp = terra::unwrap(australiaGPCP) +#' pcp_anom = anomalize(pcp) #' -#' opar <- par(mfrow = c(1,2)) -#' plot(australiaGPCP[[10]], main = "original") -#' plot(aus_anom[[10]], main = "anomalized") +#' opar = par(mfrow = c(1,2)) +#' plot(pcp[[10]], main = "original") +#' plot(pcp_anom[[10]], main = "anomalized") #' par(opar) -anomalize <- function(x, - reference = NULL, - ...) { - - if (is.null(reference)) { - mn <- raster::calc(x, fun = mean, ...) - } else { - mn <- reference + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases anomalize,RasterStackBrick-method +#' @rdname anomalize +methods::setMethod( + "anomalize" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + anomalize( + terra::rast(x) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases anomalize,SpatRaster-method +#' @rdname anomalize +methods::setMethod( + "anomalize" + , signature(x = "SpatRaster") + , function( + x + , reference = NULL + , ... + ) { + + if (is.null(reference)) { + reference = terra::app( + x + , fun = mean + , ... + ) + } + + ## early exit: 'reference' is not a raster + if (!inherits(reference, what = c("SpatRaster", "Raster"))) { + stop( + sprintf( + "Expected 'reference' to inherit from [%s], but got [%s]." + , paste0(c("SpatRaster", "Raster"), collapse = ", ") + , class(reference)[1L] + ) + , call. = FALSE + ) + } + + ## if required, convert 'reference' to `SpatRaster` + if (inherits(reference, what = "Raster")) { + reference = terra::rast(reference) + } + + ## if required, use only the first 'reference' layer + if (terra::nlyr(reference) > 1L) { + warning( + sprintf( + paste( + "Expected 'reference' to have a single layer, but got [%s]." + , "Using the first layer only." + ) + , terra::nlyr(reference) + ) + , call. = FALSE + ) + reference = reference[[1L]] + } + + return(x - reference) } - - return(x - mn) - -} \ No newline at end of file +) \ No newline at end of file diff --git a/R/calcVar.R b/R/calcVar.R index 02b3f05..d9c12c9 100644 --- a/R/calcVar.R +++ b/R/calcVar.R @@ -1,33 +1,70 @@ -#' Calculate space-time variance of a RasterStack or RasterBrick +methods::setGeneric( + "calcVar" + , function(x, ...) { + standardGeneric("calcVar") + } +) + +#' Calculate space-time variance of a raster series #' -#' @description -#' The function calculates the (optionally standardised) space-time -#' variance of a RasterStack or RasterBrick. +#' @description The function calculates the (optionally standardised) space-time +#' variance of a raster series. #' -#' @param x a RasterStack or RasterBrick -#' @param standardised logical. -#' @param ... currently not used +#' @param x A `SpatRaster` (or `Raster*`) series. +#' @param standardised `logical`, defaults to `FALSE`. +#' @param ... For `SpatRaster` input: currently not used. For `Raster*` input: +#' arguments passed to the underlying `SpatRaster` method. #' -#' @return the mean (optionally standardised) space-time variance. +#' @return The mean (optionally standardised) space-time variance as `numeric`. #' #' @export calcVar +#' @name calcVar #' -#' @examples -#' data("pacificSST") +#' @examples +#' sst = terra::unwrap(pacificSST) #' -#' calcVar(pacificSST) -calcVar <- function(x, standardised = FALSE, ...) { - - if (!standardised) { - t <- mean(apply(raster::getValues(x), 1, var, na.rm = TRUE), - na.rm = TRUE) - s <- mean(apply(raster::getValues(x), 2, var, na.rm = TRUE), - na.rm = TRUE) - vrnc <- t + s - } else { - vrnc <- var(as.vector(raster::getValues(x)), na.rm = TRUE) +#' calcVar(sst) # default non-standardised +#' calcVar(sst, standardised = TRUE) + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases calcVar,RasterStackBrick-method +#' @rdname calcVar +methods::setMethod( + "calcVar" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + calcVar( + terra::rast(x) + , ... + ) } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases calcVar,SpatRaster-method +#' @rdname calcVar +methods::setMethod( + "calcVar" + , signature(x = "SpatRaster") + , function(x, standardised = FALSE, ...) { - return(vrnc) - -} \ No newline at end of file + if (!standardised) { + # compute variance across time and space, leveraging c++ functions for + # speed (see `?terra::app` for details) + t <- mean(terra::values(terra::app(x, "sd", na.rm = TRUE)^2)[, 1L]) + s <- mean((terra::global(x, fun = "sd", na.rm = TRUE)[, 1L])^2) + vrnc <- t + s + } else { + vrnc <- var(as.vector(terra::values(x)), na.rm = TRUE) + } + + return(vrnc) + } +) \ No newline at end of file diff --git a/R/covWeight.R b/R/covWeight.R index 15b4961..3a7945f 100644 --- a/R/covWeight.R +++ b/R/covWeight.R @@ -1,19 +1,22 @@ #' Create a weighted covariance matrix #' -#' @param m a matrix (e.g. as returned by [raster::getValues()]) -#' @param weights a numeric vector of weights. For lat/lon data this -#' can be produced with [getWeights()] -#' @param ... additional arguments passed to [stats::cov.wt()] +#' @param m A `matrix`, e.g. as returned by [terra::values()]. +#' @param weights A `numeric` vector of weights. For lat/lon data this can be +#' produced with [getWeights()]. +#' @param ... Additional arguments passed to [stats::cov.wt()] #' #' @return -#' see [stats::cov.wt()] +#' See [stats::cov.wt()]. #' -#' @seealso -#' [stats::cov.wt()] -#' -#' @export covWeight +#' @export covWeight <- function(m, weights, ...) { - cov.wt(na.exclude(m), weights, cor = TRUE, ...) + # TODO: `SpatRaster` method + stats::cov.wt( + stats::na.exclude(m) + , weights + , cor = TRUE + , ... + ) } \ No newline at end of file diff --git a/R/cutStack.R b/R/cutStack.R index c22fff4..d090765 100644 --- a/R/cutStack.R +++ b/R/cutStack.R @@ -1,42 +1,83 @@ -#' Shorten a RasterStack +methods::setGeneric( + "cutStack" + , function(x, ...) { + standardGeneric("cutStack") + } +) + +#' Shorten a raster series #' -#' @description -#' The function cuts a specified number of layers off a RrasterStack in -#' order to create lagged RasterStacks. +#' @description The function cuts a specified number of layers off a raster +#' series in order to create a lagged stack. #' -#' @param x a RasterStack -#' @param tail logical. If `TRUE` the layers will be taken off +#' @param x A `SpatRaster` (or `Raster*`) series. +#' @param tail `logical`. If `TRUE` (default) the layers will be taken off #' the end of the stack. If `FALSE` layers will be taken off #' the beginning. -#' @param n the number of layers to take away. +#' @param n The number of layers to take away as `integer`. If `NULL` (default), +#' 'x' is returned unchanged. +#' @param ... Additional arguments passed to the underlying `SpatRaster` method. #' -#' @return a RasterStack shortened by \code{n} layers either from the -#' beginning or the end, depending on the specification of \code{tail} +#' @return A `SpatRaster` series shortened by 'n' layers either from the +#' beginning or the end, depending on the specification of 'tail'. #' #' @examples -#' data(australiaGPCP) +#' pcp = terra::unwrap(australiaGPCP) #' #' # 6 layers from the beginning -#' cutStack(australiaGPCP, tail = FALSE, n = 6) +#' cutStack(pcp, tail = FALSE, n = 6) #' # 8 layers from the end -#' cutStack(australiaGPCP, tail = TRUE, n = 8) +#' cutStack(pcp, tail = TRUE, n = 8) #' #' @export cutStack -cutStack <- function(x, - tail = TRUE, - n = NULL) { +#' @name cutStack + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases cutStack,RasterStackBrick-method +#' @rdname cutStack +methods::setMethod( + "cutStack" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + cutStack( + terra::rast(x) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases cutStack,SpatRaster-method +#' @rdname cutStack +methods::setMethod( + "cutStack" + , signature(x = "SpatRaster") + , function( + x + , tail = TRUE + , n = NULL + ) { - # Return unmodified RasterStack if n == NULL - if (is.null(n)) { - return(x) - } else { - # Supplied RasterStack is predictor: - if (tail) { - return(x[[1:(raster::nlayers(x)-n)]]) - # Supplied RasterStack is response: + ## return unmodified raster series if `n == NULL` + if (is.null(n)) { + return(x) + } + + ## take away layers from the end, e.g. if supplied series is predictor + idx = if (tail) { + 1:(terra::nlyr(x) - n) + ## take away layers from the start, e.g. if supplied series is response } else { - return(x[[(n+1):raster::nlayers(x)]]) + seq.int(n + 1L, terra::nlyr(x)) } + + return(x[[idx]]) } - -} +) \ No newline at end of file diff --git a/R/deg2rad.R b/R/deg2rad.R index aafdf71..8db9fa0 100644 --- a/R/deg2rad.R +++ b/R/deg2rad.R @@ -1,20 +1,34 @@ #' Convert degrees to radians #' -#' @export deg2rad +#' @description Converts a series of degrees to radians. #' -#' @param deg vector of degrees to be converted to radians +#' @param deg A series of degrees to be converted to radians. #' -#' @examples -#' data(vdendool) +#' @return A vector of radians. #' +#' @examples #' ## latitude in degrees -#' degrees <- coordinates(vdendool)[, 2] +#' gph = terra::unwrap(vdendool) +#' +#' degrees = terra::crds(gph)[, 2] #' head(degrees) #' #' ## latitude in radians -#' radians <- deg2rad(coordinates(vdendool)[, 2]) +#' radians = deg2rad(degrees) #' head(radians) #' +#' ## `SpatRaster` input, e.g. useful for topographic operations +#' tmp = gph[[1L]] +#' terra::values(tmp) = degrees +#' +#' tmp_rad = deg2rad(tmp) +#' +#' opar = par(mfrow = c(1, 2)) +#' plot(tmp, main = "lat (degrees)") +#' plot(tmp_rad, main = "lat (radians)") +#' par(opar) +#' +#' @export deg2rad <- function(deg) { radians <- deg * pi / 180 diff --git a/R/denoise.R b/R/denoise.R index 6403581..94963ea 100644 --- a/R/denoise.R +++ b/R/denoise.R @@ -1,97 +1,158 @@ +methods::setGeneric( + "denoise" + , function(x, ...) { + standardGeneric("denoise") + } +) + #' Noise filtering through principal components #' -#' Filter noise from a RasterStack by decomposing into principal components -#' and subsequent reconstruction using only a subset of components +#' @description +#' Filter noise from a raster series by decomposing into principal components +#' and subsequent reconstruction using only a subset of components. #' -#' @param x RasterStack to be filtered -#' @param k number of components to be kept for reconstruction -#' (ignored if \code{expl.var} is supplied) -#' @param expl.var minimum amount of variance to be kept after reconstruction -#' (should be set to NULL or omitted if \code{k} is supplied) +#' @param x A `SpatRaster` (or `Raster*`) series to be filtered. +#' @param k The number of components to be kept for reconstruction (ignored if +#' 'expl.var' is supplied). +#' @param expl.var Minimum amount of variance to be kept after reconstruction +#' (should be set to `NULL` or omitted if 'k' is supplied). #' @param weighted logical. If `TRUE` the covariance matrix will be -#' geographically weighted using the cosine of latitude during decomposition -#' (only important for lat/lon data) -#' @param use.cpp logical. Determines whether to use \strong{Rcpp} -#' functionality, defaults to `TRUE`. -#' @param verbose logical. If `TRUE` some details about the -#' calculation process will be output to the console -#' @param ... additional arguments passed to [stats::princomp()] +#' geographically weighted using the cosine of latitude during decomposition +#' (only important for lat/lon data). +#' @param use.cpp logical. Determines whether to use **Rcpp** functionality, +#' defaults to `TRUE`. +#' @param verbose logical. If `TRUE` some details about the calculation process +#' will be output to the console. +#' @param ... Additional arguments passed to [stats::princomp()]. +#' +#' @return A denoised `SpatRaster` series. #' -#' @return a denoised RasterStack +#' @note +#' Either 'k' or 'expl.var' must be specified. If both are supplied, 'k' will be +#' ignored. If none are supplied, an error will be thrown. #' #' @seealso #' [anomalize()], [deseason()] #' #' @export denoise +#' @name denoise #' #' @examples -#' data("vdendool") -#' vdd_dns <- denoise(vdendool, expl.var = 0.8) +#' gph = terra::unwrap(vdendool) +#' gph_dns = denoise(gph, expl.var = 0.8) #' -#' opar <- par(mfrow = c(1,2)) -#' plot(vdendool[[1]], main = "original") -#' plot(vdd_dns[[1]], main = "denoised") +#' opar = par(mfrow = c(1,2)) +#' plot(gph[[1]], main = "original") +#' plot(gph_dns[[1]], main = "denoised") #' par(opar) -denoise <- function(x, - k = NULL, - expl.var = NULL, - weighted = TRUE, - use.cpp = TRUE, - verbose = TRUE, - ...) { - - x.vals <- raster::getValues(x) - #x.vals[is.na(x.vals)] <- 0 - - # PCA - if (weighted) { - pca <- princomp(~ x.vals, covmat = covWeight(x.vals, - remote::getWeights(x)), - scores = TRUE, na.action = na.exclude, ...) - } else { - pca <- princomp(~ x.vals, scores = TRUE, na.action = na.exclude, ...) - } - - # declare reconstruction characteristics according to supplied values - if (!is.null(expl.var)) { - k <- which(cumsum(pca$sdev^2 / sum(pca$sdev^2)) >= expl.var)[1] - } else { - expl.var <- cumsum(pca$sdev^2 / sum(pca$sdev^2))[k] - } - - if (verbose) { - cat("\n", - "Using the first ", - k, - " components (of ", - raster::nlayers(x), - ") to reconstruct series...\n", - " these account for ", - expl.var, - " of variance in orig. series\n\n", - sep = "") + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases denoise,RasterStackBrick-method +#' @rdname denoise +methods::setMethod( + "denoise" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + denoise( + terra::rast(x) + , ... + ) } - - # Reconstruction - recons <- lapply(seq(nlayers(x)), function(i) { - rowSums(t(as.matrix(pca$loadings[, 1:k])[i, ] * - t(pca$scores[, 1:k]))) + pca$center[i] - }) - - # Insert reconstructed values in original data set - if (use.cpp) { - mat <- raster::as.matrix(x) - jnk <- insertReconsC(recons, mat) - rst <- raster::setValues(x, jnk) - } else { - rst <- raster::brick(lapply(seq(recons), function(i) { - tmp_rst <- x[[i]] - tmp_rst[] <- recons[[i]] - return(tmp_rst) - })) +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases denoise,SpatRaster-method +#' @rdname denoise +methods::setMethod( + "denoise" + , signature(x = "SpatRaster") + , function( + x + , k = NULL + , expl.var = NULL + , weighted = TRUE + , use.cpp = TRUE + , verbose = TRUE + , ... + ) { + + x.vals <- terra::values(x) + #x.vals[is.na(x.vals)] <- 0 + + # PCA + if (weighted) { + pca <- stats::princomp( + ~ x.vals + , covmat = covWeight( + x.vals + , getWeights(x) + ) + , scores = TRUE + , na.action = stats::na.exclude + , ... + ) + } else { + pca <- stats::princomp( + ~ x.vals + , scores = TRUE + , na.action = stats::na.exclude + , ... + ) + } + + # declare reconstruction characteristics according to supplied values + stopifnot( + "Either 'expl.var' or 'k' must be supplied." = + !is.null(expl.var) | !is.null(k) + ) + + if (!is.null(expl.var)) { + k <- which(cumsum(pca$sdev^2 / sum(pca$sdev^2)) >= expl.var)[1] + } else { + expl.var <- cumsum(pca$sdev^2 / sum(pca$sdev^2))[k] + } + + if (verbose) { + paste( + "\nUsing the first %s components (of %s) to reconstruct series..." + , "these account for %s of variance in orig. series\n\n" + , sep = "\n " + ) |> + sprintf( + k + , terra::nlyr(x) + , expl.var + ) |> + cat() + } + + # Reconstruction + recons <- lapply(seq(terra::nlyr(x)), function(i) { + rowSums(t(as.matrix(pca$loadings[, 1:k])[i, ] * + t(pca$scores[, 1:k]))) + pca$center[i] + }) + + # Insert reconstructed values in original data set + # TODO: discard 'use.cpp' option and always use {terra} for speed-up + if (use.cpp) { + jnk <- insertReconsC(recons, x.vals) + rst <- terra::setValues(x, jnk) + } else { + rst = terra::setValues( + x + , values = do.call(cbind, recons) + ) + } + + # Return denoised data set + return(rst) + } - - # Return denoised data set - return(rst) - -} +) \ No newline at end of file diff --git a/R/deseason.R b/R/deseason.R index c48917e..4ce8bce 100644 --- a/R/deseason.R +++ b/R/deseason.R @@ -1,27 +1,30 @@ -if ( !isGeneric("deseason") ) { - setGeneric("deseason", function(x, ...) - standardGeneric("deseason")) -} +methods::setGeneric( + "deseason" + , function(x, ...) { + standardGeneric("deseason") + } +) + #' Create seasonal anomalies #' -#' @description -#' The function calculates anomalies of a RasterStack by supplying a -#' suitable seasonal window. E. g. to create monthly anomalies of a -#' raster stack of 12 layers per year, use \code{cycle.window = 12}. +#' @description The function calculates anomalies of a raster object by +#' supplying a suitable seasonal window. E.g. to create monthly anomalies of a +#' raster stack of 12 layers per year, use `cycle.window = 12`. #' -#' @param x An `Raster*` object or, alternatively, a `numeric` time -#' series. -#' @param cycle.window `integer`, defaults to \code{12}. The window for the -#' creation of the anomalies. -#' @param use.cpp `logical`, defaults to `FALSE`. Determines whether -#' or not to use \strong{Rcpp} functionality. Only applies if \code{x} is a -#' `Raster*` object. +#' @param x A `SpatRaster` (or `Raster*`) object or, alternatively, a `numeric` +#' time series. +#' @param cycle.window `integer`, defaults to `12`. The window for the creation +#' of the anomalies. +#' @param use.cpp `logical`, defaults to `FALSE`. Determines whether or not to +#' use **Rcpp** functionality. Only applies if `x` is a raster object. #' @param filename `character`. Output filename (optional). -#' @param ... Additional arguments passed on to [raster::writeRaster()], only -#' considered if \code{filename} is specified. +#' @param ... For `SpatRaster` input: additional arguments passed to +#' [terra::writeRaster()], only considered if `filename` is specified. +#' For `Raster*` input: arguments passed to the underlying `SpatRaster` +#' method (e.g. `cycle.window`, `use.cpp`, `filename`). #' -#' @return If \code{x} is a `Raster*` object, a deseasoned -#' \code{RasterStack}; else a deseasoned `numeric` vector. +#' @return If `x` is a raster object, a deseasoned `SpatRaster`; else a +#' deseasoned `numeric` vector. #' #' @seealso #' [anomalize()], [denoise()] @@ -30,13 +33,12 @@ if ( !isGeneric("deseason") ) { #' @name deseason #' #' @examples -#' data("australiaGPCP") -#' -#' aus_dsn <- deseason(australiaGPCP, 12) +#' pcp = terra::unwrap(australiaGPCP) +#' pcp_dsn = deseason(pcp, cycle.window = 12) #' -#' opar <- par(mfrow = c(1,2)) -#' plot(australiaGPCP[[1]], main = "original") -#' plot(aus_dsn[[1]], main = "deseasoned") +#' opar = par(mfrow = c(1, 2)) +#' plot(pcp[[1]], main = "original") +#' plot(pcp_dsn[[1]], main = "deseasoned") #' par(opar) @@ -44,8 +46,27 @@ if ( !isGeneric("deseason") ) { ### function using 'RasterStack' or 'RasterBrick' ############################## #' @aliases deseason,RasterStackBrick-method #' @rdname deseason +setMethod( + "deseason" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + deseason( + terra::rast(x) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases deseason,SpatRaster-method +#' @rdname deseason setMethod("deseason", - signature(x = "RasterStackBrick"), + signature(x = "SpatRaster"), function(x, cycle.window = 12L, use.cpp = FALSE, @@ -54,18 +75,18 @@ setMethod("deseason", if (use.cpp) { ## raster to matrix - mat <- as.matrix(x) + mat <- terra::as.matrix(x) ## deseasoning mat_mv <- monthlyMeansC(mat, cycle.window) x_mv <- x[[1:cycle.window]] - x_mv <- setValues(x_mv, values = mat_mv) + x_mv <- terra::setValues(x_mv, values = mat_mv) } else { # Calculate layer averages based on supplied seasonal window - x_mv <- raster::stack(rep(lapply(1:cycle.window, function(i) { - raster::calc(x[[seq(i, raster::nlayers(x), cycle.window)]], + x_mv <- terra::rast(lapply(1:cycle.window, function(i) { + terra::app(x[[seq(i, terra::nlyr(x), cycle.window)]], fun = mean, na.rm = TRUE) - }), raster::nlayers(x) / cycle.window)) + })) } # Subtract monthly averages from actually measured values @@ -73,7 +94,7 @@ setMethod("deseason", # Write to file (optional) if (filename != "") - x_dsn <- raster::writeRaster(x_dsn, filename = filename, ...) + x_dsn <- terra::writeRaster(x_dsn, filename = filename, ...) # Return output return(x_dsn) diff --git a/R/eot.R b/R/eot.R index 766dd48..9be6e07 100644 --- a/R/eot.R +++ b/R/eot.R @@ -1,40 +1,42 @@ -if ( !isGeneric('eot') ) { - setGeneric('eot', function(x, ...) - standardGeneric('eot')) -} +methods::setGeneric( + "eot" + , function(x, ...) { + standardGeneric("eot") + } +) -#' EOT analysis of a predictor and (optionally) a response RasterStack +#' EOT analysis of a predictor and (optionally) a response raster series #' #' @description -#' Calculate a given number of EOT modes either internally or between -#' RasterStacks. +#' Calculate a given number of EOT modes either internally or between raster +#' series. #' -#' @param x a `Raster*` object used as predictor -#' @param y a `Raster*` object used as response. If \code{y} is -#' `NULL`, \code{x} is used as \code{y} -#' @param n the number of EOT modes to calculate -#' @param standardised logical. If `FALSE` the calculated r-squared values -#' will be multiplied by the variance -#' @param write.out logical. If `TRUE` results will be written to disk -#' using \code{path.out} -#' @param path.out the file path for writing results if \code{write.out} is `TRUE`. -#' Defaults to current working directory -#' @param prefix optional prefix to be used for naming of results if -#' \code{write.out} is `TRUE` -#' @param reduce.both logical. If `TRUE` both \code{x} and \code{y} -#' are reduced after each iteration. If `FALSE` only \code{y} is reduced -#' @param type the type of the link function. Defaults to \code{'rsq'} as in original -#' proposed method from \cite{van den Dool 2000}. If set to \code{'ioa'} index of agreement is -#' used instead -#' @param verbose logical. If `TRUE` some details about the -#' calculation process will be output to the console -#' @param ... not used at the moment +#' @param x A `SpatRaster` (or `Raster*`) object used as predictor. +#' @param y A `SpatRaster` (or `Raster*`) object used as response. If 'y' is +#' `NULL` (default), 'x' is used as response. +#' @param n The number of EOT modes to calculate as `integer`, defaults to `1`. +#' @param standardised `logical`, default `TRUE`. If `FALSE` the calculated +#' r-squared values will be multiplied by the variance. +#' @param write.out `logical`, default `FALSE`. If `TRUE` results will be +#' written to disk using 'path.out'. +#' @param path.out The file path for writing results if 'write.out' is `TRUE`. +#' Defaults to current working directory. +#' @param prefix optional prefix to be used for naming of results if 'write.out' +#' is `TRUE`. +#' @param reduce.both `logical`, default `FALSE`. If `TRUE` both 'x' and 'y' are +#' reduced after each iteration. If `FALSE` only 'y' is reduced. +#' @param type The type of the link function. Defaults to `"rsq"` as in original +#' proposed method from \cite{van den Dool 2000}. If set to `"ioa"` index of +#' agreement is used instead. +#' @param verbose `logical`, default `FALSE`. If `TRUE` some details about the +#' calculation process will be output to the console. +#' @param ... Currently not used. #' #' @details #' For a detailed description of the EOT algorithm and the mathematics behind it, #' see the References section. In brief, the algorithm works as follows: -#' First, the temporal profiles of each pixel \emph{xp} of the predictor domain -#' are regressed against the profiles of all pixels \emph{xr} in the +#' First, the temporal profiles of each pixel _xp_ of the predictor domain +#' are regressed against the profiles of all pixels _xr_ in the #' response domain. #' The calculated coefficients of determination are summed up and the pixel #' with the highest sum is identified as the 'base point' of the first/leading mode. @@ -42,32 +44,33 @@ if ( !isGeneric('eot') ) { #' Then, the residuals from the regression are taken to be the basis #' for the calculation of the next EOT, thus ensuring orthogonality #' of the identified teleconnections. This procedure is repeated until -#' a predefined amount of \emph{n} EOTs is calculated. In general, +#' a predefined amount of _n_ EOTs is calculated. In general, #' \pkg{remote} implements a 'brute force' spatial data mining approach to #' identify locations of enhanced potential to explain spatio-temporal #' variability within the same or another geographic field. #' #' @return -#' if n = 1 an \emph{EotMode}, if n > 1 an \emph{EotStack} of \code{n} -#' \emph{EotMode}s. Each \emph{EotMode} has the following components: +#' If `n = 1`` an `EotMode`, if n > 1 an `EotStack` of 'n' `EotMode`s. Each +#' `EotMode` has the following components: +#' #' \itemize{ #' \item \emph{mode} - the number of the identified mode (1 - n) #' \item \emph{eot} - the EOT (time series) at the identified base point. -#' Note, this is a simple numeric vector, not of class \code{ts} +#' Note, this is a simple numeric vector, not of class `ts` #' \item \emph{coords_bp} - the coordinates of the identified base point #' \item \emph{cell_bp} - the cell number of the indeified base point #' \item \emph{cum_exp_var} - the (cumulative) explained variance of the considered EOT -#' \item \emph{r_predictor} - the \emph{RasterLayer} of the correlation coefficients +#' \item \emph{r_predictor} - the `SpatRaster` of the correlation coefficients #' between the base point and each pixel of the predictor domain #' \item \emph{rsq_predictor} - as above but for the coefficient of determination #' \item \emph{rsq_sums_predictor} - as above but for the sums of coefficient of determination -#' \item \emph{int_predictor} - the \emph{RasterLayer} of the intercept of the +#' \item \emph{int_predictor} - the `SpatRaster` of the intercept of the #' regression equation for each pixel of the predictor domain #' \item \emph{slp_predictor} - same as above but for the slope of the #' regression equation for each pixel of the predictor domain -#' \item \emph{p_predictor} - the \emph{RasterLayer} of the significance (p-value) +#' \item \emph{p_predictor} - the `SpatRaster` of the significance (p-value) #' of the the regression equation for each pixel of the predictor domain -#' \item \emph{resid_predictor} - the \emph{RasterBrick} of the reduced data +#' \item \emph{resid_predictor} - the `SpatRaster` of the reduced data #' for the predictor domain #' } #' @@ -75,8 +78,7 @@ if ( !isGeneric('eot') ) { #' also returned for the \emph{*_response} domain, #' even if predictor and response domain are equal. This is due to that fact, #' that if not both fields are reduced after the first EOT is found, -#' these \emph{RasterLayers} will differ. -#' +#' these `SpatRaster`s will differ. #' #' @references #' \bold{Empirical Orthogonal Teleconnections}\cr @@ -93,214 +95,139 @@ if ( !isGeneric('eot') ) { #' ### EXAMPLE I #' ### a single field #' \donttest{ -#' data(vdendool) +#' gph = terra::unwrap(vdendool) #' -#' ## claculate 2 leading modes -#' nh_modes <- eot(x = vdendool, y = NULL, n = 2, +#' ## calculate 2 leading modes +#' nh_modes <- eot(x = gph, y = NULL, n = 2, #' standardised = FALSE, #' verbose = TRUE) #' #' plot(nh_modes, y = 1, show.bp = TRUE) #' plot(nh_modes, y = 2, show.bp = TRUE) #' } +#' #' @export #' @name eot -#' @rdname eot -#' @aliases eot,RasterStackBrick-method -# set methods ------------------------------------------------------------- -setMethod('eot', signature(x = 'RasterStackBrick'), - function(x, - y = NULL, - n = 1, - standardised = TRUE, - write.out = FALSE, - path.out = ".", - prefix = "remote", - reduce.both = FALSE, - type = c("rsq", "ioa"), - verbose = TRUE, - ...) { - - # Duplicate predictor set in case predictor and response are identical - if (is.null(y)) { - y <- x - } - - orig.var <- calcVar(y, standardised = standardised) - - ### EOT - - # Loop through number of desired EOTs - for (z in 1:n) { - - # Use initial response data set in case of first iteration - if (z == 1) { - - x.eot <- EotCycle(x = x, - y = y, - n = z, - type = type, - standardised = standardised, - orig.var = orig.var, - write.out = write.out, - path.out = path.out, - verbose = verbose, - prefix = prefix) - - names(x.eot) <- paste("mode_", sprintf("%02.f", z), - sep = "") - - # Use last entry of slot 'residuals' otherwise - } else if (z > 1) { - tmp.x.eot <- EotCycle( - x = if (!reduce.both) { - x - } else { - if (z == 2) { - x.eot@resid_predictor - } else { - x.eot[[z-1]]@resid_predictor - } - }, - y = if (z == 2) { - x.eot@resid_response - } else { - x.eot[[z-1]]@resid_response - }, - # y.eq.x = y.eq.x, - n = z, - type = type, - standardised = standardised, - orig.var = orig.var, - write.out = write.out, - path.out = path.out, - verbose = verbose, - prefix = prefix) - - if (z == 2) { - x.eot <- list(x.eot, tmp.x.eot) - names(x.eot) <- c(paste("mode_", sprintf("%02.f", 1), - sep = ""), - paste("mode", sprintf("%02.f", z), - sep = "_")) - } else { - tmp.names <- names(x.eot) - x.eot <- append(x.eot, list(tmp.x.eot)) - names(x.eot) <- c(tmp.names, - paste("mode", sprintf("%02.f", z), - sep = "_")) - } - } - } - - if (length(x.eot) == 1) { - out <- x.eot - } else { - out <- new('EotStack', modes = x.eot, names = names(x.eot)) - } - return(out) - } +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases eot,RasterStackBrick-method +#' @rdname eot +methods::setMethod( + "eot" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + eot( + terra::rast(x) + , ... + ) + } ) -# #' @describeIn eot -# setMethod('eot', signature(x = 'RasterBrick'), -# function(x, -# y = NULL, -# n = 1, -# standardised = TRUE, -# write.out = FALSE, -# path.out = ".", -# prefix = "remote", -# reduce.both = FALSE, -# type = c("rsq", "ioa"), -# verbose = TRUE, -# ...) { -# -# # Duplicate predictor set in case predictor and response are identical -# if (is.null(y)) { -# y <- x -# y.eq.x <- TRUE -# } else { -# y.eq.x <- FALSE -# } -# -# orig.var <- calcVar(y, standardised = standardised) -# -# ### EOT -# -# # Loop through number of desired EOTs -# for (z in seq(n)) { -# -# # Use initial response data set in case of first iteration -# if (z == 1) { -# -# x.eot <- EotCycle(x = x, -# y = y, -# y.eq.x = y.eq.x, -# n = z, -# type = type, -# standardised = standardised, -# orig.var = orig.var, -# write.out = write.out, -# path.out = path.out, -# verbose = verbose, -# prefix = prefix) -# -# names(x.eot) <- paste("mode_", sprintf("%02.f", z), -# sep = "") -# -# # Use last entry of slot 'residuals' otherwise -# } else if (z > 1) { -# tmp.x.eot <- EotCycle( -# x = if (!reduce.both) { -# x -# } else { -# if (z == 2) { -# x.eot@resid_predictor -# } else { -# x.eot[[z-1]]@resid_predictor -# } -# }, -# y = if (z == 2) { -# x.eot@resid_response -# } else { -# x.eot[[z-1]]@resid_response -# }, -# y.eq.x = y.eq.x, -# n = z, -# type = type, -# standardised = standardised, -# orig.var = orig.var, -# write.out = write.out, -# path.out = path.out, -# verbose = verbose, -# prefix = prefix) -# -# if (z == 2) { -# x.eot <- list(x.eot, tmp.x.eot) -# names(x.eot) <- c(paste("mode_", sprintf("%02.f", 1), -# sep = ""), -# paste("mode", sprintf("%02.f", z), -# sep = "_")) -# } else { -# tmp.names <- names(x.eot) -# x.eot <- append(x.eot, list(tmp.x.eot)) -# names(x.eot) <- c(tmp.names, -# paste("mode", sprintf("%02.f", z), -# sep = "_")) -# } -# } -# } -# -# if (length(x.eot) == 1) { -# out <- x.eot -# } else { -# out <- new('EotStack', modes = x.eot, names = names(x.eot)) -# } -# return(out) -# } -# ) -# +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases eot,SpatRaster-method +#' @rdname eot +methods::setMethod( + "eot" + , signature(x = "SpatRaster") + , function( + x + , y = NULL + , n = 1 + , standardised = TRUE + , write.out = FALSE + , path.out = "." + , prefix = "remote" + , reduce.both = FALSE + , type = c("rsq", "ioa") + , verbose = TRUE + , ... # TODO: pass to `writeEot()` (via `EotCycle()`), e.g. 'filetype' + ) { + + type = match.arg(type) + + ## duplicate predictor set in case predictor and response are identical + if (is.null(y)) { + y <- x + } + + orig.var <- calcVar(y, standardised = standardised) + + ## loop through number of desired eots + for (z in 1:n) { + + # use initial response data set in case of first iteration + if (z == 1) { + + x.eot <- EotCycle( + x = x + , y = y + , n = z + , type = type + , standardised = standardised + , orig.var = orig.var + , write.out = write.out + , path.out = path.out + , verbose = verbose + , prefix = prefix + ) + + names(x.eot) <- sprintf("mode_%02d", z) + next + + } + + # use last entry of slot 'residuals' otherwise + tmp.x.eot <- EotCycle( + x = if (!reduce.both) { + x + } else { + if (z == 2) { + x.eot@resid_predictor + } else { + x.eot[[z-1]]@resid_predictor + } + }, + y = if (z == 2) { + x.eot@resid_response + } else { + x.eot[[z-1]]@resid_response + }, + # y.eq.x = y.eq.x, + n = z, + type = type, + standardised = standardised, + orig.var = orig.var, + write.out = write.out, + path.out = path.out, + verbose = verbose, + prefix = prefix + ) + + if (z == 2) { + x.eot <- list(x.eot, tmp.x.eot) + names(x.eot) <- sprintf("mode_%02d", c(1, z)) + } else { + tmp.names <- names(x.eot) + x.eot <- append(x.eot, list(tmp.x.eot)) + names(x.eot) <- c( + tmp.names + , sprintf("mode_%02d", z) + ) + } + } + + if (length(x.eot) == 1) { + out <- x.eot + } else { + out <- new('EotStack', modes = x.eot, names = names(x.eot)) + } + return(out) + } +) diff --git a/R/geoWeight.R b/R/geoWeight.R index 57dd47d..5272a52 100644 --- a/R/geoWeight.R +++ b/R/geoWeight.R @@ -1,37 +1,76 @@ +methods::setGeneric( + "geoWeight" + , function(x, ...) { + standardGeneric("geoWeight") + } +) + #' Geographic weighting #' -#' @description -#' The function performs geographic weighting of non-projected long/lat -#' data. By default it uses the cosine of latitude to compensate for the -#' area distortion, though the user can supply other functions via \code{f}. +#' @description +#' The function performs geographic weighting of non-projected long/lat data. By +#' default it uses the cosine of latitude (in radians) to compensate for the +#' area distortion, though the user can supply other weighting functions. #' +#' @param x,f,... See [getWeights()]. In case of `Raster*` input, '...' funnels +#' additional arguments (i.e. 'f') to the underlying `SpatRaster` method. #' -#' @param x a Raster* object -#' @param f a function to be used to the weighting. -#' Defaults to \code{cos(x)} -#' @param ... additional arguments to be passed to f -#' -#' @return a weighted Raster* object +#' @return A weighted `SpatRaster`. #' #' @export geoWeight +#' @name geoWeight #' #' @examples -#' data(vdendool) -#' -#' wgtd <- geoWeight(vdendool) +#' gph <- terra::unwrap(vdendool) +#' wgtd <- geoWeight(gph) #' #' opar <- par(mfrow = c(1,2)) -#' plot(vdendool[[1]], main = "original") +#' plot(gph[[1]], main = "original") #' plot(wgtd[[1]], main = "weighted") #' par(opar) -geoWeight <- function(x, - f = function(x) cos(x), - ...) { - - x.vals <- x[] - rads <- deg2rad(sp::coordinates(x)[, 2]) - x.weightd <- x.vals * f(rads, ...) - x[] <- x.weightd - return(x) - -} \ No newline at end of file + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases geoWeight,RasterStackBrick-method +#' @rdname geoWeight +methods::setMethod( + "geoWeight" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + geoWeight( + terra::rast(x) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases geoWeight,SpatRaster-method +#' @rdname geoWeight +methods::setMethod( + "geoWeight" + , signature(x = "SpatRaster") + , function( + x + , f = cos + , ... + ) { + + # TODO: test for epsg:4326 + # NOTE: can't use `getWeights()` hereafter because of `na.rm = TRUE` + # (required for compatibility with `princomp()`) therein + + x.vals <- x[] + rads <- deg2rad(terra::crds(x)[, 2L]) + x.weightd <- x.vals * f(rads, ...) + x[] <- x.weightd + return(x) + + } +) \ No newline at end of file diff --git a/R/getWeights.R b/R/getWeights.R index f4964b7..ef5fd35 100644 --- a/R/getWeights.R +++ b/R/getWeights.R @@ -1,31 +1,83 @@ +methods::setGeneric( + "getWeights" + , function(x, ...) { + standardGeneric("getWeights") + } +) + #' Calculate weights from latitude #' -#' Calculate weights using the cosine of latitude to compensate for area -#' distortion of non-projected lat/lon data +#' @description Calculate weights using the cosine of latitude to compensate for +#' area distortion of non-projected lat/lon data. #' -#' @param x a Raster* object -#' @param f a function to be used to the weighting. -#' Defaults to \code{cos(x)} -#' @param ... additional arguments to be passed to f +#' @param x A non-projected `SpatRaster` (or `Raster*`) object. +#' @param f A `function` applied to the latitude (in radians) to compute +#' weights. Defaults to `cos`. +#' @param ... Additional arguments passed to 'f', or to the underlying +#' `SpatRaster` method in general for `Raster*` input. #' -#' @return a numeric vector of weights +#' @return A `numeric` vector of weights for non-`NA` cells in 'x'. #' #' @examples -#' data("australiaGPCP") -#' wghts <- getWeights(australiaGPCP) -#' wghts_rst <- australiaGPCP[[1]] -#' wghts_rst[] <- wghts +#' pcp = terra::unwrap(australiaGPCP) +#' +#' wghts = getWeights(pcp) +#' utils::head(wghts) +#' +#' wghts_rst = terra::setValues(pcp[[1]], wghts) #' -#' opar <- par(mfrow = c(1,2)) -#' plot(australiaGPCP[[1]], main = "data") +#' opar = par(mfrow = c(1,2)) +#' plot(pcp[[1]], main = "data") #' plot(wghts_rst, main = "weights") #' par(opar) #' #' @export getWeights -getWeights <- function(x, - f = function(x) cos(x), - ...) { - - f(deg2rad(sp::coordinates(x)[, 2][!is.na(x[[1]][])]), ...) +#' @name getWeights + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases getWeights,RasterStackBrick-method +#' @rdname getWeights +methods::setMethod( + "getWeights" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + getWeights( + terra::rast(x) + , ... + ) + } +) + -} \ No newline at end of file +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases getWeights,SpatRaster-method +#' @rdname getWeights +methods::setMethod( + "getWeights" + , signature(x = "SpatRaster") + , function( + x + , f = cos + , ... + ) { + + # TODO: + # * what happens in the presence of `NA` values (length of weights vector is + # not a multiple of `ncell(x)`)? Use `na.all = TRUE` to account for all + # layers in 'x', not just first? + # * test for epsg:4326 + f( + deg2rad( + terra::crds(x, na.rm = TRUE)[, 2L] + ) + , ... + ) + + } +) \ No newline at end of file diff --git a/R/lagalize.R b/R/lagalize.R index dd475bb..73c9f01 100644 --- a/R/lagalize.R +++ b/R/lagalize.R @@ -1,51 +1,94 @@ -#' Create lagged RasterStacks +methods::setGeneric( + "lagalize" + , function(x, y, ...) { + standardGeneric("lagalize") + } +) + +#' Create lagged raster series #' #' @description -#' The function is used to produce two lagged RasterStacks. The second is cut +#' The function is used to produce two lagged raster series. The second is cut #' from the beginning, the first from the tail to ensure equal output lengths #' (provided that input lengths were equal). #' -#' @param x a RasterStack (to be cut from tail) -#' @param y a RasterStack (to be cut from beginning) -#' @param lag the desired lag (in the native frequency of the RasterStack) -#' @param freq the frequency of the RasterStacks -#' @param ... currently not used +#' @param x A `SpatRaster` (or `Raster*`) series to be cut from tail. +#' @param y A `SpatRaster` (or `Raster*`) series to be cut from beginning. +#' @param lag The desired lag in the native frequency of the series passed to +#' [cutStack()]. +#' @param freq The frequency of the raster series as `integer`. +#' @param ... For `SpatRaster` input: currently not used. For `Raster*` input: +#' arguments passed to the underlying `SpatRaster` method. #' #' @return -#' a list with the two RasterStacks lagged by \code{lag} +#' A `list` with the two raster series lagged by 'lag'. #' #' @examples -#' data(pacificSST) -#' data(australiaGPCP) +#' sst = terra::unwrap(pacificSST) +#' pcp = terra::unwrap(australiaGPCP) #' #' # lag GPCP by 4 months -#' lagged <- lagalize(pacificSST, australiaGPCP, lag = 4, freq = 12) -#' lagged[[1]][[1]] #check names to see date of layer -#' lagged[[2]][[1]] #check names to see date of layer +#' lagged = lagalize(sst, pcp, lag = 4, freq = 12) +#' lagged[[1]][[1]] # check names to see date of layer +#' lagged[[2]][[1]] # -"- #' -#' @export lagalize -lagalize <- function(x, - y, - lag = NULL, - freq = 12, - ...) { - - # Return list of unmodified RasterStacks if lag == NULL - if (is.null(lag)) { - return(list(x, y)) - } else { +#' @export +#' @name lagalize + + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases lagalize,RasterStackBrick,RasterStackBrick-method +#' @rdname lagalize +methods::setMethod( + "lagalize" + , signature(x = "RasterStackBrick", y = "RasterStackBrick") + , function( + x + , y + , ... + ) { + lagalize( + x = terra::rast(x) + , y = terra::rast(y) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases lagalize,SpatRaster,SpatRaster-method +#' @rdname lagalize +methods::setMethod( + "lagalize" + , signature(x = "SpatRaster", y = "SpatRaster") + , function( + x + , y + , lag = NULL + , freq = 12L + , ... + ) { + + # Return list of unmodified RasterStacks if lag == NULL + if (is.null(lag)) { + return(list(x, y)) + } + rest <- freq - lag # Lagalize predictor stack x.lag <- cutStack(x = x, tail = TRUE, n = lag) - x.lag.adj <- x.lag[[1:(raster::nlayers(x.lag) - rest)]] + x.lag.adj <- x.lag[[1:(terra::nlyr(x.lag) - rest)]] # Lagalize response stack y.lag <- cutStack(x = y, tail = FALSE, n = lag) - y.lag.adj <- y.lag[[1:(raster::nlayers(y.lag) - rest)]] + y.lag.adj <- y.lag[[1:(terra::nlyr(y.lag) - rest)]] # Return list of lagalized stacks return(list(x.lag, y.lag)) + } - -} \ No newline at end of file +) \ No newline at end of file diff --git a/R/longtermMeans.R b/R/longtermMeans.R index 62b5d26..acd2cb1 100644 --- a/R/longtermMeans.R +++ b/R/longtermMeans.R @@ -1,17 +1,26 @@ -#' Calculate long-term means from a 'RasterStack' +methods::setGeneric( + "longtermMeans" + , function(x, ...) { + standardGeneric("longtermMeans") + } +) + +#' Calculate long-term means from a raster series #' #' @description -#' Calculate long-term means from an input 'RasterStack' (or 'RasterBrick') -#' object. Ideally, the number of input layers should be divisable by the -#' supplied \code{cycle.window}. For instance, if \code{x} consists of monthly -#' layers, \code{cycle.window} should be a multiple of 12. +#' Calculate long-term means from an input raster series. Ideally, the number of +#' input layers should be divisable by the supplied 'cycle.window'. For +#' instance, if 'x' consists of monthly layers, 'cycle.window' should be a +#' multiple of `12`. #' -#' @param x A 'RasterStack' (or 'RasterBrick') object. -#' @param cycle.window 'integer'. See [deseason()]. +#' @param x A `SpatRaster` (or `Raster*`) series. +#' @param cycle.window `integer`, defaults to `12`. See [deseason()]. +#' @param ... For `Raster*` input, arguments passed to the underlying +#' `SpatRaster` method. #' #' @return -#' If \code{cycle.window} equals \code{nlayers(x)} (which obviously doesn't make -#' much sense), a 'RasterLayer' object; else a 'RasterStack' object. +#' A `SpatRaster` with 'cycle.window' layers, each containing the mean across +#' all corresponding time steps. #' #' @author #' Florian Detsch @@ -20,24 +29,55 @@ #' [deseason()]. #' #' @examples -#' data("australiaGPCP") +#' pcp = terra::unwrap(australiaGPCP) #' -#' longtermMeans(australiaGPCP) +#' longtermMeans(pcp) #' -#' @export longtermMeans +#' @export #' @name longtermMeans -longtermMeans <- function(x, cycle.window = 12L) { - ## raster to matrix - mat <- raster::as.matrix(x) - - ## long-term means - mat_ltm <- monthlyMeansC(mat, 12) - - ## insert values - rst_ltm <- x[[1:(raster::nlayers(x) / cycle.window)]] - rst_ltm <- raster::setValues(rst_ltm, NA) - - rst_ltm <- raster::setValues(rst_ltm, mat_ltm) - return(rst_ltm) -} + +################################################################################ +### function using 'RasterStackBrick' ########################################## +#' @aliases longtermMeans,RasterStackBrick-method +#' @rdname longtermMeans +methods::setMethod( + "longtermMeans" + , signature(x = "RasterStackBrick") + , function( + x + , ... + ) { + longtermMeans( + x = terra::rast(x) + , ... + ) + } +) + + +################################################################################ +### function using 'SpatRaster' ################################################ +#' @aliases longtermMeans,SpatRaster-method +#' @rdname longtermMeans +methods::setMethod( + "longtermMeans" + , signature(x = "SpatRaster") + , function(x, cycle.window = 12L) { + + ## insert values + idx = rep( + 1:cycle.window + , times = terra::nlyr(x) / cycle.window + ) + + rst_ltm = terra::tapp( + x + , index = idx + , fun = mean + , na.rm = TRUE + ) + + return(rst_ltm) + } +) \ No newline at end of file diff --git a/R/nXplain.R b/R/nXplain.R index a756dea..24ca73f 100644 --- a/R/nXplain.R +++ b/R/nXplain.R @@ -1,29 +1,31 @@ -if (!isGeneric('nXplain')) { - setGeneric('nXplain', function(x, ...) - standardGeneric('nXplain')) -} - +methods::setGeneric( + 'nXplain' + , function(x, ...) { + standardGeneric('nXplain') + } +) #' Number of EOTs needed for variance explanation #' #' @description -#' The function identifies the number of modes needed to explain a certain amount of -#' variance within the response field. +#' The function identifies the number of modes needed to explain a certain +#' amount of variance within the response field. #' -#' @param x an \emph{EotStack} -#' @param var the minimum amount of variance to be explained by the modes +#' @param x An `EotStack`. +#' @param var The minimum amount of variance to be explained by the modes as +#' `numeric`, defaults to `0.9`. #' -#' @note This is a post-hoc function. It needs an \emph{EotStack} -#' created as returned by [eot()]. Depending on the potency -#' of the identified EOTs, it may be necessary to compute a high number of -#' modes in order to be able to explain a large enough part of the variance. +#' @note This is a post-hoc function. It needs an `EotStack` created as returned +#' by [eot()]. Depending on the potency of the identified EOTs, it may be +#' necessary to compute a high number of modes in order to be able to explain +#' a large enough part of the variance. #' -#' @return an integer denoting the number of EOTs needed to explain \code{var} +#' @return The number of EOTs needed to explain 'var'. #' #' @examples -#' data(vdendool) +#' gph <- terra::unwrap(vdendool) #' -#' nh_modes <- eot(x = vdendool, y = NULL, n = 3, +#' nh_modes <- eot(x = gph, y = NULL, n = 3, #' standardised = FALSE, #' verbose = TRUE) #' @@ -41,14 +43,24 @@ setMethod('nXplain', signature(x = 'EotStack'), x[[i]]@cum_exp_var }) - n <- min(which(var - expl.var <= 0), na.rm = TRUE) - - if (!is.finite(n)) { - stop("explained variance of EotStack is lower than: ", var, - "\n", - "maximum explained variance of this EotStack is: ", - x[[nmodes(x)]]@cum_exp_var) + idx = var - expl.var <= 0 + + if (!any(idx)) { + paste( + "explained variance of EotStack is lower than: %s" + , "maximum explained variance of this EotStack is: %s" + , sep = "\n" + ) |> + sprintf( + var + , x[[nmodes(x)]]@cum_exp_var + ) |> + stop( + call. = FALSE + ) } + + n <- min(which(idx), na.rm = TRUE) return(n) } diff --git a/R/names.R b/R/names.R index 95db145..ec0cda7 100644 --- a/R/names.R +++ b/R/names.R @@ -15,9 +15,9 @@ if (!isGeneric('names')) { #' if \code{x} is a EotMode, the name the respective mode #' #' @examples -#' data(vdendool) +#' gph <- terra::unwrap(vdendool) #' -#' nh_modes <- eot(vdendool, n = 2) +#' nh_modes <- eot(gph, n = 2) #' #' ## mode names #' names(nh_modes) diff --git a/R/nmodes.R b/R/nmodes.R index f54277d..d63148e 100644 --- a/R/nmodes.R +++ b/R/nmodes.R @@ -1,32 +1,19 @@ -if (!isGeneric('nmodes')) { - setGeneric('nmodes', function(x, ...) - standardGeneric('nmodes')) -} - -#' Number of modes of an EotStack +#' Number of modes of an `EotStack` #' -#' @param x an EotStack +#' @description +#' Retrieves the number of modes of an `EotStack`. #' -#' @details -#' retrieves the number of modes of an EotStack +#' @param x An `EotStack`. #' #' @return -#' integer +#' The number of modes as `integer`. #' #' @examples -#' data(vdendool) -#' -#' nh_modes <- eot(vdendool, n = 2) -#' +#' gph = terra::unwrap(vdendool) +#' nh_modes = eot(gph, n = 2) #' nmodes(nh_modes) #' #' @export -#' @name nmodes -#' @rdname nmodes -#' @aliases nmodes,EotStack-method - -setMethod('nmodes', signature(x = 'EotStack'), - function(x) { - length(x@modes) - } -) +nmodes = function(x) { + length(x@modes) +} diff --git a/R/plot.R b/R/plot.R index ce03764..71c828d 100644 --- a/R/plot.R +++ b/R/plot.R @@ -32,13 +32,13 @@ if ( !isGeneric('plot') ) { #' @param locations logical. If x is an EotStack, set this to TRUE to #' produce a map showing the locations of all modes. Ignored if x is an #' EotMode -#' @param ... further arguments to be passed to [raster::spplot()] +#' @param ... further arguments to be passed to [sp::spplot()] #' #' @examples -#' data(vdendool) +#' gph <- terra::unwrap(vdendool) #' #' ## claculate 2 leading modes -#' nh_modes <- eot(x = vdendool, y = NULL, n = 2, +#' nh_modes <- eot(x = gph, y = NULL, n = 2, #' standardised = FALSE, #' verbose = TRUE) #' @@ -87,6 +87,8 @@ setMethod('plot', signature(x = 'EotMode', locations = FALSE, ...) { + arrange = match.arg(arrange) + pkgs <- c("lattice", "latticeExtra", "grid", "gridExtra", "RColorBrewer", "maps") tst <- sapply(pkgs, "requireNamespace", @@ -110,9 +112,9 @@ setMethod('plot', signature(x = 'EotMode', rs <- slot(x, r.prm) if (is.null(ts.vec)) - ts.vec <- seq(raster::nlayers(x@resid_response)) + ts.vec <- seq(terra::nlyr(x@resid_response)) - xy <- raster::xyFromCell(x@rsq_predictor, + xy <- terra::xyFromCell(x@rsq_predictor, cell = x@cell_bp) mode.location.p <- lattice::xyplot(xy[1, 2] ~ xy[1, 1], @@ -133,13 +135,13 @@ setMethod('plot', signature(x = 'EotMode', x <- ifelse((x < 1) | (x > 359), NA, x) }) - if (max(extent(ps)@xmax) > 180) { + if (max(terra::xmax(ps)) > 180) { mm.pred <- mm360 } else { mm.pred <- mm180 } - if (max(extent(rs)@xmax) > 180) { + if (max(terra::xmax(rs)) > 180) { mm.resp <- mm360 } else { mm.resp <- mm180 @@ -151,8 +153,8 @@ setMethod('plot', signature(x = 'EotMode', } - px.pred <- raster::ncell(ps) - px.resp <- raster::ncell(rs) + px.pred <- terra::ncell(ps) + px.resp <- terra::ncell(rs) pred.p <- sp::spplot(ps, mm = mm.pred, maxpixels = px.pred, @@ -219,7 +221,6 @@ setMethod('plot', signature(x = 'EotMode', } ### set layout to wide or long - arrange <- arrange[1] if (arrange == "wide") ncls <- 2 else ncls <- 1 ### amalgamate pred.p and resp.p according to layout @@ -261,6 +262,8 @@ setMethod('plot', signature(x = 'EotStack', locations = FALSE, ...) { + arrange = match.arg(arrange) + if (missing(y)) y <- 1 if (!locations) { @@ -281,7 +284,7 @@ setMethod('plot', signature(x = 'EotStack', ...) } else { - remote::plotLocations(x, ...) + plotLocations(x, ...) } } @@ -313,9 +316,9 @@ setMethod('plot', signature(x = 'EotStack', # rs <- slot(x[[y]], r.prm) # # if (is.null(ts.vec)) -# ts.vec <- seq(raster::nlayers(x[[y]]@resid_response)) +# ts.vec <- seq(terra::nlyr(x[[y]]@resid_response)) # -# xy <- raster::xyFromCell(x[[y]]@rsq_predictor, +# xy <- terra::xyFromCell(x[[y]]@rsq_predictor, # cell = x[[y]]@cell_bp) # # mode.location.p <- lattice::xyplot(xy[1, 2] ~ xy[1, 1], @@ -354,8 +357,8 @@ setMethod('plot', signature(x = 'EotStack', # } # # -# px.pred <- raster::ncell(ps) -# px.resp <- raster::ncell(rs) +# px.pred <- terra::ncell(ps) +# px.resp <- terra::ncell(rs) # # pred.p <- sp::spplot(ps, # mm = mm.pred, maxpixels = px.pred, @@ -456,7 +459,7 @@ setMethod('plot', signature(x = 'EotStack', # ) -# definde function -------------------------------------------------------- +# define function -------------------------------------------------------- plotLocations <- function(x, ...) { pkgs <- c("lattice", "latticeExtra", "grid", "gridExtra", @@ -472,7 +475,7 @@ plotLocations <- function(x, ...) { ### plot function loc.df <- as.data.frame(do.call("rbind", lapply(seq(nmodes(x)), function(i) { - raster::xyFromCell( + terra::xyFromCell( x[[i]]@rsq_predictor, cell = x[[i]]@cell_bp) }))) @@ -481,7 +484,7 @@ plotLocations <- function(x, ...) { sep = "_") mm <- maps::map("world", plot = FALSE, fill = TRUE) - px.pred <- raster::ncell(x[[1]]@r_predictor) + px.pred <- terra::ncell(x[[1]]@r_predictor) pred.p <- sp::spplot(x[[1]]@rsq_predictor, mm = mm, maxpixels = px.pred, @@ -568,7 +571,8 @@ plotLocations <- function(x, ...) { } grid::upViewport(0) - + return(invisible()) + } else { stop("need packages 'gridExtra', 'latticeExtra' & 'maps' to plot locations") } diff --git a/R/predict.R b/R/predict.R index b3e2f3f..0f7fca1 100644 --- a/R/predict.R +++ b/R/predict.R @@ -1,141 +1,146 @@ -# set methods ------------------------------------------------------------- -if ( !isGeneric('predict') ) { - setGeneric('predict', function(object, ...) - standardGeneric('predict')) +if ( !methods::isGeneric('predict') ) { + methods::setGeneric( + 'predict' + , function(object, ...) { + standardGeneric('predict') + } + ) } #' EOT based spatial prediction #' #' @description #' Make spatial predictions using the fitted model returned by -#' [eot()]. A (user-defined) set of \emph{n} modes will be used to +#' [eot()]. A (user-defined) set of 'n' modes will be used to #' model the outcome using the identified link functions of the respective modes #' which are added together to produce the final prediction. #' -#' @param object an \code{Eot*} object -#' @param newdata the data to be used as predictor -#' @param n the number of modes to be used for the prediction. +#' @param object An `Eot*` object. +#' @param newdata The data to be used as predictor. +#' @param n The number of modes to be used for the prediction. #' See [nXplain()] for calculating the number of modes based #' on their explanatory power. #' @param filename `character`, output filenames (optional). If specified, -#' this must be of the same length as \code{nlayers(newdata)}. -#' @param cores `integer`. Number of cores for parallel processing. -#' @param ... further arguments passed to [raster::calc()], and hence, -#' [raster::writeRaster()]. +#' this must be of the same length as `nlayers(newdata)`. +#' @param ... Further arguments passed to [terra::app()], and hence, +#' [terra::writeRaster()]. #' #' @return -#' a \emph{RasterStack} of \code{nlayers(newdata)} +#' A `SpatRaster` of `nlayers(newdata)`. #' #' @seealso -#' [raster::calc()], [raster::writeRaster()]. +#' [terra::app()], [terra::writeRaster()]. #' #' @examples #' ### not very useful, but highlights the workflow #' \donttest{ -#' data(pacificSST) -#' data(australiaGPCP) +#' sst = terra::unwrap(pacificSST) +#' pcp = terra::unwrap(australiaGPCP) #' #' ## train data using eot() -#' train <- eot(x = pacificSST[[1:10]], -#' y = australiaGPCP[[1:10]], +#' train <- eot(x = sst[[1:10]], +#' y = pcp[[1:10]], #' n = 1) #' #' ## predict using identified model #' pred <- predict(train, -#' newdata = pacificSST[[11:20]], +#' newdata = sst[[11:20]], #' n = 1) #' #' ## compare results #' opar <- par(mfrow = c(1,2)) -#' plot(australiaGPCP[[13]], main = "original", zlim = c(0, 10)) +#' plot(pcp[[13]], main = "original", zlim = c(0, 10)) #' plot(pred[[3]], main = "predicted", zlim = c(0, 10)) #' par(opar) #' } #' @export #' @name predict -#' @rdname predict -#' @aliases predict,EotStack-method - -setMethod('predict', signature(object = 'EotStack'), - function(object, - newdata, - n = 1, - cores = 1L, - filename = '', - ...) { - - ### extract identified EOT (@cell_bp) - bps <- sapply(seq(n), function(i) object[[i]]@cell_bp) - ts.modes <- t(raster::extract(newdata, bps)) - - ### target files and parallelization - vld <- length(filename) == raster::nlayers(newdata) - filename <- if (vld) filename else rep("", nrow(ts.modes)) - - dots <- list(...) - - cl <- parallel::makePSOCKcluster(cores) - on.exit(parallel::stopCluster(cl)) - parallel::clusterExport(cl, c("ts.modes", "object", "filename", - "dots"), envir = environment()) - - ### prediction using calculated intercept, slope and values - raster::stack( - parallel::parLapply(cl, seq(nrow(ts.modes)), function(i) { - rst <- raster::stack(lapply(seq(ncol(ts.modes)), function(k) { - object[[k]]@int_response + - object[[k]]@slp_response * ts.modes[i, k] - })) - - ### summate prediction for each mode at each time step - dots_sub <- list(x = rst, fun = sum, filename = filename[i]) - dots_sub <- append(dots, dots_sub) - - do.call(raster::calc, args = dots_sub) - }) - ) - } +# set methods ------------------------------------------------------------- +#' @aliases predict,EotStack-method +#' @rdname predict +methods::setMethod( + 'predict' + , signature(object = 'EotStack') + , function( + object + , newdata + , n = 1L + , filename = '' + , ...) { + + if (n == 1L) { + return( + predict( + object[[1L]] + , newdata = newdata + , n = 1L + , filename = filename + , ... + ) + ) + } + + ### extract identified EOT (@cell_bp) + bps <- sapply(seq(n), function(i) object[[i]]@cell_bp) + ts.modes <- t(terra::extract(newdata, bps)) + + ### target files + vld <- length(filename) == terra::nlyr(newdata) + filename <- if (vld) filename else rep("", nrow(ts.modes)) + + dots <- list(...) + + ### prediction using calculated intercept, slope and values + terra::rast( + lapply(seq(nrow(ts.modes)), function(i) { + + rst <- terra::rast(lapply(seq(ncol(ts.modes)), function(k) { + object[[k]]@int_response + + object[[k]]@slp_response * ts.modes[i, k] + })) + + ### summate prediction for each mode at each time step + dots_sub <- list(x = rst, fun = sum, filename = filename[i]) + dots_sub <- append(dots, dots_sub) + + do.call(terra::app, args = dots_sub) + } + ) + ) + } ) #' @aliases predict,EotMode-method #' @rdname predict -setMethod('predict', signature(object = 'EotMode'), - function(object, - newdata, - n = 1, - cores = 1L, - filename = '', - ...) { - - ### extract identified EOT (@cell_bp) - bps <- object@cell_bp - ts.modes <- t(raster::extract(newdata, bps)) - - ### target files and parallelization - vld <- length(filename) == raster::nlayers(newdata) - filename <- if (vld) filename else rep("", nrow(ts.modes)) - - dots <- list(...) - - cl <- parallel::makePSOCKcluster(cores) - on.exit(parallel::stopCluster(cl)) - parallel::clusterExport(cl, c("ts.modes", "object", "filename", - "dots"), envir = environment()) - - ### prediction using claculated intercept, slope and values - raster::stack( - parallel::parLapply(cl, seq(nrow(ts.modes)), function(i) { - - rst <- raster::overlay(object@int_response, object@slp_response, - fun = function(x, y) x + y * ts.modes[i, ]) - - ### summate prediction for each mode at each time step - dots_sub <- list(x = rst, fun = sum, filename = filename[i]) - dots_sub <- append(dots, dots_sub) - - do.call(raster::calc, args = dots_sub) - }) - ) - } +methods::setMethod( + 'predict' + , signature(object = 'EotMode') + , function( + object + , newdata + , n = 1L + , filename = '' + , ... + ) { + + ### extract identified EOT (@cell_bp) + bps <- object@cell_bp + ts.modes <- t(terra::extract(newdata, bps)) + + ### target files + vld <- length(filename) == terra::nlyr(newdata) + filename <- if (vld) filename else rep("", nrow(ts.modes)) + + dots <- list(...) + + ### prediction using claculated intercept, slope and values + terra::rast( + lapply(seq(nrow(ts.modes)), function(i) { + + rst <- terra::app(c(object@int_response, object@slp_response), + fun = function(x) x[[1L]] + x[[2L]] * ts.modes[i, ]) + }) + ) + } ) diff --git a/R/readEot.R b/R/readEot.R index 8493e11..857808f 100644 --- a/R/readEot.R +++ b/R/readEot.R @@ -1,49 +1,71 @@ -#' Read \code{Eot}* files from disk +#' Read `Eot*` files from disk #' #' @description -#' Read \code{Eot}* related files from disk, e.g. for further use with +#' Read `Eot*` related files from disk, e.g. for further use with #' [remote::predict()] or [remote::plot()]. #' -#' @param x `character`, search path for \code{Eot}* related files passed +#' @param x `character`, search path for `Eot*` related files passed #' to [list.files()]. -#' @param prefix `character`, see [writeEot()] for details. -#' Should be the same as previously supplied to [eot()]. -#' @param suffix `character`, file extension depending on the output file -#' type of locally stored \code{Eot}* files, see [raster::writeRaster()]. +#' @param prefix `character`, see [writeEot()] for details. Should be the same +#' 'prefix' as used during file creation in [eot()] or [writeEot()]. +#' @param suffix `character`, default `.grd` for native `"RRASTER"` format. File +#' extension depending on the output file type of locally stored `Eot*` files. #' -#' @return An \code{Eot}* object. +#' @return An `EotMode` if a single mode is found on disk, or an `EotStack` +#' otherwise. #' -#' @seealso [eot()], [writeEot()], -#' [raster::writeRaster()]. +#' @seealso [writeEot()] #' #' @author Florian Detsch #' #' @examples #' \dontrun{ #' ## calculate 3 leading modes -#' data(vdendool) -#' nh_modes <- eot(x = vdendool, n = 3, standardised = FALSE, -#' write.out = TRUE, path.out = "~/data") +#' gph <- terra::unwrap(vdendool) +#' nh_modes <- eot(x = gph, n = 3, standardised = FALSE, +#' write.out = TRUE, path.out = tempdir()) #' #' ## reimport related files #' rm(nh_modes) -#' nh_modes <- readEot("~/data") +#' nh_modes <- readEot(tempdir()) #' nh_modes #' } #' -#' @export readEot -#' @name readEot -readEot <- function(x, prefix = "remote", suffix = "grd") { +#' @export +readEot <- function(x, prefix = "remote", suffix = ".grd") { ## identify available files and modes - fls_mds <- list.files(x, pattern = paste0(prefix, "_mode.*", suffix), - full.names = TRUE) - mds <- unique(sapply(strsplit(fls_mds, "_"), "[[", 3)) - + fls_mds = list.files( + x + , pattern = sprintf( + "^%s_mode_\\d+_.*_(predictor|response)%s$" + , prefix + , suffix + ) + , full.names = TRUE + ) + + mds = regmatches( + fls_mds + , m = regexpr( + "(?<=mode_)\\d+" + , fls_mds + , perl = TRUE + ) + ) |> + unique() + ## import locations and explained variance related to leading modes - dat_mds <- list.files(x, paste0(prefix, "_eot_locations.csv"), - full.names = TRUE) - dat_mds <- utils::read.csv(dat_mds) + # TODO: `.csv` file is only created in `eot()`, not in `writeEot()` + dat_mds = list.files( + x + , sprintf( + "^%s_eot_locations\\.csv$" + , prefix + ) + , full.names = TRUE + ) |> + utils::read.csv() ## loop over modes, creating 'EotMode' objects for each mode available lst_eot <- lapply(mds, function(n) { @@ -55,7 +77,7 @@ readEot <- function(x, prefix = "remote", suffix = "grd") { # import files lst <- lapply(1:length(fls), function(j) { - if (j %in% c(7, 13)) raster::brick(fls[j]) else raster::raster(fls[j]) + terra::rast(fls[j]) }) # create 'EotMode' object @@ -82,16 +104,15 @@ readEot <- function(x, prefix = "remote", suffix = "grd") { p_response = lst[[12]], resid_response = lst[[13]]) }) - - ## create 'EotStack' objects if more than one leading mode is available - if (length(mds) > 1) { - names(lst_eot) <- sapply(lst_eot, function(i) i@name) - new('EotStack', modes = lst_eot, names = names(lst_eot)) - - ## else return 'EotMode' object - } else { - lst_eot[[1]] + + ## for a single leading mode, return `EotMode` + if (length(mds) == 1L) { + return(lst_eot[[1L]]) } + + ## else create an `EotStack` if more than one leading mode is available + names(lst_eot) <- sapply(lst_eot, function(i) i@name) + new('EotStack', modes = lst_eot, names = names(lst_eot)) } # function to create patterns of required raster* files diff --git a/R/remote-package.R b/R/remote-package.R index 8989082..7fa79db 100644 --- a/R/remote-package.R +++ b/R/remote-package.R @@ -22,7 +22,7 @@ #' Oxford University Press, Oxford, New York (2007)\cr #' @seealso \pkg{remote} is built upon Raster* classes from the #' [raster::raster-package]. Please see their documentation for data preparation etc. -#' @import Rcpp raster gridExtra latticeExtra mapdata scales methods parallel +#' @import Rcpp raster gridExtra latticeExtra mapdata scales methods parallel terra #' @importFrom grDevices colorRampPalette hcl #' @importFrom stats pt var cov.wt na.exclude princomp na.omit #' @importFrom utils read.csv write.table @@ -33,15 +33,17 @@ #' #' @docType data #' @name vdendool +#' @aliases vdendool #' @title Mean seasonal (DJF) 700 mb geopotential heights #' @description NCEP/NCAR reanalysis data of mean seasonal (DJF) 700 mb geopotential heights from 1948 to 1998 #' @details NCEP/NCAR reanalysis data of mean seasonal (DJF) 700 mb geopotential heights from 1948 to 1998 -#' @format a RasterBrick with the following attributes\cr -#' \cr -#' dimensions : 14, 36, 504, 50 (nrow, ncol, ncell, nlayers)\cr -#' resolution : 10, 4.931507 (x, y)\cr -#' extent : -180, 180, 20.9589, 90 (xmin, xmax, ymin, ymax)\cr -#' coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0\cr +#' @format A `PackedSpatRaster` with the following attributes:\cr +#' ```sh +#' size : 14, 36, 50 (nrow, ncol, nlyr) +#' resolution : 10, 5 (x, y) +#' extent : -180, 180, 20, 90 (xmin, xmax, ymin, ymax) +#' coord. ref. : +proj=longlat +datum=WGS84 +no_defs +#' ``` #' @references #' The NCEP/NCAR 40-year reanalysis project\cr #' Kalnay et al. (1996)\cr @@ -50,42 +52,65 @@ #' @source #' \cr #' \emph{Original Source:} NOAA National Center for Environmental Prediction +#' +#' @usage +#' vdendool +#' +#' @examples +#' terra::unwrap(vdendool) NULL + #' #' @docType data #' @name australiaGPCP +#' @aliases australiaGPCP #' @title Monthly GPCP precipitation data for Australia #' @description Monthly Gridded Precipitation Climatology Project precipitation data #' for Australia from 1982/01 to 2010/12 #' @details Monthly Gridded Precipitation Climatology Project precipitation data #' for Australia from 1982/01 to 2010/12 -#' @format a RasterBrick with the following attributes\cr -#' \cr -#' dimensions : 12, 20, 240, 348 (nrow, ncol, ncell, nlayers)\cr -#' resolution : 2.5, 2.5 (x, y)\cr -#' extent : 110, 160, -40, -10 (xmin, xmax, ymin, ymax)\cr -#' coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs\cr +#' @format A `PackedSpatRaster` with the following attributes:\cr +#' ```sh +#' size : 12, 20, 348 (nrow, ncol, nlyr) +#' resolution : 2.5, 2.5 (x, y) +#' extent : 110, 160, -40, -10 (xmin, xmax, ymin, ymax) +#' coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs +#' ``` #' @references #' The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979 - Present)\cr #' Adler et al. (2003)\cr #' Journal of Hydrometeorology, Volume 4, Issue 6, pp. 1147 - 1167\cr #' \doi{10.1175/1525-7541(2003)004<1147:TVGPCP>2.0.CO;2} +#' +#' @usage +#' australiaGPCP +#' +#' @examples +#' terra::unwrap(australiaGPCP) NULL #' #' @docType data #' @name pacificSST +#' @aliases pacificSST #' @title Monthly SSTs for the tropical Pacific Ocean #' @description Monthly NOAA sea surface temperatures for the tropical Pacific Ocean from 1982/01 to 2010/12 #' @details Monthly NOAA sea surface temperatures for the tropical Pacific Ocean from 1982/01 to 2010/12 -#' @format a RasterBrick with the following attributes\cr -#' \cr -#' dimensions : 30, 140, 4200, 348 (nrow, ncol, ncell, nlayers)\cr -#' resolution : 1, 1 (x, y)\cr -#' extent : 150, 290, -15, 15 (xmin, xmax, ymin, ymax)\cr -#' coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs\cr +#' @format A `PackedSpatRaster` with the following attributes:\cr +#' ```sh +#' size : 30, 140, 348 (nrow, ncol, nlyr) +#' resolution : 1, 1 (x, y) +#' extent : 150, 290, -15, 15 (xmin, xmax, ymin, ymax) +#' coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs +#' ``` #' @references #' Daily High-Resolution-Blended Analyses for Sea Surface Temperature\cr #' Reynolds et al. (2007)\cr #' Journal of Climate, Volume 20, Issue 22, pp. 5473 - 5496\cr #' \doi{10.1175/2007JCLI1824.1} +#' +#' @usage +#' pacificSST +#' +#' @examples +#' terra::unwrap(pacificSST) NULL diff --git a/R/show.R b/R/show.R index 1be3747..dd7d2be 100644 --- a/R/show.R +++ b/R/show.R @@ -1,52 +1,54 @@ -setMethod ('show' , 'EotMode', +methods::setMethod ('show' , 'EotMode', function(object) { cat('class :', class(object), '\n') cat('name :', names(object), '\n') cat('base point (x, y) :', object@coords_bp, '\n') cat('cum. expl. variance :', object@cum_exp_var, '\n') cat('dimensions : ', - raster::nrow(object@r_predictor), ', ', - raster::ncol(object@r_predictor), ', ', - raster::ncell(object@r_predictor), + terra::nrow(object@r_predictor), ', ', + terra::ncol(object@r_predictor), ', ', + terra::ncell(object@r_predictor), ' (nrow, ncol, ncell)\n', sep="") cat('resolution : ', - raster::xres(object@r_predictor), ', ', - raster::yres(object@r_predictor), ' (x, y)\n', sep="") - cat('extent : ', object@r_predictor@extent@xmin, - ', ', object@r_predictor@extent@xmax, ', ', - object@r_predictor@extent@ymin, ', ', - object@r_predictor@extent@ymax, + terra::xres(object@r_predictor), ', ', + terra::yres(object@r_predictor), ' (x, y)\n', sep="") + xtnt = as.vector(terra::ext(object@r_predictor)) + cat('extent : ', xtnt['xmin'], + ', ', xtnt['xmax'], ', ', + xtnt['ymin'], ', ', xtnt['ymax'], ' (xmin, xmax, ymin, ymax)\n', sep="") cat('coord. ref. :', - raster::projection(object@r_predictor, TRUE), '\n') + terra::crs(object@r_predictor, proj = TRUE), '\n') } -) +) -setMethod ('show' , 'EotStack', +methods::setMethod ('show' , 'EotStack', function(object) { + obj1 = object[[1L]] cat('class :', class(object), '\n') cat('cum. expl. variance :', object[[nmodes(object)]]@cum_exp_var, '\n') cat('names :', names(object), '\n') cat('dimensions : ', - raster::nrow(object[[1]]@r_predictor), ', ', - raster::ncol(object[[1]]@r_predictor), ', ', - raster::ncell(object[[1]]@r_predictor), ', ', + terra::nrow(obj1@r_predictor), ', ', + terra::ncol(obj1@r_predictor), ', ', + terra::ncell(obj1@r_predictor), ', ', nmodes(object), ' (nrow, ncol, ncell, nmodes)\n', sep="") cat('resolution : ', - raster::xres(object[[1]]@r_predictor), ', ', - raster::yres(object[[1]]@r_predictor), + terra::xres(obj1@r_predictor), ', ', + terra::yres(obj1@r_predictor), ' (x, y)\n', sep="") + xtnt = as.vector(terra::ext(obj1@r_predictor)) cat('extent : ', - object[[1]]@r_predictor@extent@xmin, ', ', - object[[1]]@r_predictor@extent@xmax, ', ', - object[[1]]@r_predictor@extent@ymin, ', ', - object[[1]]@r_predictor@extent@ymax, + xtnt['xmin'], ', ', + xtnt['xmax'], ', ', + xtnt['ymin'], ', ', + xtnt['ymax'], ' (xmin, xmax, ymin, ymax)\n', sep="") cat('coord. ref. :', - raster::projection(object[[1]]@r_predictor, TRUE), '\n') + terra::crs(obj1@r_predictor, proj = TRUE), '\n') } ) \ No newline at end of file diff --git a/R/subset.R b/R/subset.R index 237fa6c..8a51ff3 100644 --- a/R/subset.R +++ b/R/subset.R @@ -3,24 +3,26 @@ if (!isGeneric('subset')) { standardGeneric('subset')) } -#' Subset modes in EotStacks +#' Subset modes in an `EotStack` #' #' @description -#' Extract a set of modes from an EotStack +#' Extract a set of modes from an `EotStack`. #' -#' @param x EotStack to be subset -#' @param subset integer or character. The modes to ectract (either by -#' integer or by their names) -#' @param drop if `TRUE` a single mode will be returned as an EotMode -#' @param ... currently not used +#' @param x `EotStack` to be subset +#' @param subset `integer` or `character`. The modes to extract (either by +#' their indexes or names). +#' @param drop If `TRUE`, a single selected mode is returned as an `EotMode`. +#' Defaults to `FALSE`, which always returns an `EotStack`. +#' @param ... Currently not used. #' #' @return -#' an Eot* object +#' An `EotMode` if a single mode is selected and `drop = TRUE`, otherwise an +#' `EotStack`. #' #' @examples -#' data(vdendool) +#' gph <- terra::unwrap(vdendool) #' -#' nh_modes <- eot(x = vdendool, y = NULL, n = 3, +#' nh_modes <- eot(x = gph, y = NULL, n = 3, #' standardised = FALSE, #' verbose = TRUE) #' @@ -33,6 +35,10 @@ if (!isGeneric('subset')) { #' subs <- subset(nh_modes, 2, drop = TRUE) #' class(subs) #' +#' ## similarly to `drop = TRUE` above: +#' nh_modes[[2L]] +#' nh_modes[["mode_02"]] +#' #' @export subset #' @name subset @@ -51,6 +57,8 @@ setMethod('subset', signature(x = 'EotStack'), } subset <- i } + # TODO: inconsistent behavior: invalid mode names are omitted in the + # presence of valid names, but invalid indexes are not subset <- as.integer(subset) if (! all(subset %in% 1:nmodes(x))) { stop('not a valid subset') @@ -65,7 +73,7 @@ setMethod('subset', signature(x = 'EotStack'), ) #' @rdname subset -#' @param i number of EotMode to be subset +#' @param i Index(es) or name(s) to be subset, delegated to [remote::subset()]. setMethod("[[", signature(x = "EotStack"), function(x, i) { diff --git a/R/writeEot.R b/R/writeEot.R index 8133f25..317dde5 100644 --- a/R/writeEot.R +++ b/R/writeEot.R @@ -1,38 +1,44 @@ -if (!isGeneric('writeEot')) { - setGeneric('writeEot', function(x, ...) - standardGeneric('writeEot')) -} +methods::setGeneric( + 'writeEot' + , function(x, ...) { + standardGeneric('writeEot') + } +) -#' Write Eot* objects to disk +#' Write `Eot*` objects to disk #' #' @description -#' Write Eot* objects to disk. This is merely a wrapper around -#' [raster::writeRaster()] so see respective help section for details. +#' Write `Eot*` objects to disk. This is merely a wrapper around +#' [terra::writeRaster()] so see respective help section for details. #' -#' @param x an Eot* object -#' @param path.out the path to the folder to write the files to -#' @param prefix a prefix to be added to the file names (see Details) -#' @param overwrite see [raster::writeRaster()]. -#' Defaults to `TRUE` in [writeEot()] -#' @param ... further arguments passed to [raster::writeRaster()] +#' @param x An `Eot*` object. +#' @param path.out The path to the folder to write the files to. +#' @param prefix A prefix to be added to the file names (see Details). +#' @param overwrite See [terra::writeRaster()]. In [writeEot()], this defaults +#' to `TRUE`. +#' @param ... Further arguments passed to [terra::writeRaster()] for `EotMode` +#' input, or to the underlying `EotMode` method for `EotStack` input. If +#' 'filetype' is not specified, the default `"RRASTER"` format (`.grd`) is +#' used. #' #' @details -#' [writeEot()] will write the results of either an EotMode or an EotStack -#' to disk. For each mode the following files will be written: +#' [writeEot()] will write the results of either an `EotMode` or an `EotStack` +#' to disk. For each mode the following files will be written: #' #' \itemize{ -#' \item \emph{pred_r} - the \emph{RasterLayer} of the correlation coefficients -#' between the base point and each pixel of the predictor domain +#' \item \emph{pred_r} - the `SpatRaster` of the correlation coefficients +#' between the base point and each pixel of the predictor domain #' \item \emph{pred_rsq} - as above but for the coefficient of determination -#' \item \emph{pred_rsq_sums} - as above but for the sums of coefficient of determination -#' \item \emph{pred_int} - the \emph{RasterLayer} of the intercept of the -#' regression equation for each pixel of the predictor domain +#' \item \emph{pred_rsq_sums} - as above but for the sums of coefficient of +#' determination +#' \item \emph{pred_int} - the `SpatRaster` of the intercept of the +#' regression equation for each pixel of the predictor domain #' \item \emph{pred_slp} - same as above but for the slope of the -#' regression equation for each pixel of the predictor domain -#' \item \emph{pred_p} - the \emph{RasterLayer} of the significance (p-value) -#' of the the regression equation for each pixel of the predictor domain -#' \item \emph{pred_resid} - the \emph{RasterBrick} of the reduced data -#' for the predictor domain +#' regression equation for each pixel of the predictor domain +#' \item \emph{pred_p} - the `SpatRaster` of the significance (p-value) +#' of the the regression equation for each pixel of the predictor domain +#' \item \emph{pred_resid} - the `SpatRaster` of the reduced data +#' for the predictor domain #' } #' #' Apart from \emph{pred_rsq_sums}, all these files are also created for @@ -42,16 +48,16 @@ if (!isGeneric('writeEot')) { #' #' \emph{prefix_mode_n_pred_r.grd} #' -#' for the \emph{RasterLayer} of the predictor correlation coefficient -#' of mode n using the standard \emph{raster} file type (.grd). +#' for the `SpatRaster` of the predictor correlation coefficient +#' of mode n using the standard \pkg{raster} file type (`.grd`). #' -#' @seealso [raster::writeRaster()] +#' @seealso [terra::writeRaster()] #' #' @examples #' \dontrun{ -#' data(vdendool) +#' gph <- terra::unwrap(vdendool) #' -#' nh_modes <- eot(x = vdendool, y = NULL, n = 2, +#' nh_modes <- eot(x = gph, y = NULL, n = 2, #' standardised = FALSE, #' verbose = TRUE) #' @@ -68,48 +74,83 @@ if (!isGeneric('writeEot')) { # set methods ------------------------------------------------------------- -setMethod('writeEot', signature(x = 'EotMode'), - function(x, - path.out = ".", - prefix = "remote", - overwrite = TRUE, - ...) { - - out.name <- lapply(slotNames(x)[7:19], - function(i) { - paste(prefix, "mode", sprintf("%02.f", - x@mode), - i, sep = "_") - }) - - out.object <- lapply(slotNames(x)[7:19], function(j) { - slot(x, j) - }) - - a <- unlist(out.object) - b <- unlist(out.name) - - lapply(1:length(a), function(i) { - raster::writeRaster(a[[i]], paste(path.out, b[i], sep = "/"), - overwrite = overwrite, ...) - }) - } +methods::setMethod( + 'writeEot' + , signature(x = 'EotMode') + , function( + x + , path.out = "." + , prefix = "remote" + , overwrite = TRUE + , ... + ) { + + ## get slots of interest + slots = grep( + "^(r|rsq|int|slp|p|resid).*_(predictor|response)$" + , slotNames(x) + , value = TRUE + ) + + ## imitate `raster::writeRaster()` default behavior for file type and + ## extension if 'filetype' is not specified + dots = list(...) + fext = "" + + if (!"filetype" %in% names(dots)) { + dots$filetype = "RRASTER" + fext = ".grd" + } + + ## construct file names + out.name = sprintf( + "%s_mode_%02.f_%s%s" + , prefix + , x@mode + , slots + , fext + ) + + ## get objects to write + out.object = lapply(slots, slot, object = x) + + ## write objects to disk + Map( + \(oo, on) { + do.call( + terra::writeRaster + , args = c( + list( + x = oo + , filename = file.path(path.out, on) + , overwrite = overwrite + ) + , dots + ) + ) + } + , out.object + , out.name + ) + } ) #' @describeIn writeEot EotStack #' @aliases writeEot,EotStack-method -setMethod('writeEot', signature(x = 'EotStack'), - function(x, - path.out = ".", - prefix, - ...) { - - for (i in seq(nmodes(x))) { - writeEot(x[[i]], - path.out = path.out, - prefix = prefix, - ...) - } - } -) +methods::setMethod( + 'writeEot' + , signature(x = 'EotStack') + , function( + x + , ... + ) { + + for (i in seq(nmodes(x))) { + writeEot( + x[[i]] + , ... + ) + } + } +) \ No newline at end of file diff --git a/data/australiaGPCP.RData b/data/australiaGPCP.RData deleted file mode 100644 index 73175fb..0000000 Binary files a/data/australiaGPCP.RData and /dev/null differ diff --git a/data/australiaGPCP.rda b/data/australiaGPCP.rda new file mode 100644 index 0000000..e27fd4d Binary files /dev/null and b/data/australiaGPCP.rda differ diff --git a/data/pacificSST.RData b/data/pacificSST.RData deleted file mode 100644 index 4b9c521..0000000 Binary files a/data/pacificSST.RData and /dev/null differ diff --git a/data/pacificSST.rda b/data/pacificSST.rda new file mode 100644 index 0000000..e4783f9 Binary files /dev/null and b/data/pacificSST.rda differ diff --git a/data/vdendool.RData b/data/vdendool.RData deleted file mode 100755 index ea1eecf..0000000 Binary files a/data/vdendool.RData and /dev/null differ diff --git a/data/vdendool.rda b/data/vdendool.rda new file mode 100644 index 0000000..e290423 Binary files /dev/null and b/data/vdendool.rda differ diff --git a/dev_history.R b/dev_history.R index 6b6f169..0af078a 100644 --- a/dev_history.R +++ b/dev_history.R @@ -19,3 +19,202 @@ txt = utils::bibentry( , "doi:10.18637/jss.v065.i10 ." ) ) + + +# 2026-06-21 ==== + +## document, check and build package +devtools::check( + document = TRUE # `devtools::document()` + , build_args = "--resave-data=best" + , manual = FALSE + , cran = TRUE + , run_dont_test = TRUE +) +pak::local_install(ask = FALSE) + +## bump version +if (!requireNamespace("oiseasy", quietly = TRUE)) { + remotes::install_git( + "https://codeberg.org/tim-salabim/oiseasy.git" + , Ncpus = 4L + ) +} + +oiseasy::bumpDevVersion() + + +# 2026-07-07 ==== + +## BUILT-IN DATA ==== + +## if missing, create target folders +for (i in c("data-raw", "inst/extdata")) { + dir.create(i, showWarnings = FALSE) +} + +to_tif = FALSE +to_rda = TRUE + +## process built-in datasets +for (rdata_name in c("australiaGPCP", "pacificSST", "vdendool")) { + + # # debug: + # rdata_name = "australiaGPCP" + + # construct file names + rdata_file = paste0( + rdata_name + , ".RData" + ) + + rdata_file_raw = file.path("data-raw", rdata_file) + rdata_file_data = file.path("data", rdata_file) + + # if applicable, move existing `.RData` file to `data-raw/` + if (file.exists(rdata_file_data)) { + file.rename( + rdata_file_data + , rdata_file_raw + ) + } + + # load built-in data into memory + load(rdata_file_raw) + + get(rdata_name) |> + terra::rast() |> + assign( + rdata_name + , value = _ + ) + + # write `.tif` to `inst/extdata` + if (to_tif) { + tif_file_ext = sprintf("inst/extdata/%s.tif", rdata_name) + dir.create(dirname(tif_file_ext), showWarnings = FALSE, recursive = TRUE) + + terra::writeRaster( + get(rdata_name) + , filename = tif_file_ext + , overwrite = TRUE + , gdal = c("COMPRESS=LZW") + ) + } + + # write `.rda` to `data/` + if (to_rda) { + rda_file_data = sprintf("data/%s.rda", rdata_name) + dir.create(dirname(rda_file_data), showWarnings = FALSE) + + get(rdata_name) |> + terra::wrap() |> + assign( + rdata_name + , value = _ + ) + + save( + list = rdata_name + , envir = environment() + , file = rda_file_data + , compress = "bzip2" + ) + } +} + +## verify new built-in data +load("data/australiaGPCP.rda") +terra::unwrap(australiaGPCP) + + +# 2026-07-13 ==== + +## {tinytest} SETUP ==== + +tinytest::setup_tinytest(pkgdir = ".") +tinytest::run_test_dir() + +covr::report() + + +# 2026-07-14 ==== + +## `calcVar()` SPEED-UP ==== + +x = terra::unwrap(pacificSST) + +## across time +microbenchmark::microbenchmark( + ref = vls <- mean(apply(terra::values(x), 1, var, na.rm = TRUE)) + , new = clc <- mean(terra::values(terra::app(x, "sd", na.rm = TRUE)^2)[, 1L]) + , times = 25L +) +# Unit: milliseconds +# expr min lq mean median uq max neval +# ref 46.488107 53.404365 93.521187 59.092277 62.971696 367.62212 25 +# new 6.915151 7.785462 8.564513 8.633283 9.261476 10.22619 25 + +tinytest::expect_equal( + clc + , target = vls +) + +## across space +microbenchmark::microbenchmark( + ref = vls <- mean(apply(terra::values(x), 2, var, na.rm = TRUE), na.rm = TRUE) + , new = glbl <- mean((terra::global(x, fun = "sd", na.rm = TRUE)[, 1L])^2) + , times = 25L +) +# Unit: milliseconds +# expr min lq mean median uq max neval +# ref 21.047098 26.973126 67.45445 29.584694 36.766237 345.3942 25 +# new 5.703815 6.513882 7.13055 7.042239 7.507108 11.4029 25 + +tinytest::expect_equal( + glbl + , target = vls +) + + +# 2026-07-15 ==== + +## `denoise()` BENCHMARKING ==== + +gph = terra::unwrap(vdendool) + +microbenchmark::microbenchmark( + terra = gph_dns <- denoise(gph, expl.var = 0.8, use.cpp = FALSE, verbose = F) + , rcpp = gph_dns_cpp <- denoise(gph, expl.var = 0.8, verbose = FALSE) +) +# Unit: milliseconds +# expr min lq mean median uq max neval +# terra 7.954648 8.832743 10.63336 9.91166 10.84664 19.66279 100 +# rcpp 8.732450 9.636186 11.87453 10.55196 11.59424 28.42883 100 + +tinytest::expect_equal( + gph_dns + , target = gph_dns_cpp +) + + +# 2026-07-16 ==== + +## `denoise()` / `covWeight()` WITH MISSING VALUES ==== + +pcp = terra::unwrap(australiaGPCP) + +set.seed(1899L) +idx = terra::spatSample(pcp, size = n, values = FALSE, cells = TRUE) + +pcp[idx] = NA_real_ +rng = global(pcp[[1L]], fun = "range", na.rm = TRUE) + +pcp_dns = denoise(pcp, expl.var = 0.8) + +opar = par(mfrow = c(1,2)) +plot(pcp[[1L]], main = "original") +plot(pcp_dns[[1L]], main = "denoised") +par(opar) + +plot(pcp[[1L]] - pcp_dns[[1L]], main = "residuals") diff --git a/inst/tinytest/test_anomalize.R b/inst/tinytest/test_anomalize.R new file mode 100644 index 0000000..1d34731 --- /dev/null +++ b/inst/tinytest/test_anomalize.R @@ -0,0 +1,45 @@ +pcp = terra::unwrap(australiaGPCP) +anm = anomalize(pcp) + +expect_inherits( + anm + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input" +) + +expect_identical( + dim(anm) + , target = dim(pcp) + , info = "returns a `SpatRaster` of the same dimensions as the input" +) + +expect_equal( + anomalize(as(pcp, "Raster")) + , target = anm + , info = "returns same result as `Raster*` input" +) + +## errors +expect_error( + anomalize(pcp, reference = "frozen_hearts") + , pattern = "^Expected 'reference' to inherit from .* but got" + , info = "throws error if 'reference' is not a raster" +) + +## warnings +ref = terra::app(pcp, fun = mean, na.rm = TRUE) + +expect_warning( + anm1 <- anomalize( + pcp + , reference = as(terra::rast(replicate(2L, ref)), "Raster") + ) + , pattern = "to have a single layer, but got .* Using the first layer only" + , info = "throws warning if 'reference' has more than one layer" +) + +expect_equal( + anm1 + , target = anm + , info = "returns same result as single-layer 'reference'" +) diff --git a/inst/tinytest/test_calcVar.R b/inst/tinytest/test_calcVar.R new file mode 100644 index 0000000..ee94e52 --- /dev/null +++ b/inst/tinytest/test_calcVar.R @@ -0,0 +1,36 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +sst = terra::unwrap(pacificSST) + +## non-standardized +var_ns = calcVar(sst) + +expect_number( + var_ns + , lower = 0 + , finite = TRUE + , info = "returns a single `numeric` value" +) + +## standardized +var_st = calcVar(sst, standardised = TRUE) + +expect_number( + var_st + , lower = 0 + , finite = TRUE + , info = "returns a single `numeric` value" +) + +expect_true( + var_st != var_ns +) + +## `Raster*` input +expect_identical( + calcVar(as(sst, "Raster")) + , target = var_ns + , info = "returns same result as `Raster*` input" +) diff --git a/inst/tinytest/test_covWeight.R b/inst/tinytest/test_covWeight.R new file mode 100644 index 0000000..d8e4694 --- /dev/null +++ b/inst/tinytest/test_covWeight.R @@ -0,0 +1,36 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +pcp = terra::unwrap(australiaGPCP) + +cov_wghts = covWeight( + terra::values(pcp) + , weights = getWeights(pcp) +) + +expect_list( + cov_wghts + , types = c( + "matrix" + , "numeric" + , "integer" + ) + , any.missing = FALSE + , len = 5L + , unique = TRUE + , names = "named" + , info = "returns a `list` of length 5 with named elements (see `?cov.wt`)" +) + +expect_names( + names(cov_wghts) + , identical.to = c( + "cov" + , "center" + , "n.obs" + , "wt" + , "cor" + ) +) + diff --git a/inst/tinytest/test_cutStack.R b/inst/tinytest/test_cutStack.R new file mode 100644 index 0000000..37c9446 --- /dev/null +++ b/inst/tinytest/test_cutStack.R @@ -0,0 +1,55 @@ +pcp = terra::unwrap(australiaGPCP) + +## `n = NULL` (default) +cut_dpl = cutStack(pcp) + +expect_equal( + cut_dpl + , target = pcp + , info = "returns unmodified raster series if `n == NULL`" +) + +## take away from end, e.g. if supplied series is predictor +n_nd = 8L +cut_nd = cutStack(pcp, tail = TRUE, n = n_nd) + +expect_inherits( + cut_nd + , class = "SpatRaster" + , info = "returns a 'SpatRaster' object" +) + +expect_true( + terra::nlyr(cut_nd) == (terra::nlyr(pcp) - n_nd) + , info = "removes 'n' layers from input (tail)" +) + +expect_equal( + cut_nd[[terra::nlyr(cut_nd)]] + , target = pcp[[terra::nlyr(pcp) - n_nd]] + , info = "last layer of cut series is correct" +) + +## take away from start, e.g. if supplied series is response +n_st = 6L +cut_st = cutStack(pcp, tail = FALSE, n = n_st) + +expect_true( + terra::nlyr(cut_st) == (terra::nlyr(pcp) - n_st) + , info = "removes 'n' layers from input (head)" +) + +expect_equal( + cut_st[[1L]] + , target = pcp[[n_st + 1L]] + , info = "first layer of cut series is correct" +) + +## same with `Raster*` input +cut_st_rst = cutStack(as(pcp, "Raster"), tail = FALSE, n = n_st) + +expect_equal( + cut_st_rst + , target = cut_st + , info = "`Raster*` input returns identical result" +) diff --git a/inst/tinytest/test_data.R b/inst/tinytest/test_data.R new file mode 100644 index 0000000..e0145d1 --- /dev/null +++ b/inst/tinytest/test_data.R @@ -0,0 +1,97 @@ +### `australiaGPCP` ---- + +n = 348L +prj = "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs" + +data( + "australiaGPCP" + , package = "remote" +) + +expect_inherits( + australiaGPCP + , "PackedSpatRaster" +) + +pcp = terra::unwrap(australiaGPCP) + +expect_true( + nrow(pcp) == 12L +) + +expect_true( + ncol(pcp) == 20L +) + +expect_true( + terra::nlyr(pcp) == n +) + +expect_identical( + terra::crs(pcp, proj = TRUE) + , target = prj +) + + +### `pacificSST` ---- + +data( + "pacificSST" + , package = "remote" +) + +expect_inherits( + pacificSST + , "PackedSpatRaster" +) + +sst = terra::unwrap(pacificSST) + +expect_true( + nrow(sst) == 30L +) + +expect_true( + ncol(sst) == 140L +) + +expect_true( + terra::nlyr(sst) == n +) + +expect_identical( + terra::crs(sst, proj = TRUE) + , target = prj +) + + +### `vdendool` ---- + +data( + "vdendool" + , package = "remote" +) + +expect_inherits( + vdendool + , "PackedSpatRaster" +) + +vdd = terra::unwrap(vdendool) + +expect_true( + nrow(vdd) == 14L +) + +expect_true( + ncol(vdd) == 36L +) + +expect_true( + terra::nlyr(vdd) == 50L +) + +expect_identical( + terra::crs(vdd, proj = TRUE) + , target = "+proj=longlat +datum=WGS84 +no_defs" +) diff --git a/inst/tinytest/test_deg2rad.R b/inst/tinytest/test_deg2rad.R new file mode 100644 index 0000000..2d46e39 --- /dev/null +++ b/inst/tinytest/test_deg2rad.R @@ -0,0 +1,43 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +gph = terra::unwrap(vdendool) + +## `numeric` input +degrees = terra::crds(gph)[, 2L] +radians = deg2rad(degrees) + +expect_numeric( + radians + , lower = -pi / 2 + , upper = pi / 2 + , finite = TRUE + , any.missing = FALSE + , len = length(degrees) + , info = "returns a `numeric` vector in case of `numeric` input" +) + +## `SpatRaster` input +tmp = gph[[1L]] +terra::values(tmp) = degrees + +tmp_rad = deg2rad(tmp) + +expect_inherits( + tmp_rad + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input" +) + +expect_identical( + dim(tmp_rad) + , target = dim(tmp) + , info = "returns a `SpatRaster` with the same dimensions as the input" +) + +expect_identical( + terra::values(tmp_rad)[, 1L] + , target = radians + , info = "returns same values as `numeric` input" +) diff --git a/inst/tinytest/test_denoise.R b/inst/tinytest/test_denoise.R new file mode 100644 index 0000000..4e53a3c --- /dev/null +++ b/inst/tinytest/test_denoise.R @@ -0,0 +1,57 @@ +pcp = terra::unwrap(australiaGPCP) +dns = denoise(pcp, expl.var = 0.8, verbose = FALSE) + +expect_inherits( + dns + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input" +) + +expect_identical( + dim(dns) + , target = dim(pcp) + , info = "returns a `SpatRaster` of the same dimensions as the input" +) + +expect_stdout( + dns_rst <- denoise(as(pcp, "Raster"), expl.var = 0.8, use.cpp = FALSE) + , pattern = paste( + "Using the first \\d+ components .* to reconstruct series" + , "these account for .* of variance in orig. series" + , sep = ".*" + ) + , info = "prints a message about the # of components and explained variance" +) + +expect_equal( + dns_rst + , target = dns + , info = "returns same result with `Raster*` input and `use.cpp = FALSE`" +) + +## unweighted, with 'k' instead of 'expl.var' +sst = terra::unwrap(pacificSST) + +expect_stdout( + dns_no_wghts <- denoise(sst, k = 5L, weighted = FALSE) + , pattern = "Using the first 5 components" + , info = "prints a message about the specified # of components" +) + +expect_inherits( + dns_no_wghts + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input (unweighted)" +) + +expect_identical( + dim(dns_no_wghts) + , target = dim(sst) + , info = "returns a `SpatRaster` of the same dimensions (unweighted)" +) + +## errors +expect_error( + denoise(pcp) + , pattern = "^Either 'expl.var' or 'k' must be supplied\\.$" +) diff --git a/inst/tinytest/test_deseason.R b/inst/tinytest/test_deseason.R new file mode 100644 index 0000000..be69f6f --- /dev/null +++ b/inst/tinytest/test_deseason.R @@ -0,0 +1,42 @@ +### `SpatRaster` ---- + +pcp = terra::unwrap(australiaGPCP) + +dsn = deseason(pcp, cycle.window = 12L) + +expect_inherits( + dsn + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input" +) + +expect_identical( + dim(dsn) + , target = dim(pcp) + , info = "returns a `SpatRaster` of the same dimensions as the input" +) + +## `Raster*` input +expect_equal( + deseason(as(pcp, "Raster"), cycle.window = 12L) + , target = dsn + , info = "returns same result as `Raster*` input" +) + +## with cpp +dsn_cpp = deseason(pcp, cycle.window = 12L, use.cpp = TRUE) + +expect_equal( + dsn_cpp + , target = dsn + , info = "returns same result as non-cpp implementation" +) + +## with file output +tmp = tempfile(fileext = ".tif") +jnk = deseason(pcp, cycle.window = 12L, filename = tmp) + +expect_true( + file.exists(tmp) + , info = "creates a file if 'filename' is specified" +) diff --git a/inst/tinytest/test_eot.R b/inst/tinytest/test_eot.R new file mode 100644 index 0000000..4cd0b80 --- /dev/null +++ b/inst/tinytest/test_eot.R @@ -0,0 +1,356 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +gph = terra::unwrap(vdendool) +n = 3L + + +### `eot()` ---- + +expect_stdout( + nh_modes <- eot( + x = gph + , y = NULL + , n = n + , standardised = FALSE + , verbose = TRUE + ) + , pattern = paste( + "Calculating linear model" + , "Locating \\d+\\. EOT" + , "Location" + , "Cum\\. expl\\. variance \\(\\%\\)" + , sep = ".*" + ) + , info = "prints progress messages to console" +) + +expect_inherits( + nh_modes + , class = "EotStack" + , info = "returns an `EotStack` object" +) + +expect_inherits( + nh_modes[[1L]] + , class = "EotMode" + , info = "returns an `EotMode` object per mode" +) + +## single-mode output with `write.out = TRUE` +nh_mode = eot( + x = gph + , y = NULL + , n = 1L + , standardised = FALSE + , write.out = TRUE + , path.out = tempdir() + , prefix = "dooln1" + , verbose = FALSE +) + +expect_equal( + nh_mode + , target = nh_modes[[1L]] + , info = "returns the same result for single-mode output" +) + +## `reduce.both = TRUE` with `write.out = TRUE` +nh_modes_rb <- eot( + x = gph + , y = NULL + , n = n + , standardised = FALSE + , write.out = TRUE + , path.out = tempdir() + , prefix = "dool" + , reduce.both = TRUE + , verbose = FALSE +) + +expect_true( + nh_modes_rb[[n]]@cum_exp_var != nh_modes[[n]]@cum_exp_var + , info = "returns different results for `reduce.both = TRUE`" +) + +nh_modes_rb_fls = list.files( + tempdir() + , pattern = "^dool_.*\\.grd$" +) + +expect_true( + length(nh_modes_rb_fls) > 0L && + length(nh_modes_rb_fls) %% n == 0L + , info = "writes EOT results to disk intrinsically if `write.out = TRUE`" +) + +## `Raster*` input +nh_modes_rst = eot( + x = as(gph, "Raster") + , y = NULL + , n = n + , standardised = FALSE + , verbose = FALSE +) + +expect_equal( + nh_modes_rst + , target = nh_modes + , info = "returns the same result for `Raster*` input" +) + + +### `readEot()` ---- + +nh_modes_rb_reimport = readEot( + x = tempdir() + , prefix = "dool" + , suffix = ".grd" +) + +expect_inherits( + nh_modes_rb_reimport + , class = "EotStack" + , info = "reimports EOT results with 2+ leading modes as `EotStack`" +) + +expect_equivalent( + nh_modes_rb_reimport + , target = nh_modes_rb + , info = "reimports EOT results from disk" +) + +## single leading mode +nh_mode_reimport = readEot( + x = tempdir() + , prefix = "dooln1" +) + +expect_inherits( + nh_mode_reimport + , class = "EotMode" + , info = "reimports EOT results with a single leading mode as `EotMode`" +) + + +### `names()` ---- + +## `EotStack` method +expect_character( + names(nh_modes) + , pattern = "^mode_\\d+$" + , any.missing = FALSE + , len = n + , unique = TRUE + , sorted = TRUE + , info = "sets default mode names" +) + +## `EotMode` method +expect_identical( + names(nh_modes[[1L]]) + , target = names(nh_modes)[1L] + , info = "returns the name of the respective mode" +) + +## discard mode names, i.e. set '' +names(nh_modes) = NULL + +expect_true( + all( + !nzchar(names(nh_modes)) + ) + , info = "accepts `NULL` to set empty mode names, i.e. ''" +) + +## set custom mode names +names(nh_modes) = paste0( + "vdendool" + , seq(n) +) + +expect_character( + names(nh_modes) + , pattern = "^vdendool\\d+$" + , any.missing = FALSE + , len = n + , unique = TRUE + , sorted = TRUE + , info = "accepts custom mode names" +) + +## error: lengths don't match +expect_error( + names(nh_modes) <- paste0( + "vdendool" + , seq(n + 1L) + ) + , pattern = "incorrect number of mode names" + , info = "throws an error if the number of mode names and modes don't match" +) + + +### `nmodes()` ---- + +expect_identical( + nmodes(nh_modes) + , target = n + , info = "returns the correct number of modes" +) + + +### `nXplain()` ---- + +expect_number( + nXplain(nh_modes, 0.25) + , lower = 1L + , upper = 3L + , finite = TRUE + , info = "returns a number between 1 (min. EOTs) and 3 ('n')" +) + +expect_error( + nXplain(nh_modes, 0.8) + , pattern = "explained variance of EotStack is lower than" + , info = "throws an error if the explained variance is lower than 'var'" +) + + +### `plot()` ---- + +expect_null( + plot(nh_modes, show.bp = TRUE) + , info = "returns `NULL` invisibly" +) + +expect_null( + plot(nh_modes, show.bp = TRUE, locations = TRUE) + , info = "returns `NULL` invisibly for plot of mode locations" +) + + +### `print()` ---- + +## `EotStack` +expect_stdout( + print(nh_modes) + , pattern = paste( + "^class .* EotStack" + , "cum. expl. variance" + , "names" + , "dimensions" + , "resolution" + , "extent" + , "coord. ref" + , sep = ".*" + ) + , info = "prints a summary of the `EotStack` object to console" +) + +## `EotMode` +expect_stdout( + print(nh_mode) + , pattern = "^class .* EotMode" + , info = "prints a summary of the `EotMode` object to console" +) + + +### `writeEot()` ---- + +## `filetype = "RRASTER"` +writeEot( + nh_modes + , prefix = "vdendool" + , path.out = tempdir() +) + +ofl = list.files( + tempdir() + , pattern = "^vdendool_mode_\\d+_.*_(predictor|response)\\.(grd|gri)$" + , full.names = TRUE +) + +expect_true( + length(ofl) > 0L && + length(ofl) %% n == 0L + , info = "writes EOT results to disk (default `.grd`)" +) + +## other 'filetype' +writeEot( + nh_modes + , prefix = "vdendool1" + , path.out = tempdir() + , filetype = "GTIFF" +) + +ofl1 = list.files( + tempdir() + , pattern = "^vdendool1_mode_\\d+_.*_(predictor|response)$" # no extension + , full.names = TRUE +) + +expect_true( + length(ofl1) > 0L && + length(ofl1) %% n == 0L + , info = "writes EOT results to disk (custom 'filetype')" +) + + +### `subset()` ---- + +## indexes +nh_modes_s2 = subset( + nh_modes + , subset = 2:3 +) + +expect_inherits( + nh_modes_s2 + , class = "EotStack" + , info = "returns a subset `EotStack`" +) + +expect_true( + remote::nmodes(nh_modes_s2) == 2L + , info = "returns a subset `EotStack` with the specified # of modes" +) + +## names +nh_modes_s1 = subset( + nh_modes + , subset = names(nh_modes)[3L] + , drop = TRUE +) + +expect_inherits( + nh_modes_s1 + , class = "EotMode" + , info = "returns a single subset `EotMode`" +) + +## errors and warnings +expect_error( + subset( + nh_modes + , subset = "flying_waters" + ) + , pattern = "invalid mode names" +) + +expect_warning( + subset( + nh_modes + , subset = c(names(nh_modes)[1L], "flying_waters") + ) + , pattern = "invalid mode names omitted" +) + +expect_error( + subset( + nh_modes + , subset = n:(n + 1L) + ) + , pattern = "not a valid subset" +) diff --git a/inst/tinytest/test_geoWeight.R b/inst/tinytest/test_geoWeight.R new file mode 100644 index 0000000..b9eca64 --- /dev/null +++ b/inst/tinytest/test_geoWeight.R @@ -0,0 +1,35 @@ +gph = terra::unwrap(vdendool) +gph_wghts = geoWeight(gph) + +expect_inherits( + gph_wghts + , class = "SpatRaster" + , info = "returns a `SpatRaster` in case of `SpatRaster` input" +) + +expect_identical( + dim(gph_wghts) + , target = dim(gph) + , info = "returns a `SpatRaster` of the same dimensions as the input" +) + +## different weighting function 'f' +gph_wghts_sqrt_cos = geoWeight( + gph + , f = \(x) sqrt(cos(x)) # downweights high latitudes less aggressively +) + +expect_true( + terra::global( + gph_wghts_sqrt_cos[[1L]] != gph_wghts[[1L]] + , fun = all + )[[1L]] + , info = "returns different results for different weighting function" +) + +## `Raster*` input +expect_equal( + geoWeight(as(gph, "Raster")) + , target = gph_wghts + , info = "returns same result as `Raster*` input" +) diff --git a/inst/tinytest/test_getWeights.R b/inst/tinytest/test_getWeights.R new file mode 100644 index 0000000..9737136 --- /dev/null +++ b/inst/tinytest/test_getWeights.R @@ -0,0 +1,39 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +pcp = terra::unwrap(australiaGPCP) + +wghts = getWeights(pcp) + +expect_numeric( + wghts + , lower = 0 + , upper = 1 + , finite = TRUE + , any.missing = FALSE + , len = terra::ncell(pcp) + , info = "returns a `numeric` vector of the same length as `ncell(x)`" +) + +## with `NA` +n = 5L + +set.seed(1899L) +idx = terra::spatSample(pcp, size = n, values = FALSE, cells = TRUE) + +pcp1 = pcp +pcp1[idx] = NA_real_ + +wghts_na = getWeights(as(pcp1, "Raster")) + +expect_true( + length(wghts_na) == terra::ncell(pcp) - n + , info = "returns a shorter weights vector in the presence of `NA` values" +) + +expect_identical( + wghts_na + , target = wghts[-idx] + , info = "returns the same weights as on original data, minus the `NA` values" +) diff --git a/inst/tinytest/test_lagalize.R b/inst/tinytest/test_lagalize.R new file mode 100644 index 0000000..6cb3b3d --- /dev/null +++ b/inst/tinytest/test_lagalize.R @@ -0,0 +1,49 @@ +library("tinytest") +library("checkmate") +using("checkmate") + +sst = terra::unwrap(pacificSST) +pcp = terra::unwrap(australiaGPCP) + +lagged = lagalize(sst, pcp, lag = 4L, freq = 12L) + +expect_list( + lagged + , types = "SpatRaster" + , any.missing = FALSE + , len = 2L + , unique = TRUE + , info = "returns a list of two `SpatRaster` objects" +) + +expect_true( + terra::nlyr(lagged[[1L]]) == terra::nlyr(lagged[[2L]]), + info = "returns two `SpatRaster` objects of equal length" +) + +expect_match( + names(lagged)[[1L]][1L] + , pattern = "_01" + , info = "leaves the 1st layer in 'x' unchanged" +) + +expect_match( + names(lagged)[[2L]][1L] + , pattern = "_05" + , info = "shifts the original 'y' series by 4 lags" +) + +## `Raster*` input +lagged_rst = lagalize(as(sst, "Raster"), as(pcp, "Raster"), freq = 12L) + +expect_equal( + lagged_rst[[1L]] + , target = sst + , info = "returns the original 'x' series when 'lag' is `NULL`" +) + +expect_equal( + lagged_rst[[2L]] + , target = pcp + , info = "returns the original 'y' series when 'lag' is `NULL`" +) diff --git a/inst/tinytest/test_longtermMeans.R b/inst/tinytest/test_longtermMeans.R new file mode 100644 index 0000000..21de188 --- /dev/null +++ b/inst/tinytest/test_longtermMeans.R @@ -0,0 +1,27 @@ +pcp = terra::unwrap(australiaGPCP) + +pcp_mn = longtermMeans(pcp, cycle.window = 12L) + +expect_inherits( + pcp_mn + , class = "SpatRaster" + , info = "returns a `SpatRaster` object" +) + +expect_true( + terra::nlyr(pcp_mn) == 12L + , info = "returns the correct number of layers" +) + +expect_identical( + dim(pcp_mn)[1:2] + , target = dim(pcp)[1:2] # nrows, ncols + , info = "returns the correct number of rows and columns" +) + +pcp_mn1 = longtermMeans(as(pcp, "Raster"), cycle.window = 1L) + +expect_true( + terra::nlyr(pcp_mn1) == 1L + , info = "returns the correct number of layers for overall mean" +) diff --git a/inst/tinytest/test_predict.R b/inst/tinytest/test_predict.R new file mode 100644 index 0000000..875707d --- /dev/null +++ b/inst/tinytest/test_predict.R @@ -0,0 +1,70 @@ +sst = terra::unwrap(pacificSST) +pcp = terra::unwrap(australiaGPCP) + +trn_idx = 1:10 +tst_idx = 11:20 +n = 3L + +trn = eot( + x = sst[[trn_idx]] + , y = pcp[[trn_idx]] + , n = n + , verbose = FALSE +) + +## `EotStack` input +prd_st = predict( + trn + , newdata = sst[[tst_idx]] + , n = n +) + +expect_inherits( + prd_st + , class = "SpatRaster" + , info = "returns a `SpatRaster` (`EotStack` input)" +) + +expect_identical( + dim(prd_st) + , target = dim(pcp[[tst_idx]]) + , info = "returns same dimensions as ground truth data (`EotStack` input)" +) + +## `EotMode` input +prd_md = predict( + trn[[1L]] + , newdata = sst[[tst_idx]] +) + +expect_inherits( + prd_md + , class = "SpatRaster" + , info = "returns a `SpatRaster` (`EotMode` input)" +) + +expect_identical( + dim(prd_md) + , target = dim(pcp[[tst_idx]]) + , info = "returns same dimensions as ground truth data (`EotMode` input)" +) + +prd_n1 = predict( + trn + , newdata = sst[[tst_idx]] + , n = 1L +) + +expect_equal( + prd_n1 + , target = prd_md + , info = "returns the same result with `EotStack` input and `n = 1L`" +) + +expect_true( + terra::global( + prd_md[[1L]] != prd_st[[1L]] + , fun = "min" + ) == 1L + , info = "returns different results for `n = 3L` vs. `n = 1L`" +) diff --git a/man/EotCycle.Rd b/man/EotCycle.Rd index ac2e6c5..069b137 100644 --- a/man/EotCycle.Rd +++ b/man/EotCycle.Rd @@ -2,9 +2,13 @@ % Please edit documentation in R/EotCycle.R \name{EotCycle} \alias{EotCycle} +\alias{EotCycle,RasterStackBrick-method} +\alias{EotCycle,SpatRaster-method} \title{Calculate a single EOT} \usage{ -EotCycle( +\S4method{EotCycle}{RasterStackBrick}(x, ...) + +\S4method{EotCycle}{SpatRaster}( x, y, n = 1, @@ -13,43 +17,19 @@ EotCycle( write.out, path.out, prefix, - type, + type = c("rsq", "ioa"), verbose, ... ) } \arguments{ -\item{x}{a ratser stack used as predictor} - -\item{y}{a RasterStack used as response. If \code{y} is \code{NULL}, -\code{x} is used as \code{y}} - -\item{n}{the number of EOT modes to calculate} +\item{x, y, n, standardised, write.out, path.out, prefix, type, verbose}{See \code{\link[=eot]{eot()}}.} -\item{standardised}{logical. If \code{FALSE} the calculated r-squared values -will be multiplied by the variance} +\item{...}{Further arguments passed to \code{\link[=writeEot]{writeEot()}} if \code{write.out = TRUE}, +or to the underlying \code{SpatRaster} method in general for \verb{Raster*} input.} \item{orig.var}{original variance of the response domain} - -\item{write.out}{logical. If \code{TRUE} results will be written to disk -using \code{path.out}} - -\item{path.out}{the file path for writing results if \code{write.out} is \code{TRUE}. -Defaults to current working directory} - -\item{prefix}{optional prefix to be used for naming of results if -\code{write.out} is \code{TRUE}} - -\item{type}{the type of the link function. Defaults to \code{'rsq'} as in original -proposed method from \cite{Dool2000}. If set to \code{'ioa'} index of agreement is -used instead} - -\item{verbose}{logical. If \code{TRUE} some details about the -calculation process will be output to the console} - -\item{...}{If \code{write.out = TRUE}, further arguments passed to -\code{\link[=writeEot]{writeEot()}}.} } \description{ -EotCycle() calculates a single EOT and is controlled by the main eot() function +Calculates a single EOT and is controlled by the main \code{\link[=eot]{eot()}} function. } diff --git a/man/EotMode-class.Rd b/man/EotMode-class.Rd index b32475f..a022c0f 100644 --- a/man/EotMode-class.Rd +++ b/man/EotMode-class.Rd @@ -22,30 +22,30 @@ Class EotMode \item{\code{cum_exp_var}}{the cumulative explained variance of the considered EOT mode} -\item{\code{r_predictor}}{the RasterLayer of the correlation coefficients between the base point and each pixel of the predictor domain} +\item{\code{r_predictor}}{\code{SpatRaster} of the correlation coefficients between the base point and each pixel of the predictor domain} \item{\code{rsq_predictor}}{as above but for the coefficient of determination of the predictor domain} \item{\code{rsq_sums_predictor}}{as above but for the sums of coefficient of determination of the predictor domain} -\item{\code{int_predictor}}{the RasterLayer of the intercept of the regression equation for each pixel of the predictor domain} +\item{\code{int_predictor}}{\code{SpatRaster} of the intercept of the regression equation for each pixel of the predictor domain} -\item{\code{slp_predictor}}{same as above but for the slope of the regression equation for each pixel of the predictor domain} +\item{\code{slp_predictor}}{as above but for the slope of the regression equation for each pixel of the predictor domain} -\item{\code{p_predictor}}{the RasterLayer of the significance (p-value) of the the regression equation for each pixel of the predictor domain} +\item{\code{p_predictor}}{\code{SpatRaster} of the significance (p-value) of the regression equation for each pixel of the predictor domain} -\item{\code{resid_predictor}}{the RasterBrick of the reduced data for the predictor domain} +\item{\code{resid_predictor}}{\code{SpatRaster} of the reduced data for the predictor domain} -\item{\code{r_response}}{the RasterLayer of the correlation coefficients between the base point and each pixel of the response domain} +\item{\code{r_response}}{\code{SpatRaster} of the correlation coefficients between the base point and each pixel of the response domain} \item{\code{rsq_response}}{as above but for the coefficient of determination of the response domain} -\item{\code{int_response}}{the RasterLayer of the intercept of the regression equation for each pixel of the response domain} +\item{\code{int_response}}{\code{SpatRaster} of the intercept of the regression equation for each pixel of the response domain} \item{\code{slp_response}}{as above but for the slope of the regression equation for each pixel of the response domain} -\item{\code{p_response}}{same the RasterLayer of the significance (p-value) of the the regression equation for each pixel of the response domain} +\item{\code{p_response}}{\code{SpatRaster} of the significance (p-value) of the regression equation for each pixel of the response domain} -\item{\code{resid_response}}{the RasterBrick of the reduced data for the response domain} +\item{\code{resid_response}}{\code{SpatRaster} of the reduced data for the response domain} }} diff --git a/man/anomalize.Rd b/man/anomalize.Rd index d0fbfc9..01b9b81 100644 --- a/man/anomalize.Rd +++ b/man/anomalize.Rd @@ -2,36 +2,42 @@ % Please edit documentation in R/anomalize.R \name{anomalize} \alias{anomalize} -\title{Create an anomaly RasterStack} +\alias{anomalize,RasterStackBrick-method} +\alias{anomalize,SpatRaster-method} +\title{Create an anomaly raster series} \usage{ -anomalize(x, reference = NULL, ...) +\S4method{anomalize}{RasterStackBrick}(x, ...) + +\S4method{anomalize}{SpatRaster}(x, reference = NULL, ...) } \arguments{ -\item{x}{a RasterStack} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} -\item{reference}{an optional RasterLayer to be used as the reference} +\item{...}{Additional arguments passed to \code{\link[terra:app]{terra::app()}} (e.g. 'cores', +'filename') to calculate the overall mean if 'reference' is \code{NULL}, or to +the underlying \code{SpatRaster} method for \verb{Raster*} input.} -\item{...}{additional arguments passed to \code{\link[raster:calc]{raster::calc()}} (and, in turn, -\code{\link[raster:writeRaster]{raster::writeRaster()}}) which is used under the hood} +\item{reference}{An optional single-layer \code{SpatRaster} (or \code{RasterLayer}) to +be used as the reference. Uses the overall mean of the original series if +\code{NULL} (default).} } \value{ -an anomaly RasterStack +An anomaly \code{SpatRaster} series. } \description{ -The function creates an anomaly RasterStack either based on the -overall mean of the original stack, or a supplied reference RasterLayer. +The function creates an anomaly raster series either based on +the overall mean of the original series, or a supplied reference raster. For the creation of seasonal anomalies use \code{\link[=deseason]{deseason()}}. } \examples{ -data(australiaGPCP) - -aus_anom <- anomalize(australiaGPCP) +pcp = terra::unwrap(australiaGPCP) +pcp_anom = anomalize(pcp) -opar <- par(mfrow = c(1,2)) -plot(australiaGPCP[[10]], main = "original") -plot(aus_anom[[10]], main = "anomalized") +opar = par(mfrow = c(1,2)) +plot(pcp[[10]], main = "original") +plot(pcp_anom[[10]], main = "anomalized") par(opar) } \seealso{ -\code{\link[=deseason]{deseason()}}, \code{\link[=denoise]{denoise()}}, \code{\link[raster:calc]{raster::calc()}} +\code{\link[=deseason]{deseason()}}, \code{\link[=denoise]{denoise()}} } diff --git a/man/australiaGPCP.Rd b/man/australiaGPCP.Rd index 17ae9bd..f5617d2 100644 --- a/man/australiaGPCP.Rd +++ b/man/australiaGPCP.Rd @@ -5,12 +5,16 @@ \alias{australiaGPCP} \title{Monthly GPCP precipitation data for Australia} \format{ -a RasterBrick with the following attributes\cr -\cr -dimensions : 12, 20, 240, 348 (nrow, ncol, ncell, nlayers)\cr -resolution : 2.5, 2.5 (x, y)\cr -extent : 110, 160, -40, -10 (xmin, xmax, ymin, ymax)\cr -coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs\cr +A \code{PackedSpatRaster} with the following attributes:\cr + +\if{html}{\out{
}}\preformatted{size : 12, 20, 348 (nrow, ncol, nlyr) +resolution : 2.5, 2.5 (x, y) +extent : 110, 160, -40, -10 (xmin, xmax, ymin, ymax) +coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs +}\if{html}{\out{
}} +} +\usage{ +australiaGPCP } \description{ Monthly Gridded Precipitation Climatology Project precipitation data @@ -20,6 +24,9 @@ for Australia from 1982/01 to 2010/12 Monthly Gridded Precipitation Climatology Project precipitation data for Australia from 1982/01 to 2010/12 } +\examples{ +terra::unwrap(australiaGPCP) +} \references{ The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979 - Present)\cr Adler et al. (2003)\cr diff --git a/man/calcVar.Rd b/man/calcVar.Rd index 31980ba..c31deb0 100644 --- a/man/calcVar.Rd +++ b/man/calcVar.Rd @@ -2,26 +2,32 @@ % Please edit documentation in R/calcVar.R \name{calcVar} \alias{calcVar} -\title{Calculate space-time variance of a RasterStack or RasterBrick} +\alias{calcVar,RasterStackBrick-method} +\alias{calcVar,SpatRaster-method} +\title{Calculate space-time variance of a raster series} \usage{ -calcVar(x, standardised = FALSE, ...) +\S4method{calcVar}{RasterStackBrick}(x, ...) + +\S4method{calcVar}{SpatRaster}(x, standardised = FALSE, ...) } \arguments{ -\item{x}{a RasterStack or RasterBrick} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} -\item{standardised}{logical.} +\item{...}{For \code{SpatRaster} input: currently not used. For \verb{Raster*} input: +arguments passed to the underlying \code{SpatRaster} method.} -\item{...}{currently not used} +\item{standardised}{\code{logical}, defaults to \code{FALSE}.} } \value{ -the mean (optionally standardised) space-time variance. +The mean (optionally standardised) space-time variance as \code{numeric}. } \description{ The function calculates the (optionally standardised) space-time -variance of a RasterStack or RasterBrick. +variance of a raster series. } \examples{ -data("pacificSST") +sst = terra::unwrap(pacificSST) -calcVar(pacificSST) +calcVar(sst) # default non-standardised +calcVar(sst, standardised = TRUE) } diff --git a/man/covWeight.Rd b/man/covWeight.Rd index 470512f..2de7706 100644 --- a/man/covWeight.Rd +++ b/man/covWeight.Rd @@ -7,19 +7,16 @@ covWeight(m, weights, ...) } \arguments{ -\item{m}{a matrix (e.g. as returned by \code{\link[raster:getValues]{raster::getValues()}})} +\item{m}{A \code{matrix}, e.g. as returned by \code{\link[terra:values]{terra::values()}}.} -\item{weights}{a numeric vector of weights. For lat/lon data this -can be produced with \code{\link[=getWeights]{getWeights()}}} +\item{weights}{A \code{numeric} vector of weights. For lat/lon data this can be +produced with \code{\link[=getWeights]{getWeights()}}.} -\item{...}{additional arguments passed to \code{\link[stats:cov.wt]{stats::cov.wt()}}} +\item{...}{Additional arguments passed to \code{\link[stats:cov.wt]{stats::cov.wt()}}} } \value{ -see \code{\link[stats:cov.wt]{stats::cov.wt()}} +See \code{\link[stats:cov.wt]{stats::cov.wt()}}. } \description{ Create a weighted covariance matrix } -\seealso{ -\code{\link[stats:cov.wt]{stats::cov.wt()}} -} diff --git a/man/cutStack.Rd b/man/cutStack.Rd index e523b6a..c34ef6d 100644 --- a/man/cutStack.Rd +++ b/man/cutStack.Rd @@ -2,33 +2,40 @@ % Please edit documentation in R/cutStack.R \name{cutStack} \alias{cutStack} -\title{Shorten a RasterStack} +\alias{cutStack,RasterStackBrick-method} +\alias{cutStack,SpatRaster-method} +\title{Shorten a raster series} \usage{ -cutStack(x, tail = TRUE, n = NULL) +\S4method{cutStack}{RasterStackBrick}(x, ...) + +\S4method{cutStack}{SpatRaster}(x, tail = TRUE, n = NULL) } \arguments{ -\item{x}{a RasterStack} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} + +\item{...}{Additional arguments passed to the underlying \code{SpatRaster} method.} -\item{tail}{logical. If \code{TRUE} the layers will be taken off +\item{tail}{\code{logical}. If \code{TRUE} (default) the layers will be taken off the end of the stack. If \code{FALSE} layers will be taken off the beginning.} -\item{n}{the number of layers to take away.} +\item{n}{The number of layers to take away as \code{integer}. If \code{NULL} (default), +'x' is returned unchanged.} } \value{ -a RasterStack shortened by \code{n} layers either from the -beginning or the end, depending on the specification of \code{tail} +A \code{SpatRaster} series shortened by 'n' layers either from the +beginning or the end, depending on the specification of 'tail'. } \description{ -The function cuts a specified number of layers off a RrasterStack in -order to create lagged RasterStacks. +The function cuts a specified number of layers off a raster +series in order to create a lagged stack. } \examples{ -data(australiaGPCP) +pcp = terra::unwrap(australiaGPCP) # 6 layers from the beginning -cutStack(australiaGPCP, tail = FALSE, n = 6) +cutStack(pcp, tail = FALSE, n = 6) # 8 layers from the end -cutStack(australiaGPCP, tail = TRUE, n = 8) +cutStack(pcp, tail = TRUE, n = 8) } diff --git a/man/deg2rad.Rd b/man/deg2rad.Rd index d69d0e4..78cb920 100644 --- a/man/deg2rad.Rd +++ b/man/deg2rad.Rd @@ -7,20 +7,34 @@ deg2rad(deg) } \arguments{ -\item{deg}{vector of degrees to be converted to radians} +\item{deg}{A series of degrees to be converted to radians.} +} +\value{ +A vector of radians. } \description{ -Convert degrees to radians +Converts a series of degrees to radians. } \examples{ -data(vdendool) - ## latitude in degrees -degrees <- coordinates(vdendool)[, 2] +gph = terra::unwrap(vdendool) + +degrees = terra::crds(gph)[, 2] head(degrees) ## latitude in radians -radians <- deg2rad(coordinates(vdendool)[, 2]) +radians = deg2rad(degrees) head(radians) +## `SpatRaster` input, e.g. useful for topographic operations +tmp = gph[[1L]] +terra::values(tmp) = degrees + +tmp_rad = deg2rad(tmp) + +opar = par(mfrow = c(1, 2)) +plot(tmp, main = "lat (degrees)") +plot(tmp_rad, main = "lat (radians)") +par(opar) + } diff --git a/man/denoise.Rd b/man/denoise.Rd index e9e5027..3968760 100644 --- a/man/denoise.Rd +++ b/man/denoise.Rd @@ -2,9 +2,13 @@ % Please edit documentation in R/denoise.R \name{denoise} \alias{denoise} +\alias{denoise,RasterStackBrick-method} +\alias{denoise,SpatRaster-method} \title{Noise filtering through principal components} \usage{ -denoise( +\S4method{denoise}{RasterStackBrick}(x, ...) + +\S4method{denoise}{SpatRaster}( x, k = NULL, expl.var = NULL, @@ -15,40 +19,44 @@ denoise( ) } \arguments{ -\item{x}{RasterStack to be filtered} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series to be filtered.} + +\item{...}{Additional arguments passed to \code{\link[stats:princomp]{stats::princomp()}}.} -\item{k}{number of components to be kept for reconstruction -(ignored if \code{expl.var} is supplied)} +\item{k}{The number of components to be kept for reconstruction (ignored if +'expl.var' is supplied).} -\item{expl.var}{minimum amount of variance to be kept after reconstruction -(should be set to NULL or omitted if \code{k} is supplied)} +\item{expl.var}{Minimum amount of variance to be kept after reconstruction +(should be set to \code{NULL} or omitted if 'k' is supplied).} \item{weighted}{logical. If \code{TRUE} the covariance matrix will be geographically weighted using the cosine of latitude during decomposition -(only important for lat/lon data)} +(only important for lat/lon data).} -\item{use.cpp}{logical. Determines whether to use \strong{Rcpp} -functionality, defaults to \code{TRUE}.} +\item{use.cpp}{logical. Determines whether to use \strong{Rcpp} functionality, +defaults to \code{TRUE}.} -\item{verbose}{logical. If \code{TRUE} some details about the -calculation process will be output to the console} - -\item{...}{additional arguments passed to \code{\link[stats:princomp]{stats::princomp()}}} +\item{verbose}{logical. If \code{TRUE} some details about the calculation process +will be output to the console.} } \value{ -a denoised RasterStack +A denoised \code{SpatRaster} series. } \description{ -Filter noise from a RasterStack by decomposing into principal components -and subsequent reconstruction using only a subset of components +Filter noise from a raster series by decomposing into principal components +and subsequent reconstruction using only a subset of components. +} +\note{ +Either 'k' or 'expl.var' must be specified. If both are supplied, 'k' will be +ignored. If none are supplied, an error will be thrown. } \examples{ -data("vdendool") -vdd_dns <- denoise(vdendool, expl.var = 0.8) +gph = terra::unwrap(vdendool) +gph_dns = denoise(gph, expl.var = 0.8) -opar <- par(mfrow = c(1,2)) -plot(vdendool[[1]], main = "original") -plot(vdd_dns[[1]], main = "denoised") +opar = par(mfrow = c(1,2)) +plot(gph[[1]], main = "original") +plot(gph_dns[[1]], main = "denoised") par(opar) } \seealso{ diff --git a/man/deseason.Rd b/man/deseason.Rd index 17174b4..a685b95 100644 --- a/man/deseason.Rd +++ b/man/deseason.Rd @@ -3,46 +3,49 @@ \name{deseason} \alias{deseason} \alias{deseason,RasterStackBrick-method} +\alias{deseason,SpatRaster-method} \alias{deseason,numeric-method} \title{Create seasonal anomalies} \usage{ -\S4method{deseason}{RasterStackBrick}(x, cycle.window = 12L, use.cpp = FALSE, filename = "", ...) +\S4method{deseason}{RasterStackBrick}(x, ...) + +\S4method{deseason}{SpatRaster}(x, cycle.window = 12L, use.cpp = FALSE, filename = "", ...) \S4method{deseason}{numeric}(x, cycle.window = 12L) } \arguments{ -\item{x}{An \verb{Raster*} object or, alternatively, a \code{numeric} time -series.} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) object or, alternatively, a \code{numeric} +time series.} -\item{cycle.window}{\code{integer}, defaults to \code{12}. The window for the -creation of the anomalies.} +\item{...}{For \code{SpatRaster} input: additional arguments passed to +\code{\link[terra:writeRaster]{terra::writeRaster()}}, only considered if \code{filename} is specified. +For \verb{Raster*} input: arguments passed to the underlying \code{SpatRaster} +method (e.g. \code{cycle.window}, \code{use.cpp}, \code{filename}).} -\item{use.cpp}{\code{logical}, defaults to \code{FALSE}. Determines whether -or not to use \strong{Rcpp} functionality. Only applies if \code{x} is a -\verb{Raster*} object.} +\item{cycle.window}{\code{integer}, defaults to \code{12}. The window for the creation +of the anomalies.} -\item{filename}{\code{character}. Output filename (optional).} +\item{use.cpp}{\code{logical}, defaults to \code{FALSE}. Determines whether or not to +use \strong{Rcpp} functionality. Only applies if \code{x} is a raster object.} -\item{...}{Additional arguments passed on to \code{\link[raster:writeRaster]{raster::writeRaster()}}, only -considered if \code{filename} is specified.} +\item{filename}{\code{character}. Output filename (optional).} } \value{ -If \code{x} is a \verb{Raster*} object, a deseasoned -\code{RasterStack}; else a deseasoned \code{numeric} vector. +If \code{x} is a raster object, a deseasoned \code{SpatRaster}; else a +deseasoned \code{numeric} vector. } \description{ -The function calculates anomalies of a RasterStack by supplying a -suitable seasonal window. E. g. to create monthly anomalies of a +The function calculates anomalies of a raster object by +supplying a suitable seasonal window. E.g. to create monthly anomalies of a raster stack of 12 layers per year, use \code{cycle.window = 12}. } \examples{ -data("australiaGPCP") - -aus_dsn <- deseason(australiaGPCP, 12) +pcp = terra::unwrap(australiaGPCP) +pcp_dsn = deseason(pcp, cycle.window = 12) -opar <- par(mfrow = c(1,2)) -plot(australiaGPCP[[1]], main = "original") -plot(aus_dsn[[1]], main = "deseasoned") +opar = par(mfrow = c(1, 2)) +plot(pcp[[1]], main = "original") +plot(pcp_dsn[[1]], main = "deseasoned") par(opar) } \seealso{ diff --git a/man/eot.Rd b/man/eot.Rd index 7df66db..9962306 100644 --- a/man/eot.Rd +++ b/man/eot.Rd @@ -3,9 +3,12 @@ \name{eot} \alias{eot} \alias{eot,RasterStackBrick-method} -\title{EOT analysis of a predictor and (optionally) a response RasterStack} +\alias{eot,SpatRaster-method} +\title{EOT analysis of a predictor and (optionally) a response raster series} \usage{ -\S4method{eot}{RasterStackBrick}( +\S4method{eot}{RasterStackBrick}(x, ...) + +\S4method{eot}{SpatRaster}( x, y = NULL, n = 1, @@ -20,40 +23,40 @@ ) } \arguments{ -\item{x}{a \verb{Raster*} object used as predictor} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) object used as predictor.} -\item{y}{a \verb{Raster*} object used as response. If \code{y} is -\code{NULL}, \code{x} is used as \code{y}} +\item{...}{Currently not used.} -\item{n}{the number of EOT modes to calculate} +\item{y}{A \code{SpatRaster} (or \verb{Raster*}) object used as response. If 'y' is +\code{NULL} (default), 'x' is used as response.} -\item{standardised}{logical. If \code{FALSE} the calculated r-squared values -will be multiplied by the variance} +\item{n}{The number of EOT modes to calculate as \code{integer}, defaults to \code{1}.} -\item{write.out}{logical. If \code{TRUE} results will be written to disk -using \code{path.out}} +\item{standardised}{\code{logical}, default \code{TRUE}. If \code{FALSE} the calculated +r-squared values will be multiplied by the variance.} -\item{path.out}{the file path for writing results if \code{write.out} is \code{TRUE}. -Defaults to current working directory} +\item{write.out}{\code{logical}, default \code{FALSE}. If \code{TRUE} results will be +written to disk using 'path.out'.} -\item{prefix}{optional prefix to be used for naming of results if -\code{write.out} is \code{TRUE}} +\item{path.out}{The file path for writing results if 'write.out' is \code{TRUE}. +Defaults to current working directory.} -\item{reduce.both}{logical. If \code{TRUE} both \code{x} and \code{y} -are reduced after each iteration. If \code{FALSE} only \code{y} is reduced} +\item{prefix}{optional prefix to be used for naming of results if 'write.out' +is \code{TRUE}.} -\item{type}{the type of the link function. Defaults to \code{'rsq'} as in original -proposed method from \cite{van den Dool 2000}. If set to \code{'ioa'} index of agreement is -used instead} +\item{reduce.both}{\code{logical}, default \code{FALSE}. If \code{TRUE} both 'x' and 'y' are +reduced after each iteration. If \code{FALSE} only 'y' is reduced.} -\item{verbose}{logical. If \code{TRUE} some details about the -calculation process will be output to the console} +\item{type}{The type of the link function. Defaults to \code{"rsq"} as in original +proposed method from \cite{van den Dool 2000}. If set to \code{"ioa"} index of +agreement is used instead.} -\item{...}{not used at the moment} +\item{verbose}{\code{logical}, default \code{FALSE}. If \code{TRUE} some details about the +calculation process will be output to the console.} } \value{ -if n = 1 an \emph{EotMode}, if n > 1 an \emph{EotStack} of \code{n} -\emph{EotMode}s. Each \emph{EotMode} has the following components: +If \verb{n = 1`` an }EotMode\verb{, if n > 1 an }EotStack\verb{of 'n'}EotMode\verb{s. Each }EotMode` has the following components: + \itemize{ \item \emph{mode} - the number of the identified mode (1 - n) \item \emph{eot} - the EOT (time series) at the identified base point. @@ -61,17 +64,17 @@ Note, this is a simple numeric vector, not of class \code{ts} \item \emph{coords_bp} - the coordinates of the identified base point \item \emph{cell_bp} - the cell number of the indeified base point \item \emph{cum_exp_var} - the (cumulative) explained variance of the considered EOT -\item \emph{r_predictor} - the \emph{RasterLayer} of the correlation coefficients +\item \emph{r_predictor} - the \code{SpatRaster} of the correlation coefficients between the base point and each pixel of the predictor domain \item \emph{rsq_predictor} - as above but for the coefficient of determination \item \emph{rsq_sums_predictor} - as above but for the sums of coefficient of determination -\item \emph{int_predictor} - the \emph{RasterLayer} of the intercept of the +\item \emph{int_predictor} - the \code{SpatRaster} of the intercept of the regression equation for each pixel of the predictor domain \item \emph{slp_predictor} - same as above but for the slope of the regression equation for each pixel of the predictor domain -\item \emph{p_predictor} - the \emph{RasterLayer} of the significance (p-value) +\item \emph{p_predictor} - the \code{SpatRaster} of the significance (p-value) of the the regression equation for each pixel of the predictor domain -\item \emph{resid_predictor} - the \emph{RasterBrick} of the reduced data +\item \emph{resid_predictor} - the \code{SpatRaster} of the reduced data for the predictor domain } @@ -79,11 +82,11 @@ Apart from \emph{rsq_sums_predictor}, all \emph{*_predictor} fields are also returned for the \emph{*_response} domain, even if predictor and response domain are equal. This is due to that fact, that if not both fields are reduced after the first EOT is found, -these \emph{RasterLayers} will differ. +these \code{SpatRaster}s will differ. } \description{ -Calculate a given number of EOT modes either internally or between -RasterStacks. +Calculate a given number of EOT modes either internally or between raster +series. } \details{ For a detailed description of the EOT algorithm and the mathematics behind it, @@ -106,16 +109,17 @@ variability within the same or another geographic field. ### EXAMPLE I ### a single field \donttest{ -data(vdendool) +gph = terra::unwrap(vdendool) -## claculate 2 leading modes -nh_modes <- eot(x = vdendool, y = NULL, n = 2, +## calculate 2 leading modes +nh_modes <- eot(x = gph, y = NULL, n = 2, standardised = FALSE, verbose = TRUE) plot(nh_modes, y = 1, show.bp = TRUE) plot(nh_modes, y = 2, show.bp = TRUE) } + } \references{ \bold{Empirical Orthogonal Teleconnections}\cr diff --git a/man/geoWeight.Rd b/man/geoWeight.Rd index d78d9cc..2ac7763 100644 --- a/man/geoWeight.Rd +++ b/man/geoWeight.Rd @@ -2,33 +2,32 @@ % Please edit documentation in R/geoWeight.R \name{geoWeight} \alias{geoWeight} +\alias{geoWeight,RasterStackBrick-method} +\alias{geoWeight,SpatRaster-method} \title{Geographic weighting} \usage{ -geoWeight(x, f = function(x) cos(x), ...) +\S4method{geoWeight}{RasterStackBrick}(x, ...) + +\S4method{geoWeight}{SpatRaster}(x, f = cos, ...) } \arguments{ -\item{x}{a Raster* object} - -\item{f}{a function to be used to the weighting. -Defaults to \code{cos(x)}} - -\item{...}{additional arguments to be passed to f} +\item{x, f, ...}{See \code{\link[=getWeights]{getWeights()}}. In case of \verb{Raster*} input, '...' funnels +additional arguments (i.e. 'f') to the underlying \code{SpatRaster} method.} } \value{ -a weighted Raster* object +A weighted \code{SpatRaster}. } \description{ -The function performs geographic weighting of non-projected long/lat -data. By default it uses the cosine of latitude to compensate for the -area distortion, though the user can supply other functions via \code{f}. +The function performs geographic weighting of non-projected long/lat data. By +default it uses the cosine of latitude (in radians) to compensate for the +area distortion, though the user can supply other weighting functions. } \examples{ -data(vdendool) - -wgtd <- geoWeight(vdendool) +gph <- terra::unwrap(vdendool) +wgtd <- geoWeight(gph) opar <- par(mfrow = c(1,2)) -plot(vdendool[[1]], main = "original") +plot(gph[[1]], main = "original") plot(wgtd[[1]], main = "weighted") par(opar) } diff --git a/man/getWeights.Rd b/man/getWeights.Rd index 40f9a79..da67772 100644 --- a/man/getWeights.Rd +++ b/man/getWeights.Rd @@ -2,33 +2,40 @@ % Please edit documentation in R/getWeights.R \name{getWeights} \alias{getWeights} +\alias{getWeights,RasterStackBrick-method} +\alias{getWeights,SpatRaster-method} \title{Calculate weights from latitude} \usage{ -getWeights(x, f = function(x) cos(x), ...) +\S4method{getWeights}{RasterStackBrick}(x, ...) + +\S4method{getWeights}{SpatRaster}(x, f = cos, ...) } \arguments{ -\item{x}{a Raster* object} +\item{x}{A non-projected \code{SpatRaster} (or \verb{Raster*}) object.} -\item{f}{a function to be used to the weighting. -Defaults to \code{cos(x)}} +\item{...}{Additional arguments passed to 'f', or to the underlying +\code{SpatRaster} method in general for \verb{Raster*} input.} -\item{...}{additional arguments to be passed to f} +\item{f}{A \code{function} applied to the latitude (in radians) to compute +weights. Defaults to \code{cos}.} } \value{ -a numeric vector of weights +A \code{numeric} vector of weights for non-\code{NA} cells in 'x'. } \description{ -Calculate weights using the cosine of latitude to compensate for area -distortion of non-projected lat/lon data +Calculate weights using the cosine of latitude to compensate for +area distortion of non-projected lat/lon data. } \examples{ -data("australiaGPCP") -wghts <- getWeights(australiaGPCP) -wghts_rst <- australiaGPCP[[1]] -wghts_rst[] <- wghts +pcp = terra::unwrap(australiaGPCP) + +wghts = getWeights(pcp) +utils::head(wghts) + +wghts_rst = terra::setValues(pcp[[1]], wghts) -opar <- par(mfrow = c(1,2)) -plot(australiaGPCP[[1]], main = "data") +opar = par(mfrow = c(1,2)) +plot(pcp[[1]], main = "data") plot(wghts_rst, main = "weights") par(opar) diff --git a/man/lagalize.Rd b/man/lagalize.Rd index 72dc9c2..50fd22c 100644 --- a/man/lagalize.Rd +++ b/man/lagalize.Rd @@ -2,36 +2,42 @@ % Please edit documentation in R/lagalize.R \name{lagalize} \alias{lagalize} -\title{Create lagged RasterStacks} +\alias{lagalize,RasterStackBrick,RasterStackBrick-method} +\alias{lagalize,SpatRaster,SpatRaster-method} +\title{Create lagged raster series} \usage{ -lagalize(x, y, lag = NULL, freq = 12, ...) +\S4method{lagalize}{RasterStackBrick,RasterStackBrick}(x, y, ...) + +\S4method{lagalize}{SpatRaster,SpatRaster}(x, y, lag = NULL, freq = 12L, ...) } \arguments{ -\item{x}{a RasterStack (to be cut from tail)} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series to be cut from tail.} -\item{y}{a RasterStack (to be cut from beginning)} +\item{y}{A \code{SpatRaster} (or \verb{Raster*}) series to be cut from beginning.} -\item{lag}{the desired lag (in the native frequency of the RasterStack)} +\item{...}{For \code{SpatRaster} input: currently not used. For \verb{Raster*} input: +arguments passed to the underlying \code{SpatRaster} method.} -\item{freq}{the frequency of the RasterStacks} +\item{lag}{The desired lag in the native frequency of the series passed to +\code{\link[=cutStack]{cutStack()}}.} -\item{...}{currently not used} +\item{freq}{The frequency of the raster series as \code{integer}.} } \value{ -a list with the two RasterStacks lagged by \code{lag} +A \code{list} with the two raster series lagged by 'lag'. } \description{ -The function is used to produce two lagged RasterStacks. The second is cut +The function is used to produce two lagged raster series. The second is cut from the beginning, the first from the tail to ensure equal output lengths (provided that input lengths were equal). } \examples{ -data(pacificSST) -data(australiaGPCP) +sst = terra::unwrap(pacificSST) +pcp = terra::unwrap(australiaGPCP) # lag GPCP by 4 months -lagged <- lagalize(pacificSST, australiaGPCP, lag = 4, freq = 12) -lagged[[1]][[1]] #check names to see date of layer -lagged[[2]][[1]] #check names to see date of layer +lagged = lagalize(sst, pcp, lag = 4, freq = 12) +lagged[[1]][[1]] # check names to see date of layer +lagged[[2]][[1]] # -"- } diff --git a/man/longtermMeans.Rd b/man/longtermMeans.Rd index 20384e3..23882c8 100644 --- a/man/longtermMeans.Rd +++ b/man/longtermMeans.Rd @@ -2,29 +2,36 @@ % Please edit documentation in R/longtermMeans.R \name{longtermMeans} \alias{longtermMeans} -\title{Calculate long-term means from a 'RasterStack'} +\alias{longtermMeans,RasterStackBrick-method} +\alias{longtermMeans,SpatRaster-method} +\title{Calculate long-term means from a raster series} \usage{ -longtermMeans(x, cycle.window = 12L) +\S4method{longtermMeans}{RasterStackBrick}(x, ...) + +\S4method{longtermMeans}{SpatRaster}(x, cycle.window = 12L) } \arguments{ -\item{x}{A 'RasterStack' (or 'RasterBrick') object.} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} + +\item{...}{For \verb{Raster*} input, arguments passed to the underlying +\code{SpatRaster} method.} -\item{cycle.window}{'integer'. See \code{\link[=deseason]{deseason()}}.} +\item{cycle.window}{\code{integer}, defaults to \code{12}. See \code{\link[=deseason]{deseason()}}.} } \value{ -If \code{cycle.window} equals \code{nlayers(x)} (which obviously doesn't make -much sense), a 'RasterLayer' object; else a 'RasterStack' object. +A \code{SpatRaster} with 'cycle.window' layers, each containing the mean across +all corresponding time steps. } \description{ -Calculate long-term means from an input 'RasterStack' (or 'RasterBrick') -object. Ideally, the number of input layers should be divisable by the -supplied \code{cycle.window}. For instance, if \code{x} consists of monthly -layers, \code{cycle.window} should be a multiple of 12. +Calculate long-term means from an input raster series. Ideally, the number of +input layers should be divisable by the supplied 'cycle.window'. For +instance, if 'x' consists of monthly layers, 'cycle.window' should be a +multiple of \code{12}. } \examples{ -data("australiaGPCP") +pcp = terra::unwrap(australiaGPCP) -longtermMeans(australiaGPCP) +longtermMeans(pcp) } \seealso{ diff --git a/man/nXplain.Rd b/man/nXplain.Rd index 9a551dd..8e064c8 100644 --- a/man/nXplain.Rd +++ b/man/nXplain.Rd @@ -8,27 +8,28 @@ \S4method{nXplain}{EotStack}(x, var = 0.9) } \arguments{ -\item{x}{an \emph{EotStack}} +\item{x}{An \code{EotStack}.} -\item{var}{the minimum amount of variance to be explained by the modes} +\item{var}{The minimum amount of variance to be explained by the modes as +\code{numeric}, defaults to \code{0.9}.} } \value{ -an integer denoting the number of EOTs needed to explain \code{var} +The number of EOTs needed to explain 'var'. } \description{ -The function identifies the number of modes needed to explain a certain amount of -variance within the response field. +The function identifies the number of modes needed to explain a certain +amount of variance within the response field. } \note{ -This is a post-hoc function. It needs an \emph{EotStack} -created as returned by \code{\link[=eot]{eot()}}. Depending on the potency -of the identified EOTs, it may be necessary to compute a high number of -modes in order to be able to explain a large enough part of the variance. +This is a post-hoc function. It needs an \code{EotStack} created as returned +by \code{\link[=eot]{eot()}}. Depending on the potency of the identified EOTs, it may be +necessary to compute a high number of modes in order to be able to explain +a large enough part of the variance. } \examples{ -data(vdendool) +gph <- terra::unwrap(vdendool) -nh_modes <- eot(x = vdendool, y = NULL, n = 3, +nh_modes <- eot(x = gph, y = NULL, n = 3, standardised = FALSE, verbose = TRUE) diff --git a/man/names.Rd b/man/names.Rd index 7c1425c..280c761 100644 --- a/man/names.Rd +++ b/man/names.Rd @@ -30,9 +30,9 @@ if \code{x} is a EotMode, the name the respective mode Get or set names of Eot* objects } \examples{ -data(vdendool) +gph <- terra::unwrap(vdendool) -nh_modes <- eot(vdendool, n = 2) +nh_modes <- eot(gph, n = 2) ## mode names names(nh_modes) diff --git a/man/nmodes.Rd b/man/nmodes.Rd index 5ff6b79..f702f66 100644 --- a/man/nmodes.Rd +++ b/man/nmodes.Rd @@ -2,28 +2,22 @@ % Please edit documentation in R/nmodes.R \name{nmodes} \alias{nmodes} -\alias{nmodes,EotStack-method} -\title{Number of modes of an EotStack} +\title{Number of modes of an \code{EotStack}} \usage{ -\S4method{nmodes}{EotStack}(x) +nmodes(x) } \arguments{ -\item{x}{an EotStack} +\item{x}{An \code{EotStack}.} } \value{ -integer +The number of modes as \code{integer}. } \description{ -Number of modes of an EotStack -} -\details{ -retrieves the number of modes of an EotStack +Retrieves the number of modes of an \code{EotStack}. } \examples{ -data(vdendool) - -nh_modes <- eot(vdendool, n = 2) - +gph = terra::unwrap(vdendool) +nh_modes = eot(gph, n = 2) nmodes(nh_modes) } diff --git a/man/pacificSST.Rd b/man/pacificSST.Rd index 46d1afe..dc7438a 100644 --- a/man/pacificSST.Rd +++ b/man/pacificSST.Rd @@ -5,12 +5,16 @@ \alias{pacificSST} \title{Monthly SSTs for the tropical Pacific Ocean} \format{ -a RasterBrick with the following attributes\cr -\cr -dimensions : 30, 140, 4200, 348 (nrow, ncol, ncell, nlayers)\cr -resolution : 1, 1 (x, y)\cr -extent : 150, 290, -15, 15 (xmin, xmax, ymin, ymax)\cr -coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs\cr +A \code{PackedSpatRaster} with the following attributes:\cr + +\if{html}{\out{
}}\preformatted{size : 30, 140, 348 (nrow, ncol, nlyr) +resolution : 1, 1 (x, y) +extent : 150, 290, -15, 15 (xmin, xmax, ymin, ymax) +coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs +}\if{html}{\out{
}} +} +\usage{ +pacificSST } \description{ Monthly NOAA sea surface temperatures for the tropical Pacific Ocean from 1982/01 to 2010/12 @@ -18,6 +22,9 @@ Monthly NOAA sea surface temperatures for the tropical Pacific Ocean from 1982/0 \details{ Monthly NOAA sea surface temperatures for the tropical Pacific Ocean from 1982/01 to 2010/12 } +\examples{ +terra::unwrap(pacificSST) +} \references{ Daily High-Resolution-Blended Analyses for Sea Surface Temperature\cr Reynolds et al. (2007)\cr diff --git a/man/plot.Rd b/man/plot.Rd index 6582d81..3bbda6e 100644 --- a/man/plot.Rd +++ b/man/plot.Rd @@ -69,7 +69,7 @@ predictor and response fields} produce a map showing the locations of all modes. Ignored if x is an EotMode} -\item{...}{further arguments to be passed to \code{\link[raster:spplot]{raster::spplot()}}} +\item{...}{further arguments to be passed to \code{\link[sp:spplot]{sp::spplot()}}} } \description{ This is the standard plotting routine for the results of \code{\link[=eot]{eot()}}. @@ -82,10 +82,10 @@ domain, iii) the time series at the identified base point }} \examples{ -data(vdendool) +gph <- terra::unwrap(vdendool) ## claculate 2 leading modes -nh_modes <- eot(x = vdendool, y = NULL, n = 2, +nh_modes <- eot(x = gph, y = NULL, n = 2, standardised = FALSE, verbose = TRUE) diff --git a/man/predict.Rd b/man/predict.Rd index 1bed0a3..165b632 100644 --- a/man/predict.Rd +++ b/man/predict.Rd @@ -6,59 +6,57 @@ \alias{predict,EotMode-method} \title{EOT based spatial prediction} \usage{ -\S4method{predict}{EotStack}(object, newdata, n = 1, cores = 1L, filename = "", ...) +\S4method{predict}{EotStack}(object, newdata, n = 1L, filename = "", ...) -\S4method{predict}{EotMode}(object, newdata, n = 1, cores = 1L, filename = "", ...) +\S4method{predict}{EotMode}(object, newdata, n = 1L, filename = "", ...) } \arguments{ -\item{object}{an \code{Eot*} object} +\item{object}{An \verb{Eot*} object.} -\item{newdata}{the data to be used as predictor} +\item{newdata}{The data to be used as predictor.} -\item{n}{the number of modes to be used for the prediction. +\item{n}{The number of modes to be used for the prediction. See \code{\link[=nXplain]{nXplain()}} for calculating the number of modes based on their explanatory power.} -\item{cores}{\code{integer}. Number of cores for parallel processing.} - \item{filename}{\code{character}, output filenames (optional). If specified, this must be of the same length as \code{nlayers(newdata)}.} -\item{...}{further arguments passed to \code{\link[raster:calc]{raster::calc()}}, and hence, -\code{\link[raster:writeRaster]{raster::writeRaster()}}.} +\item{...}{Further arguments passed to \code{\link[terra:app]{terra::app()}}, and hence, +\code{\link[terra:writeRaster]{terra::writeRaster()}}.} } \value{ -a \emph{RasterStack} of \code{nlayers(newdata)} +A \code{SpatRaster} of \code{nlayers(newdata)}. } \description{ Make spatial predictions using the fitted model returned by -\code{\link[=eot]{eot()}}. A (user-defined) set of \emph{n} modes will be used to +\code{\link[=eot]{eot()}}. A (user-defined) set of 'n' modes will be used to model the outcome using the identified link functions of the respective modes which are added together to produce the final prediction. } \examples{ ### not very useful, but highlights the workflow \donttest{ -data(pacificSST) -data(australiaGPCP) +sst = terra::unwrap(pacificSST) +pcp = terra::unwrap(australiaGPCP) ## train data using eot() -train <- eot(x = pacificSST[[1:10]], - y = australiaGPCP[[1:10]], +train <- eot(x = sst[[1:10]], + y = pcp[[1:10]], n = 1) ## predict using identified model pred <- predict(train, - newdata = pacificSST[[11:20]], + newdata = sst[[11:20]], n = 1) ## compare results opar <- par(mfrow = c(1,2)) -plot(australiaGPCP[[13]], main = "original", zlim = c(0, 10)) +plot(pcp[[13]], main = "original", zlim = c(0, 10)) plot(pred[[3]], main = "predicted", zlim = c(0, 10)) par(opar) } } \seealso{ -\code{\link[raster:calc]{raster::calc()}}, \code{\link[raster:writeRaster]{raster::writeRaster()}}. +\code{\link[terra:app]{terra::app()}}, \code{\link[terra:writeRaster]{terra::writeRaster()}}. } diff --git a/man/readEot.Rd b/man/readEot.Rd index 5bc1b77..814495d 100644 --- a/man/readEot.Rd +++ b/man/readEot.Rd @@ -2,44 +2,44 @@ % Please edit documentation in R/readEot.R \name{readEot} \alias{readEot} -\title{Read \code{Eot}* files from disk} +\title{Read \verb{Eot*} files from disk} \usage{ -readEot(x, prefix = "remote", suffix = "grd") +readEot(x, prefix = "remote", suffix = ".grd") } \arguments{ -\item{x}{\code{character}, search path for \code{Eot}* related files passed +\item{x}{\code{character}, search path for \verb{Eot*} related files passed to \code{\link[=list.files]{list.files()}}.} -\item{prefix}{\code{character}, see \code{\link[=writeEot]{writeEot()}} for details. -Should be the same as previously supplied to \code{\link[=eot]{eot()}}.} +\item{prefix}{\code{character}, see \code{\link[=writeEot]{writeEot()}} for details. Should be the same +'prefix' as used during file creation in \code{\link[=eot]{eot()}} or \code{\link[=writeEot]{writeEot()}}.} -\item{suffix}{\code{character}, file extension depending on the output file -type of locally stored \code{Eot}* files, see \code{\link[raster:writeRaster]{raster::writeRaster()}}.} +\item{suffix}{\code{character}, default \code{.grd} for native \code{"RRASTER"} format. File +extension depending on the output file type of locally stored \verb{Eot*} files.} } \value{ -An \code{Eot}* object. +An \code{EotMode} if a single mode is found on disk, or an \code{EotStack} +otherwise. } \description{ -Read \code{Eot}* related files from disk, e.g. for further use with +Read \verb{Eot*} related files from disk, e.g. for further use with \code{\link[=predict]{predict()}} or \code{\link[=plot]{plot()}}. } \examples{ \dontrun{ ## calculate 3 leading modes -data(vdendool) -nh_modes <- eot(x = vdendool, n = 3, standardised = FALSE, - write.out = TRUE, path.out = "~/data") +gph <- terra::unwrap(vdendool) +nh_modes <- eot(x = gph, n = 3, standardised = FALSE, + write.out = TRUE, path.out = tempdir()) ## reimport related files rm(nh_modes) -nh_modes <- readEot("~/data") +nh_modes <- readEot(tempdir()) nh_modes } } \seealso{ -\code{\link[=eot]{eot()}}, \code{\link[=writeEot]{writeEot()}}, -\code{\link[raster:writeRaster]{raster::writeRaster()}}. +\code{\link[=writeEot]{writeEot()}} } \author{ Florian Detsch diff --git a/man/subset.Rd b/man/subset.Rd index 6101567..59b25c6 100644 --- a/man/subset.Rd +++ b/man/subset.Rd @@ -4,34 +4,36 @@ \alias{subset} \alias{subset,EotStack-method} \alias{[[,EotStack,ANY,ANY-method} -\title{Subset modes in EotStacks} +\title{Subset modes in an \code{EotStack}} \usage{ \S4method{subset}{EotStack}(x, subset, drop = FALSE, ...) \S4method{[[}{EotStack,ANY,ANY}(x, i) } \arguments{ -\item{x}{EotStack to be subset} +\item{x}{\code{EotStack} to be subset} -\item{subset}{integer or character. The modes to ectract (either by -integer or by their names)} +\item{subset}{\code{integer} or \code{character}. The modes to extract (either by +their indexes or names).} -\item{drop}{if \code{TRUE} a single mode will be returned as an EotMode} +\item{drop}{If \code{TRUE}, a single selected mode is returned as an \code{EotMode}. +Defaults to \code{FALSE}, which always returns an \code{EotStack}.} -\item{...}{currently not used} +\item{...}{Currently not used.} -\item{i}{number of EotMode to be subset} +\item{i}{Index(es) or name(s) to be subset, delegated to \code{\link[=subset]{subset()}}.} } \value{ -an Eot* object +An \code{EotMode} if a single mode is selected and \code{drop = TRUE}, otherwise an +\code{EotStack}. } \description{ -Extract a set of modes from an EotStack +Extract a set of modes from an \code{EotStack}. } \examples{ -data(vdendool) +gph <- terra::unwrap(vdendool) -nh_modes <- eot(x = vdendool, y = NULL, n = 3, +nh_modes <- eot(x = gph, y = NULL, n = 3, standardised = FALSE, verbose = TRUE) @@ -44,4 +46,8 @@ class(subs) subs <- subset(nh_modes, 2, drop = TRUE) class(subs) +## similarly to `drop = TRUE` above: +nh_modes[[2L]] +nh_modes[["mode_02"]] + } diff --git a/man/vdendool.Rd b/man/vdendool.Rd index c72e21b..60df976 100644 --- a/man/vdendool.Rd +++ b/man/vdendool.Rd @@ -5,23 +5,30 @@ \alias{vdendool} \title{Mean seasonal (DJF) 700 mb geopotential heights} \format{ -a RasterBrick with the following attributes\cr -\cr -dimensions : 14, 36, 504, 50 (nrow, ncol, ncell, nlayers)\cr -resolution : 10, 4.931507 (x, y)\cr -extent : -180, 180, 20.9589, 90 (xmin, xmax, ymin, ymax)\cr -coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0\cr +A \code{PackedSpatRaster} with the following attributes:\cr + +\if{html}{\out{
}}\preformatted{size : 14, 36, 50 (nrow, ncol, nlyr) +resolution : 10, 5 (x, y) +extent : -180, 180, 20, 90 (xmin, xmax, ymin, ymax) +coord. ref. : +proj=longlat +datum=WGS84 +no_defs +}\if{html}{\out{
}} } \source{ \url{https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.derived.pressure.html}\cr \emph{Original Source:} NOAA National Center for Environmental Prediction } +\usage{ +vdendool +} \description{ NCEP/NCAR reanalysis data of mean seasonal (DJF) 700 mb geopotential heights from 1948 to 1998 } \details{ NCEP/NCAR reanalysis data of mean seasonal (DJF) 700 mb geopotential heights from 1948 to 1998 } +\examples{ +terra::unwrap(vdendool) +} \references{ The NCEP/NCAR 40-year reanalysis project\cr Kalnay et al. (1996)\cr diff --git a/man/writeEot.Rd b/man/writeEot.Rd index 18d399d..4142512 100644 --- a/man/writeEot.Rd +++ b/man/writeEot.Rd @@ -4,44 +4,48 @@ \alias{writeEot} \alias{writeEot,EotMode-method} \alias{writeEot,EotStack-method} -\title{Write Eot* objects to disk} +\title{Write \verb{Eot*} objects to disk} \usage{ \S4method{writeEot}{EotMode}(x, path.out = ".", prefix = "remote", overwrite = TRUE, ...) -\S4method{writeEot}{EotStack}(x, path.out = ".", prefix, ...) +\S4method{writeEot}{EotStack}(x, ...) } \arguments{ -\item{x}{an Eot* object} +\item{x}{An \verb{Eot*} object.} -\item{path.out}{the path to the folder to write the files to} +\item{path.out}{The path to the folder to write the files to.} -\item{prefix}{a prefix to be added to the file names (see Details)} +\item{prefix}{A prefix to be added to the file names (see Details).} -\item{overwrite}{see \code{\link[raster:writeRaster]{raster::writeRaster()}}. -Defaults to \code{TRUE} in \code{\link[=writeEot]{writeEot()}}} +\item{overwrite}{See \code{\link[terra:writeRaster]{terra::writeRaster()}}. In \code{\link[=writeEot]{writeEot()}}, this defaults +to \code{TRUE}.} -\item{...}{further arguments passed to \code{\link[raster:writeRaster]{raster::writeRaster()}}} +\item{...}{Further arguments passed to \code{\link[terra:writeRaster]{terra::writeRaster()}} for \code{EotMode} +input, or to the underlying \code{EotMode} method for \code{EotStack} input. If +'filetype' is not specified, the default \code{"RRASTER"} format (\code{.grd}) is +used.} } \description{ -Write Eot* objects to disk. This is merely a wrapper around -\code{\link[raster:writeRaster]{raster::writeRaster()}} so see respective help section for details. +Write \verb{Eot*} objects to disk. This is merely a wrapper around +\code{\link[terra:writeRaster]{terra::writeRaster()}} so see respective help section for details. } \details{ -\code{\link[=writeEot]{writeEot()}} will write the results of either an EotMode or an EotStack +\code{\link[=writeEot]{writeEot()}} will write the results of either an \code{EotMode} or an \code{EotStack} to disk. For each mode the following files will be written: \itemize{ -\item \emph{pred_r} - the \emph{RasterLayer} of the correlation coefficients +\item \emph{pred_r} - the \code{SpatRaster} of the correlation coefficients between the base point and each pixel of the predictor domain \item \emph{pred_rsq} - as above but for the coefficient of determination -\item \emph{pred_rsq_sums} - as above but for the sums of coefficient of determination -\item \emph{pred_int} - the \emph{RasterLayer} of the intercept of the +\item \emph{pred_rsq_sums} - as above but for the sums of coefficient of +determination +\item \emph{pred_int} - the \code{SpatRaster} of the intercept of the regression equation for each pixel of the predictor domain \item \emph{pred_slp} - same as above but for the slope of the regression equation for each pixel of the predictor domain -\item \emph{pred_p} - the \emph{RasterLayer} of the significance (p-value) +\item \emph{pred_p} - the \code{SpatRaster} of the significance (p-value) of the the regression equation for each pixel of the predictor domain -\item \emph{pred_resid} - the \emph{RasterBrick} of the reduced data +\item \emph{pred_resid} - the \code{SpatRaster} of the reduced data for the predictor domain } @@ -52,19 +56,19 @@ look like, e.g.: \emph{prefix_mode_n_pred_r.grd} -for the \emph{RasterLayer} of the predictor correlation coefficient -of mode n using the standard \emph{raster} file type (.grd). +for the \code{SpatRaster} of the predictor correlation coefficient +of mode n using the standard \pkg{raster} file type (\code{.grd}). } -\section{Methods (by class)}{ +\section{Functions}{ \itemize{ \item \code{writeEot(EotStack)}: EotStack }} \examples{ \dontrun{ -data(vdendool) +gph <- terra::unwrap(vdendool) -nh_modes <- eot(x = vdendool, y = NULL, n = 2, +nh_modes <- eot(x = gph, y = NULL, n = 2, standardised = FALSE, verbose = TRUE) @@ -76,5 +80,5 @@ writeEot(nh_modes[[2]], prefix = "vdendool") } } \seealso{ -\code{\link[raster:writeRaster]{raster::writeRaster()}} +\code{\link[terra:writeRaster]{terra::writeRaster()}} } diff --git a/tests/tinytest.R b/tests/tinytest.R new file mode 100644 index 0000000..cca18f5 --- /dev/null +++ b/tests/tinytest.R @@ -0,0 +1,5 @@ + +if (requireNamespace("tinytest", quietly = TRUE)) { + tinytest::test_package("remote") +} +