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..e6558db 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.9010 +Authors@R: c(person(given = "Tim", family = "Appelhans", role = c("cre", "aut"), @@ -25,24 +25,28 @@ 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, + raster, + 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..e972b4a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,19 +8,19 @@ export(cutStack) export(deg2rad) export(denoise) export(deseason) +export(eot) export(geoWeight) export(getWeights) export(lagalize) export(longtermMeans) +export(nmodes) export(readEot) export(subset) exportClasses(EotMode) exportClasses(EotStack) exportMethods("names<-") -exportMethods(eot) exportMethods(nXplain) exportMethods(names) -exportMethods(nmodes) exportMethods(plot) exportMethods(predict) exportMethods(writeEot) @@ -30,14 +30,11 @@ import(latticeExtra) import(mapdata) import(methods) import(parallel) -import(raster) -import(scales) +import(terra) importFrom(grDevices,colorRampPalette) importFrom(grDevices,hcl) importFrom(stats,cov.wt) importFrom(stats,na.exclude) -importFrom(stats,na.omit) -importFrom(stats,princomp) importFrom(stats,pt) importFrom(stats,var) importFrom(utils,read.csv) diff --git a/NEWS.md b/NEWS.md index 90bc6d0..c981778 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,18 @@ +# remote 1.2.3.9010 (2026-07-28) + +#### ✨ 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..cb2a48d 100644 --- a/R/EotCycle.R +++ b/R/EotCycle.R @@ -1,60 +1,56 @@ #' 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`. #' -#' @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] +#' @export +EotCycle = function( + x + , y + , n = 1 + , standardised + , orig.var + , write.out + , path.out + , prefix + , type = c("rsq", "ioa") + , verbose + , ... +) { + + type <- match.arg(type) + + x = asSpatRaster(x) + y = asSpatRaster(y) + + nl = terra::nlyr(x) + + ### Identification of the most explanatory pred pixel + + # Extract pixel entries from raster objects + x.vals <- terra::values(x) + y.vals <- terra::values(y) # 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) + a <- predRsquaredSum( + pred_vals = x.vals + , resp_vals = y.vals + , standardised = standardised + ) } else { - a <- iodaSumC(pred_vals = x.vals, resp_vals = y.vals) + a <- iodaSumC( + pred_vals = x.vals + , resp_vals = y.vals + ) } # Identify pred pixel with highest sum of r.squared @@ -64,27 +60,81 @@ EotCycle <- function(x, 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") + "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 = " ") + 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) { + 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 'y' raster template + brck.y.temp <- terra::setValues(y, values = NA_real_) + + y.lst = Map( + function(idx, len) { + v = vapply(y.lm.param.p, `[[`, numeric(len), idx) + + args = if (is.matrix(v)) { + list(x = brck.y.temp, values = t(v)) + } else { + list(x = brck.y.temp[[1L]], values = v) + } + + do.call( + terra::setValues + , args = args + ) + } + , 1:5 # parameter indexes + , c(1, 1, 1, nl, 1) # expected vector lengths + ) + + rst.y.r = y.lst[[1L]] # r + rst.y.rsq = rst.y.r ^ 2 # r-squared + rst.y.intercept = y.lst[[2L]] # intercept + rst.y.slp = y.lst[[3L]] # slope + brck.y.resids = y.lst[[4L]] # residuals + rst.y.p = y.lst[[5L]] # p + + # 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 + x.lm.param.t <- respLmParam(x.vals, x.vals, maxxy - 1) # C++ starts at 0! + # Calculate p value from T statistics + x.lm.param.p <- lapply(x.lm.param.t, function(i) { tmp <- i tmp[[5]] <- 2 * pt(-abs(tmp[[5]]), df = tmp[[6]]) @@ -94,107 +144,51 @@ EotCycle <- function(x, ## 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)) + # multi-layer 'x' raster template + brck.x.temp <- terra::setValues(x, values = NA_real_) - rst.y.r <- rst.y.rsq <- rst.y.intercept <- - rst.y.slp <- rst.y.p <- rst.y.template + x.lst = Map( + function(idx, len) { + v = vapply(x.lm.param.p, `[[`, numeric(len), idx) - # 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 + args = if (is.matrix(v)) { + list(x = brck.x.temp, values = t(v)) + } else { + list(x = brck.x.temp[[1L]], values = v) + } - ## Fit lm - - # lm(x.vals[i, ] ~ x.vals[maxxy, ]) with T statistics - x.lm.param.t <- respLmParam(x.vals, x.vals, maxxy - 1) # C++ starts at 0! - # Calculate p value from T statistics - x.lm.param.p <- lapply(x.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.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)) - - 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)) - - # R - rst.x.r[] <- sapply(x.lm.param.p, "[[", 1) - # R-squared - rst.x.rsq[] <- sapply(x.lm.param.p, "[[", 1) ^ 2 - # R-squared sums - rst.x.rsq.sums[] <- a - # Intercept - rst.x.intercept[] <- sapply(x.lm.param.p, "[[", 2) - # Slope - rst.x.slp[] <- sapply(x.lm.param.p, "[[", 3) - # P value - 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) -# } + do.call( + terra::setValues + , args = args + ) + } + , 1:5 # parameter indexes + , c(1, 1, 1, nl, 1) # expected vector lengths + ) + + rst.x.r = x.lst[[1L]] # r + rst.x.rsq = rst.x.r ^ 2 # r-squared + rst.x.intercept = x.lst[[2L]] # intercept + rst.x.slp = x.lst[[3L]] # slope + brck.x.resids = x.lst[[4L]] # residuals + rst.x.p = x.lst[[5L]] # p + + # R-squared sums + rst.x.rsq.sums = terra::setValues( + brck.x.temp[[1L]] + , values = a + ) + + # #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) @@ -204,82 +198,92 @@ EotCycle <- function(x, cat("Cum. expl. variance (%):", cum.expl.var * 100, "\n", sep = " ") } - xy <- raster::xyFromCell(x, maxxy) + 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) + 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 + ### 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) - - # 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 = ",") - } else { - 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")) - gc() - + 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 = "," + ) + } else { + 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")) + gc() + + } + return(out) } diff --git a/R/anomalize.R b/R/anomalize.R index 2268956..47565fa 100644 --- a/R/anomalize.R +++ b/R/anomalize.R @@ -1,41 +1,76 @@ -#' Create an anomaly RasterStack +#' 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 ... If `reference = NULL`, additional arguments passed to +#' [terra::app()] when calculating the overall mean (except for 'fun'). #' -#' @return an anomaly RasterStack +#' @return An anomaly `SpatRaster` series. #' #' @seealso -#' [deseason()], [denoise()], [raster::calc()] -#' -#' @export anomalize +#' [deseason()], [denoise()] #' #' @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, - ...) { +#' +#' @export +anomalize = function( + x + , reference = NULL + , ... +) { + + x = asSpatRaster(x) if (is.null(reference)) { - mn <- raster::calc(x, fun = mean, ...) + reference = terra::app( + x + , fun = mean + , ... + ) } else { - mn <- reference - } - return(x - mn) + ## early exit: 'reference' is neither `NULL` nor 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 + ) + } + + reference = asSpatRaster(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]] + } -} \ No newline at end of file + return(x - reference) +} diff --git a/R/calcVar.R b/R/calcVar.R index 02b3f05..32d4e5a 100644 --- a/R/calcVar.R +++ b/R/calcVar.R @@ -1,33 +1,34 @@ -#' Calculate space-time variance of a RasterStack or RasterBrick +#' 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 ... Currently not used. #' -#' @return the mean (optionally standardised) space-time variance. +#' @return The mean (optionally standardised) space-time variance as `numeric`. #' -#' @export calcVar +#' @examples +#' sst = terra::unwrap(pacificSST) #' -#' @examples -#' data("pacificSST") +#' calcVar(sst) # default non-standardised +#' calcVar(sst, standardised = TRUE) #' -#' calcVar(pacificSST) -calcVar <- function(x, standardised = FALSE, ...) { +#' @export +calcVar = function(x, standardised = FALSE, ...) { + + x = asSpatRaster(x) 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) + # 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(raster::getValues(x)), na.rm = TRUE) + 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..aff25ee 100644 --- a/R/cutStack.R +++ b/R/cutStack.R @@ -1,42 +1,47 @@ -#' Shorten a RasterStack +#' 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. #' -#' @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) { +#' @export +cutStack = function( + x + , tail = TRUE + , n = NULL +) { - # Return unmodified RasterStack if n == NULL + x = asSpatRaster(x) + + ## 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 { - # Supplied RasterStack is predictor: - if (tail) { - return(x[[1:(raster::nlayers(x)-n)]]) - # Supplied RasterStack 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..dd6b31d 100644 --- a/R/denoise.R +++ b/R/denoise.R @@ -1,58 +1,84 @@ #' 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 RasterStack +#' @return A denoised `SpatRaster` series. +#' +#' @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 -#' #' @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, - ...) { +#' +#' @export +denoise = function( + x + , k = NULL + , expl.var = NULL + , weighted = TRUE + , use.cpp = TRUE + , verbose = TRUE + , ... +) { - x.vals <- raster::getValues(x) + x = asSpatRaster(x) + + x.vals <- terra::values(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, ...) + pca <- princomp( + ~ x.vals + , covmat = covWeight( + x.vals + , getWeights(x) + ) + , scores = TRUE + , na.action = stats::na.exclude + , ... + ) } else { - pca <- princomp(~ x.vals, scores = TRUE, na.action = na.exclude, ...) + pca <- 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 { @@ -60,38 +86,41 @@ denoise <- function(x, } 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 = "") + fmt = paste( + "\nUsing the first %s components (of %s) to reconstruct series..." + , "these account for %s of variance in orig. series\n\n" + , sep = "\n " + ) + + txt = sprintf( + fmt + , k + , terra::nlyr(x) + , expl.var + ) + + cat(txt) } # Reconstruction - recons <- lapply(seq(nlayers(x)), function(i) { + 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] + 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) { - mat <- raster::as.matrix(x) - jnk <- insertReconsC(recons, mat) - rst <- raster::setValues(x, jnk) + jnk <- insertReconsC(recons, x.vals) + rst <- terra::setValues(x, jnk) } else { - rst <- raster::brick(lapply(seq(recons), function(i) { - tmp_rst <- x[[i]] - tmp_rst[] <- recons[[i]] - return(tmp_rst) - })) + rst = terra::setValues( + x + , values = do.call(cbind, recons) + ) } # Return denoised data set return(rst) -} +} \ No newline at end of file diff --git a/R/deseason.R b/R/deseason.R index c48917e..d294e5b 100644 --- a/R/deseason.R +++ b/R/deseason.R @@ -1,27 +1,28 @@ -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 ... Additional arguments passed to [terra::writeRaster()]. Only +#' considered if 'filename' is specified. #' -#' @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,74 +31,93 @@ 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) ################################################################################ -### function using 'RasterStack' or 'RasterBrick' ############################## -#' @aliases deseason,RasterStackBrick-method +### function using 'SpatRaster' ################################################ +#' @aliases deseason,SpatRaster-method #' @rdname deseason -setMethod("deseason", - signature(x = "RasterStackBrick"), - function(x, - cycle.window = 12L, - use.cpp = FALSE, - filename = "", - ...) { - - if (use.cpp) { - ## raster to matrix - mat <- as.matrix(x) - - ## deseasoning - mat_mv <- monthlyMeansC(mat, cycle.window) - x_mv <- x[[1:cycle.window]] - x_mv <- 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)]], - fun = mean, na.rm = TRUE) - }), raster::nlayers(x) / cycle.window)) - } - - # Subtract monthly averages from actually measured values - x_dsn <- x - x_mv - - # Write to file (optional) - if (filename != "") - x_dsn <- raster::writeRaster(x_dsn, filename = filename, ...) - - # Return output - return(x_dsn) - }) +methods::setMethod( + "deseason" + , signature(x = "SpatRaster") + , function( + x + , cycle.window = 12L + , use.cpp = FALSE + , filename = "" + , ... + ) { + + x = asSpatRaster(x) + + if (use.cpp) { + ## raster to matrix + mat <- terra::as.matrix(x) + + ## deseasoning + mat_mv <- monthlyMeansC(mat, cycle.window) + x_mv <- x[[1:cycle.window]] + x_mv <- terra::setValues(x_mv, values = mat_mv) + } else { + # Calculate layer averages based on supplied seasonal 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) + })) + } + + # Subtract monthly averages from actually measured values + x_dsn <- x - x_mv + + # Write to file (optional) + if (filename != "") + x_dsn <- terra::writeRaster(x_dsn, filename = filename, ...) + + # Return output + return(x_dsn) + } +) ################################################################################ ### function using 'numeric' ################################################### #' @aliases deseason,numeric-method #' @rdname deseason -setMethod("deseason", - signature(x = "numeric"), - function(x, - cycle.window = 12L) { - - ## calculate long-term mean values - x_mv <- sapply(1:cycle.window, function(i) { - val <- x[seq(i, length(x), cycle.window)] - mean(val, na.rm = TRUE) - }) - - ## create anomalies - x_dsn <- x - x_mv - - # Return output - return(x_dsn) - }) +methods::setMethod( + "deseason" + , signature(x = "numeric") + , function( + x + , cycle.window = 12L + ) { + + ## construct index for splitting time series into seasonal windows + idx = rep( + 1:cycle.window + , length.out = length(x) + ) + + ## calculate long-term mean values + x_mv = vapply( + split( + x + , f = idx + ) + , FUN = mean + , FUN.VALUE = numeric(1L) + , USE.NAMES = FALSE + ) + + ## create anomalies + x_dsn = x - x_mv + + return(x_dsn) + } +) diff --git a/R/eot.R b/R/eot.R index 766dd48..9d40511 100644 --- a/R/eot.R +++ b/R/eot.R @@ -1,40 +1,35 @@ -if ( !isGeneric('eot') ) { - 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 +37,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 +71,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 +88,112 @@ 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) +# 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) - } -) - -# #' @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) -# } -# ) -# +eot = 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) + + x = asSpatRaster(x) + y = asSpatRaster(y) # returns `NULL` unaltered + + ## 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..03431d0 100644 --- a/R/geoWeight.R +++ b/R/geoWeight.R @@ -1,37 +1,39 @@ #' 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()]. #' -#' @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 -#' -#' @export geoWeight +#' @return A weighted `SpatRaster`. #' #' @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), - ...) { +#' +#' @export +geoWeight = 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 = asSpatRaster(x) x.vals <- x[] - rads <- deg2rad(sp::coordinates(x)[, 2]) + 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..8821351 100644 --- a/R/getWeights.R +++ b/R/getWeights.R @@ -1,31 +1,47 @@ #' 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'. #' -#' @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) #' -#' opar <- par(mfrow = c(1,2)) -#' plot(australiaGPCP[[1]], main = "data") +#' wghts = getWeights(pcp) +#' utils::head(wghts) +#' +#' wghts_rst = terra::setValues(pcp[[1]], wghts) +#' +#' 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), - ...) { +#' @export +getWeights = function( + x + , f = cos + , ... +) { - f(deg2rad(sp::coordinates(x)[, 2][!is.na(x[[1]][])]), ...) + # 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 + x = asSpatRaster(x) -} \ No newline at end of file + f( + deg2rad( + terra::crds(x, na.rm = TRUE)[, 2L] + ) + , ... + ) + +} diff --git a/R/lagalize.R b/R/lagalize.R index dd475bb..f41a7d0 100644 --- a/R/lagalize.R +++ b/R/lagalize.R @@ -1,51 +1,57 @@ -#' Create lagged RasterStacks +#' 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 ... Currently not used. #' #' @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, - ...) { +#' @export +lagalize = function( + x + , y + , lag = NULL + , freq = 12L + , ... +) { - # Return list of unmodified RasterStacks if lag == NULL + x = asSpatRaster(x) + y = asSpatRaster(y) + + # Return list of unmodified raster series if lag == NULL if (is.null(lag)) { return(list(x, y)) - } else { - 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)]] - - # Lagalize response stack - y.lag <- cutStack(x = y, tail = FALSE, n = lag) - y.lag.adj <- y.lag[[1:(raster::nlayers(y.lag) - rest)]] - - # Return list of lagalized stacks - return(list(x.lag, y.lag)) } -} \ No newline at end of file + rest <- freq - lag + + # Lagalize predictor stack + x.lag <- cutStack(x = x, tail = TRUE, n = lag) + 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:(terra::nlyr(y.lag) - rest)]] + + # Return list of lagalized stacks + return(list(x.lag, y.lag)) + +} diff --git a/R/longtermMeans.R b/R/longtermMeans.R index 62b5d26..0eebb14 100644 --- a/R/longtermMeans.R +++ b/R/longtermMeans.R @@ -1,17 +1,17 @@ -#' Calculate long-term means from a 'RasterStack' +#' 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()]. #' #' @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 +20,27 @@ #' [deseason()]. #' #' @examples -#' data("australiaGPCP") +#' pcp = terra::unwrap(australiaGPCP) #' -#' longtermMeans(australiaGPCP) +#' longtermMeans(pcp) #' -#' @export longtermMeans -#' @name longtermMeans -longtermMeans <- function(x, cycle.window = 12L) { - - ## raster to matrix - mat <- raster::as.matrix(x) +#' @export +longtermMeans = function(x, cycle.window = 12L) { - ## long-term means - mat_ltm <- monthlyMeansC(mat, 12) + x = asSpatRaster(x) ## insert values - rst_ltm <- x[[1:(raster::nlayers(x) / cycle.window)]] - rst_ltm <- raster::setValues(rst_ltm, NA) + idx = rep( + 1:cycle.window + , times = terra::nlyr(x) / cycle.window + ) + + rst_ltm = terra::tapp( + x + , index = idx + , fun = mean + , na.rm = TRUE + ) - rst_ltm <- raster::setValues(rst_ltm, mat_ltm) return(rst_ltm) } diff --git a/R/nXplain.R b/R/nXplain.R index a756dea..5c33d69 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) #' @@ -35,21 +37,44 @@ if (!isGeneric('nXplain')) { #' @rdname nXplain #' @aliases nXplain,EotStack-method -setMethod('nXplain', signature(x = 'EotStack'), - function(x, var = 0.9) { - expl.var <- sapply(seq(nmodes(x)), function(i) { - 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) - } - - return(n) - } +methods::setMethod( + 'nXplain' + , signature(x = 'EotStack') + , function( + x + , var = 0.9 + ) { + + expl.var = vapply( + x@modes + , FUN = slot + , FUN.VALUE = numeric(1L) + , name = "cum_exp_var" + ) + + idx = var - expl.var <= 0 + + if (!any(idx)) { + fmt = paste( + "explained variance of EotStack is lower than: %s" + , "maximum explained variance of this EotStack is: %s" + , sep = "\n" + ) + + txt = sprintf( + fmt + , var + , x[[nmodes(x)]]@cum_exp_var + ) + + stop( + txt + , 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..2cea260 100644 --- a/R/plot.R +++ b/R/plot.R @@ -32,38 +32,38 @@ 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, +#' ## calculate 2 leading modes +#' nh_modes <- eot(x = gph, y = NULL, n = 2, #' standardised = FALSE, #' verbose = TRUE) #' -#' ## default settings -#' plot(nh_modes, y = 1) # is equivalent to -#' -#' \dontrun{ -#' plot(nh_modes[[1]]) -#' -#' plot(nh_modes, y = 2) # shows variance explained by mode 2 only -#' plot(nh_modes[[2]]) # shows cumulative variance explained by modes 1 & 2 -#' -#' ## showing the loction of the mode -#' plot(nh_modes, y = 1, show.bp = TRUE) -#' -#' ## changing parameters -#' plot(nh_modes, y = 1, show.bp = TRUE, -#' pred.prm = "r", resp.prm = "p") -#' -#' ## change plot arrangement -#' plot(nh_modes, y = 1, show.bp = TRUE, arrange = "long") -#' -#' ## plot locations of all base points -#' plot(nh_modes, locations = TRUE) -#' } +# ## default settings +# plot(nh_modes, y = 1) # is equivalent to +# +# \dontrun{ +# plot(nh_modes[[1]]) +# +# plot(nh_modes, y = 2) # shows variance explained by mode 2 only +# plot(nh_modes[[2]]) # shows cumulative variance explained by modes 1 & 2 +# +# ## showing the loction of the mode +# plot(nh_modes, y = 1, show.bp = TRUE) +# +# ## changing parameters +# plot(nh_modes, y = 1, show.bp = TRUE, +# pred.prm = "r", resp.prm = "p") +# +# ## change plot arrangement +# plot(nh_modes, y = 1, show.bp = TRUE, arrange = "long") +# +# ## plot locations of all base points +# plot(nh_modes, locations = TRUE) +# } #' #' @export #' @name plot @@ -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..cbd44ef 100644 --- a/R/predict.R +++ b/R/predict.R @@ -1,141 +1,151 @@ -# 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 = vapply( + object@modes[1:n] + , FUN = slot + , FUN.VALUE = integer(1L) + , name = "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/print.R b/R/print.R index 6cecd4a..e2c37bb 100644 --- a/R/print.R +++ b/R/print.R @@ -1,18 +1,14 @@ -setMethod ('print', 'EotMode', +methods::setMethod ('print', 'EotMode', function(x, ...) { - if (inherits(x, 'EotMode')) { - show(x) - } + show(x) } ) -setMethod ('print', 'EotStack', +methods::setMethod ('print', 'EotStack', function(x, ...) { - if (inherits(x, 'EotStack')) { - show(x) - } + show(x) } ) \ No newline at end of file diff --git a/R/readEot.R b/R/readEot.R index 8493e11..d8cec43 100644 --- a/R/readEot.R +++ b/R/readEot.R @@ -1,61 +1,95 @@ -#' 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 + ) + ) + + mds = unique(mds) + ## 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()` + fl_mds = list.files( + x + , sprintf( + "^%s_eot_locations\\.csv$" + , prefix + ) + , full.names = TRUE + ) + + dat_mds = utils::read.csv(fl_mds) ## loop over modes, creating 'EotMode' objects for each mode available lst_eot <- lapply(mds, function(n) { # track and reorder files related to current mode - fls <- fls_mds[grep(paste0("mode_", n), basename(fls_mds))] - ids <- sapply(eotLayerNames(), function(j) grep(j, fls)) - fls <- fls[ids] + pttrns = sprintf( + "mode_%s%s" + , n + , eotLayerNames() + ) + + fls = vapply( + pttrns + , FUN = grep + , FUN.VALUE = character(1L) + , x = fls_mds + , value = TRUE + ) # 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 +116,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) <- vapply(lst_eot, slot, character(1L), "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..c6cdb52 100644 --- a/R/remote-package.R +++ b/R/remote-package.R @@ -1,8 +1,9 @@ #' R EMpirical Orthogonal TEleconnections #' -#' A collection of functions to facilitate empirical orthogonal teleconnection analysis. -#' Some handy functions for preprocessing, such as deseasoning, denoising, lagging -#' are readily available for ease of usage. +#' @description +#' A collection of functions to facilitate empirical orthogonal teleconnection +#' analysis. Some handy functions for preprocessing, such as deseasoning, +#' denoising, lagging are readily available for ease of usage. #' #' @name remote-package #' @aliases remote @@ -20,11 +21,17 @@ #' Empirical methods in short-term climate prediction\cr #' H. M. van den Dool (2007)\cr #' 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 +#' +#' @seealso \pkg{remote} is built upon the \pkg{terra} package, which is the +#' direct successor of the \pkg{raster} package. The \pkg{terra} package is used +#' for raster data handling and processing, while \pkg{raster} is still +#' supported but not actively developed anymore. For more information on raster +#' data handling, please refer to the documentation of the \pkg{terra} package +#' and its functions. +#' +#' @import Rcpp gridExtra latticeExtra mapdata methods parallel terra #' @importFrom grDevices colorRampPalette hcl -#' @importFrom stats pt var cov.wt na.exclude princomp na.omit +#' @importFrom stats pt var cov.wt na.exclude #' @importFrom utils read.csv write.table #' @useDynLib remote #' @@ -33,15 +40,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 +59,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/utils.R b/R/utils.R new file mode 100644 index 0000000..3d496cc --- /dev/null +++ b/R/utils.R @@ -0,0 +1,9 @@ +asSpatRaster = function(x) { + + ## if applicable, convert `Raster*` input to `SpatRaster` + if (inherits(x, what = "Raster")) { + x = terra::rast(x) + } + + return(x) +} \ No newline at end of file diff --git a/R/writeEot.R b/R/writeEot.R index 8133f25..c138b8d 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( + function(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/README.md b/README.md index bbefd2f..2ca85e5 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ ### R EMpirical Orthogonal TEleconnections -[![R-CMD-check](https://github.com/environmentalinformatics-marburg/remote/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/environmentalinformatics-marburg/remote/actions/workflows/R-CMD-check.yaml) +[![R-CMD-check](https://github.com/r-spatial/remote/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/r-spatial/remote/actions/workflows/R-CMD-check.yaml) [![cranchecks](https://badges.cranchecks.info/worst/remote.svg)](https://cran.r-project.org/web/checks/check_results_remote.html) ![monthly](https://cranlogs.r-pkg.org/badges/remote) ![total](https://cranlogs.r-pkg.org/badges/grand-total/remote) 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..a4fe1f5 100644 --- a/dev_history.R +++ b/dev_history.R @@ -19,3 +19,215 @@ 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") + + +# 2026-07-28 ==== + +## RESOLVE CHECK ISSUES ==== + +## eliminate "Compilation used the following non-portable flag(s): +## ‘-mno-omit-leaf-frame-pointer’" note +usethis::edit_r_makevars() + +## add these lines: +## CFLAGS = -g -O2 -Wall -pedantic +## CXXFLAGS = -g -O2 -Wall -pedantic \ No newline at end of file diff --git a/inst/tinytest/test_anomalize.R b/inst/tinytest/test_anomalize.R new file mode 100644 index 0000000..61995c0 --- /dev/null +++ b/inst/tinytest/test_anomalize.R @@ -0,0 +1,52 @@ +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" +) + +## custom 'reference' input +ref_md = terra::app( + pcp + , fun = "median" +) + +expect_false( + all( + (anm - ref_md)[[1L]][] == 0L + ) + , info = "returns different result if 'reference' is not overall mean" +) + +## 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 = terra::rast(replicate(2L, ref)) + ) + , 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..ae76ae1 --- /dev/null +++ b/inst/tinytest/test_calcVar.R @@ -0,0 +1,29 @@ +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 +) 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..db04f5d --- /dev/null +++ b/inst/tinytest/test_cutStack.R @@ -0,0 +1,46 @@ +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" +) 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..119ff86 --- /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(pcp, 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 `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..c2e3528 --- /dev/null +++ b/inst/tinytest/test_deseason.R @@ -0,0 +1,45 @@ +### `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" +) + +## 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" +) + + +### `numeric` ---- + +pcp_vls = terra::global( + pcp + , fun = "mean" +)[, 1L] + +deseason(pcp_vls) diff --git a/inst/tinytest/test_eot.R b/inst/tinytest/test_eot.R new file mode 100644 index 0000000..6494de9 --- /dev/null +++ b/inst/tinytest/test_eot.R @@ -0,0 +1,365 @@ +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`" +) + +## `type = "ioa"` +nh_modes_ioa <- eot( + x = gph + , y = NULL + , n = n + , standardised = FALSE + , type = "ioa" + , verbose = FALSE +) + +expect_true( + nh_modes_ioa[[n]]@cum_exp_var != nh_modes[[n]]@cum_exp_var + , info = "returns different results for `type = 'ioa'`" +) + +## errors +expect_error( + eot( + x = gph + , type = "rmse" + ) + , pattern = "'arg' should be one of" +) + + +### `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..e7e40e9 --- /dev/null +++ b/inst/tinytest/test_geoWeight.R @@ -0,0 +1,28 @@ +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" +) diff --git a/inst/tinytest/test_getWeights.R b/inst/tinytest/test_getWeights.R new file mode 100644 index 0000000..e5c91a0 --- /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(pcp1) + +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..a58ac12 --- /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_l0 = lagalize(sst, pcp, freq = 12L) + +expect_equal( + lagged_l0[[1L]] + , target = sst + , info = "returns the original 'x' series when 'lag' is `NULL`" +) + +expect_equal( + lagged_l0[[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..94cd159 --- /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(pcp, 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/inst/tinytest/test_utils.R b/inst/tinytest/test_utils.R new file mode 100644 index 0000000..0e9164d --- /dev/null +++ b/inst/tinytest/test_utils.R @@ -0,0 +1,41 @@ +### `asSpatRaster()` ----- + +## `Raster*` input +r = system.file( + "external/rlogo.grd" + , package = "raster" +) |> + raster::brick() + +s = remote:::asSpatRaster(r) + +expect_inherits( + s + , class = "SpatRaster" + , info = "converts `Raster*` to `SpatRaster`" +) + +expect_equal( + dim(s) + , dim(r) + , info = "preserves dimensions when converting `Raster*` to `SpatRaster`" +) + +## `SpatRaster` input +s1 = system.file( + "ex/elev.tif" + , package = "terra" +) |> + terra::rast() + +expect_identical( + remote:::asSpatRaster(s1) + , target = s1 + , info = "returns `SpatRaster` input unchanged" +) + +## `NULL` input +expect_null( + remote:::asSpatRaster(NULL) + , info = "returns `NULL` input unchanged" +) diff --git a/man/EotCycle.Rd b/man/EotCycle.Rd index ac2e6c5..4484ceb 100644 --- a/man/EotCycle.Rd +++ b/man/EotCycle.Rd @@ -13,43 +13,18 @@ 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{standardised}{logical. If \code{FALSE} the calculated r-squared values -will be multiplied by the variance} +\item{x, y, n, standardised, write.out, path.out, prefix, type, verbose}{See \code{\link[=eot]{eot()}}.} \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()}}.} +\item{...}{Further arguments passed to \code{\link[=writeEot]{writeEot()}} if \code{write.out = TRUE}.} } \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..6823f27 100644 --- a/man/anomalize.Rd +++ b/man/anomalize.Rd @@ -2,36 +2,38 @@ % Please edit documentation in R/anomalize.R \name{anomalize} \alias{anomalize} -\title{Create an anomaly RasterStack} +\title{Create an anomaly raster series} \usage{ anomalize(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{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).} -\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{...}{If \code{reference = NULL}, additional arguments passed to +\code{\link[terra:app]{terra::app()}} when calculating the overall mean (except for 'fun').} } \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) +pcp = terra::unwrap(australiaGPCP) +pcp_anom = anomalize(pcp) -aus_anom <- anomalize(australiaGPCP) - -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..072b23e 100644 --- a/man/calcVar.Rd +++ b/man/calcVar.Rd @@ -2,26 +2,28 @@ % Please edit documentation in R/calcVar.R \name{calcVar} \alias{calcVar} -\title{Calculate space-time variance of a RasterStack or RasterBrick} +\title{Calculate space-time variance of a raster series} \usage{ calcVar(x, standardised = FALSE, ...) } \arguments{ -\item{x}{a RasterStack or RasterBrick} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} -\item{standardised}{logical.} +\item{standardised}{\code{logical}, defaults to \code{FALSE}.} -\item{...}{currently not used} +\item{...}{Currently not used.} } \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(sst) # default non-standardised +calcVar(sst, standardised = TRUE) -calcVar(pacificSST) } 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..416410f 100644 --- a/man/cutStack.Rd +++ b/man/cutStack.Rd @@ -2,33 +2,34 @@ % Please edit documentation in R/cutStack.R \name{cutStack} \alias{cutStack} -\title{Shorten a RasterStack} +\title{Shorten a raster series} \usage{ cutStack(x, tail = TRUE, n = NULL) } \arguments{ -\item{x}{a RasterStack} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} -\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..caf8d6e 100644 --- a/man/denoise.Rd +++ b/man/denoise.Rd @@ -15,41 +15,46 @@ denoise( ) } \arguments{ -\item{x}{RasterStack to be filtered} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series to be filtered.} -\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{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{...}{Additional arguments passed to \code{\link[stats:princomp]{stats::princomp()}}.} } \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{ \code{\link[=anomalize]{anomalize()}}, \code{\link[=deseason]{deseason()}} diff --git a/man/deseason.Rd b/man/deseason.Rd index 17174b4..a0b821a 100644 --- a/man/deseason.Rd +++ b/man/deseason.Rd @@ -2,47 +2,45 @@ % Please edit documentation in R/deseason.R \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}{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{cycle.window}{\code{integer}, defaults to \code{12}. The window for the creation +of the anomalies.} -\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{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{filename}{\code{character}. Output filename (optional).} -\item{...}{Additional arguments passed on to \code{\link[raster:writeRaster]{raster::writeRaster()}}, only -considered if \code{filename} is specified.} +\item{...}{Additional arguments passed to \code{\link[terra:writeRaster]{terra::writeRaster()}}. Only +considered if 'filename' is specified.} } \value{ -If \code{x} is a \verb{Raster*} object, a deseasoned -\code{RasterStack}; else a deseasoned \code{numeric} vector. +If '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") +pcp = terra::unwrap(australiaGPCP) +pcp_dsn = deseason(pcp, cycle.window = 12) -aus_dsn <- deseason(australiaGPCP, 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..1bf45a7 100644 --- a/man/eot.Rd +++ b/man/eot.Rd @@ -2,10 +2,9 @@ % Please edit documentation in R/eot.R \name{eot} \alias{eot} -\alias{eot,RasterStackBrick-method} -\title{EOT analysis of a predictor and (optionally) a response RasterStack} +\title{EOT analysis of a predictor and (optionally) a response raster series} \usage{ -\S4method{eot}{RasterStackBrick}( +eot( x, y = NULL, n = 1, @@ -20,40 +19,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{y}{A \code{SpatRaster} (or \verb{Raster*}) object used as response. If 'y' is +\code{NULL} (default), 'x' is used as response.} -\item{n}{the number of EOT modes to calculate} +\item{n}{The number of EOT modes to calculate as \code{integer}, defaults to \code{1}.} -\item{standardised}{logical. If \code{FALSE} the calculated r-squared values -will be multiplied by the variance} +\item{standardised}{\code{logical}, default \code{TRUE}. If \code{FALSE} the calculated +r-squared values will be multiplied by the variance.} -\item{write.out}{logical. If \code{TRUE} results will be written to disk -using \code{path.out}} +\item{write.out}{\code{logical}, default \code{FALSE}. If \code{TRUE} results will be +written to disk using 'path.out'.} -\item{path.out}{the file path for writing results if \code{write.out} is \code{TRUE}. -Defaults to current working directory} +\item{path.out}{The file path for writing results if '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{prefix}{optional prefix to be used for naming of results if 'write.out' +is \code{TRUE}.} -\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{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{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{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{verbose}{logical. If \code{TRUE} some details about the -calculation process will be output to the console} +\item{verbose}{\code{logical}, default \code{FALSE}. If \code{TRUE} some details about the +calculation process will be output to the console.} -\item{...}{not used at the moment} +\item{...}{Currently not used.} } \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 +60,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 +78,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 +105,15 @@ 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..bda6ce6 100644 --- a/man/geoWeight.Rd +++ b/man/geoWeight.Rd @@ -4,31 +4,26 @@ \alias{geoWeight} \title{Geographic weighting} \usage{ -geoWeight(x, f = function(x) cos(x), ...) +geoWeight(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()}}.} } \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..cefbd8e 100644 --- a/man/getWeights.Rd +++ b/man/getWeights.Rd @@ -4,31 +4,33 @@ \alias{getWeights} \title{Calculate weights from latitude} \usage{ -getWeights(x, f = function(x) cos(x), ...) +getWeights(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{f}{A \code{function} applied to the latitude (in radians) to compute +weights. Defaults to \code{cos}.} -\item{...}{additional arguments to be passed to f} +\item{...}{Additional arguments passed to 'f'.} } \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) -opar <- par(mfrow = c(1,2)) -plot(australiaGPCP[[1]], main = "data") +wghts = getWeights(pcp) +utils::head(wghts) + +wghts_rst = terra::setValues(pcp[[1]], wghts) + +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..b4e79f0 100644 --- a/man/lagalize.Rd +++ b/man/lagalize.Rd @@ -2,36 +2,37 @@ % Please edit documentation in R/lagalize.R \name{lagalize} \alias{lagalize} -\title{Create lagged RasterStacks} +\title{Create lagged raster series} \usage{ -lagalize(x, y, lag = NULL, freq = 12, ...) +lagalize(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{lag}{The desired lag in the native frequency of the series passed to +\code{\link[=cutStack]{cutStack()}}.} -\item{freq}{the frequency of the RasterStacks} +\item{freq}{The frequency of the raster series as \code{integer}.} -\item{...}{currently not used} +\item{...}{Currently not used.} } \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..bfa4c6c 100644 --- a/man/longtermMeans.Rd +++ b/man/longtermMeans.Rd @@ -2,29 +2,29 @@ % Please edit documentation in R/longtermMeans.R \name{longtermMeans} \alias{longtermMeans} -\title{Calculate long-term means from a 'RasterStack'} +\title{Calculate long-term means from a raster series} \usage{ longtermMeans(x, cycle.window = 12L) } \arguments{ -\item{x}{A 'RasterStack' (or 'RasterBrick') object.} +\item{x}{A \code{SpatRaster} (or \verb{Raster*}) series.} -\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..3de6778 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,34 +82,12 @@ 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, +## calculate 2 leading modes +nh_modes <- eot(x = gph, y = NULL, n = 2, standardised = FALSE, verbose = TRUE) -## default settings -plot(nh_modes, y = 1) # is equivalent to - -\dontrun{ -plot(nh_modes[[1]]) - -plot(nh_modes, y = 2) # shows variance explained by mode 2 only -plot(nh_modes[[2]]) # shows cumulative variance explained by modes 1 & 2 - -## showing the loction of the mode -plot(nh_modes, y = 1, show.bp = TRUE) - -## changing parameters -plot(nh_modes, y = 1, show.bp = TRUE, - pred.prm = "r", resp.prm = "p") - -## change plot arrangement -plot(nh_modes, y = 1, show.bp = TRUE, arrange = "long") - -## plot locations of all base points -plot(nh_modes, locations = 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/remote-package.Rd b/man/remote-package.Rd index b6b9232..ea80c15 100644 --- a/man/remote-package.Rd +++ b/man/remote-package.Rd @@ -6,12 +6,12 @@ \alias{remote} \title{R EMpirical Orthogonal TEleconnections} \description{ -R EMpirical Orthogonal TEleconnections +A collection of functions to facilitate empirical orthogonal teleconnection +analysis. Some handy functions for preprocessing, such as deseasoning, +denoising, lagging are readily available for ease of usage. } \details{ -A collection of functions to facilitate empirical orthogonal teleconnection analysis. -Some handy functions for preprocessing, such as deseasoning, denoising, lagging -are readily available for ease of usage. +R EMpirical Orthogonal TEleconnections } \references{ Empirical Orthogonal Teleconnections\cr @@ -23,8 +23,12 @@ H. M. van den Dool (2007)\cr Oxford University Press, Oxford, New York (2007)\cr } \seealso{ -\pkg{remote} is built upon Raster* classes from the -\link[raster:raster-package]{raster::raster-package}. Please see their documentation for data preparation etc. +\pkg{remote} is built upon the \pkg{terra} package, which is the +direct successor of the \pkg{raster} package. The \pkg{terra} package is used +for raster data handling and processing, while \pkg{raster} is still +supported but not actively developed anymore. For more information on raster +data handling, please refer to the documentation of the \pkg{terra} package +and its functions. } \author{ Tim Appelhans, Florian Detsch, Thomas Nauss\cr 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/remote.Rproj b/remote.Rproj index b2980c2..0fa49a8 100644 --- a/remote.Rproj +++ b/remote.Rproj @@ -14,9 +14,8 @@ RnwWeave: Sweave LaTeX: pdfLaTeX BuildType: Package -PackageUseDevtools: Yes PackageCleanBeforeInstall: No PackageInstallArgs: --no-multiarch --with-keep.source -PackageBuildArgs: --resave-data -PackageCheckArgs: --as-cran --no-manual +PackageBuildArgs: --resave-data=best +PackageCheckArgs: --as-cran --no-manual --run-donttest PackageRoxygenize: rd,collate,namespace diff --git a/src/Makevars b/src/Makevars deleted file mode 100644 index c880859..0000000 --- a/src/Makevars +++ /dev/null @@ -1,27 +0,0 @@ -## Use the R_HOME indirection to support installations of multiple R version -PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` - -## As an alternative, one can also add this code in a file 'configure' -## -## PKG_LIBS=`${R_HOME}/bin/Rscript -e "Rcpp:::LdFlags()"` -## -## sed -e "s|@PKG_LIBS@|${PKG_LIBS}|" \ -## src/Makevars.in > src/Makevars -## -## which together with the following file 'src/Makevars.in' -## -## PKG_LIBS = @PKG_LIBS@ -## -## can be used to create src/Makevars dynamically. This scheme is more -## powerful and can be expanded to also check for and link with other -## libraries. It should be complemented by a file 'cleanup' -## -## rm src/Makevars -## -## which removes the autogenerated file src/Makevars. -## -## Of course, autoconf can also be used to write configure files. This is -## done by a number of packages, but recommended only for more advanced users -## comfortable with autoconf and its related tools. - - diff --git a/src/Makevars.win b/src/Makevars.win deleted file mode 100644 index 18d49b4..0000000 --- a/src/Makevars.win +++ /dev/null @@ -1,3 +0,0 @@ - -## Use the R_HOME indirection to support installations of multiple R version -PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()") 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") +} +