diff --git a/.Rbuildignore b/.Rbuildignore index e66ed82d..b9cd0abd 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -12,13 +12,12 @@ CODE_OF_CONDUCT.md ^\.vscode ^README\.Rmd$ ^.lintr$ -^fasterRaster.code-workspace ^_pkgdown\.yml$ ^docs$ ^pkgdown$ ^\.github$ ^\vignettes\^junk junk -fasterRaster_workspace.code-workspace +workspace_fasterRaster.code-workspace ^doc$ ^Meta$ diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths index e69de29b..8a7abe28 100644 --- a/.Rproj.user/shared/notebooks/paths +++ b/.Rproj.user/shared/notebooks/paths @@ -0,0 +1,2 @@ +C:/Kaji/R/fasterRaster/DESCRIPTION="BE5F8589" +C:/Kaji/R/fasterRaster/NEWS.md="B2E56A81" diff --git a/DESCRIPTION b/DESCRIPTION index f57c234c..e3a7f2f1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: fasterRaster Type: Package Title: Faster Raster and Spatial Vector Processing Using 'GRASS' -Version: 8.4.1.0 -Date: 2025-06-19 +Version: 8.4.1.1 +Date: 2025-11-18 Authors@R: c( person( @@ -47,4 +47,4 @@ VignetteBuilder: knitr Encoding: UTF-8 LazyData: true LazyLoad: yes -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 diff --git a/NEWS.md b/NEWS.md index 4c0fd09e..75ba301d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +# fasterRaster 8.4.1.1 (2025-11-18) +### New functions and functionality +o `+` can now combine any two `GVector`s with the same geometry type (points, lines, polygons). +o `installAddon()` checks to see if addon is already installed, and if so, uninstalls it before re-installing. +o Minor documentation updates. + +### Bug fixes +o `fast()` correctly adds names to `GRaster`s loaded from disk when a `_names.csv` file exists. +o `geomorphons()` now works when `mode = '2'`. +o `rbind()` combines `GVector` data tables when each vecto has the same number of rows in its table as geometries. +o `scale()` and `scalepop()` work when `center` and/or `scale` are numeric. +o `sf` example data objects (`madCoast0`, `madCoast4`, and `madRivers`) have had their CRSs harmonized with the example raster objects. This had caused some of the examples not to work. + +### Other +o Added error checking in `fragmentation()` related to limits on size of window `w` for calculation of fragmentation indices. + # fasterRaster 8.4.1.0 (2025-06-17) ### Code-breaking changes diff --git a/R/04_arithmetic.r b/R/04_arithmetic.r index 5a1fc93b..01cdc02d 100644 --- a/R/04_arithmetic.r +++ b/R/04_arithmetic.r @@ -8,7 +8,7 @@ #' `*` operator: Same as [intersect()]\cr #' `/` operator: Same as [xor()]\cr #' -#' @param e1,e2 `GRaster`s, `numeric`s, `integer`s, or `logical`s. +#' @param e1,e2 `GRaster`s, `GVector`s, `numeric`s, `integer`s, or `logical`s. #' #' @return A `GRaster`. #' @@ -428,13 +428,21 @@ methods::setMethod( signature(e1 = "GVector", e2 = "GVector"), function(e1, e2) { + gtype1 <- geomtype(e1) + gtype2 <- geomtype(e2) + if (gtype1 != gtype2) stop("Vectors must be of the same type (polygons, lines, or points).") + compareGeom(e1, e2) .locationRestore(e1) oper <- as.vector(.Generic)[1L] if (oper == "+") { - out <- union(e1, e2) + if (gtype1 == "polygons" & gtype2 == "polygons") { + out <- union(e1, e2) + } else { + out <- rbind(e1, e2) + } } else if (oper == "-") { out <- erase(e1, e2) } else if (oper == "*") { diff --git a/R/addons.r b/R/addons.r index 5dbb362c..550726f3 100644 --- a/R/addons.r +++ b/R/addons.r @@ -36,15 +36,24 @@ addons <- function(x = NULL) { installAddon <- function(x, check = TRUE) { if (check) { + + # is it available? avails <- rgrass::execGRASS("g.extension", flags = c("l", .quiet()), intern = TRUE) if (!(x %in% avails)) { + warning("The addon is not available on the official GRASS addon repository.") return(invisible(FALSE)) + } } + # remove addon before re-installing + if (addons(x)) removeAddon(x, check = FALSE) + rgrass::execGRASS("g.extension", operation = "add", extension = x, flags = .quiet()) - invisible(TRUE) + success <- addons(x) + if (!success) warning("Addon could not be installed.") + invisible(success) } diff --git a/R/backdoor.r b/R/backdoor.r index e8eac95a..cf2ab1c1 100644 --- a/R/backdoor.r +++ b/R/backdoor.r @@ -2,22 +2,32 @@ #' #' This is a secret function to be used for faster development of **fasterRaster**. It calls [faster()] to set the install directory for **GRASS**, increases default memory, and number of cores. The function assumes development is on a Windows machine. #' +#' @param start Logical: If `TRUE`, start the **GRASS** session by creating the `madElev` `GRaster`. #' @param ver Character: **GRASS**: e.g., "83" or "84". #' #' @returns `TRUE` (invisibly). #' #' @keywords internal -.backdoor <- function(ver = "84") { +.backdoor <- function(start = FALSE, ver = "84") { verNice <- paste0(substr(ver, 1L, 1L), ".", substr(ver, 2L, 2L)) faster( grassDir = paste0("C:/Program Files/GRASS GIS ", verNice), - addonsDir = 'C:/Users/adame/AppData/Roaming/GRASS8/addons', + addonsDir = "C:/Users/adame/AppData/Roaming/GRASS8/addons", memory = 1024 * 8, cores = 2, useDataTable = TRUE, verbose = TRUE ) - invisible(TRUE) + + if (start) { + madElev <- fastData("madElev") + elev <- fast(madElev) + out <- elev + } else { + out <- TRUE + } + + invisible(out) } diff --git a/R/fast.r b/R/fast.r index 7a09fa81..c7a4b7f6 100644 --- a/R/fast.r +++ b/R/fast.r @@ -179,7 +179,16 @@ methods::setMethod( xRast <- terra::rast(x) nLayers <- terra::nlyr(xRast) + xNames <- names(xRast) + fileExtension <- .fileExt(x) + nc <- nchar(fileExtension) + 1 + namesFile <- substr(x, 1, nchar(x) - nc) + namesFile <- paste0(namesFile, '_names.csv') + if (file.exists(namesFile)) { + xNames <- utils::read.csv(namesFile) + xNames <- xNames$name + } location <- .locationFind(xRast, match = "crs") diff --git a/R/fastData.r b/R/fastData.r index 244ab4a9..7b7c04d1 100644 --- a/R/fastData.r +++ b/R/fastData.r @@ -5,24 +5,24 @@ #' @param x The name of the raster or spatial vector to get. All of these represent a portion of the eastern coast of Madagascar. #' #' Spatial vectors (objects of class `sf` from the **sf** package): -#' * [madCoast0]: Outline of the region (polygon) -#' * [madCoast4]: Outlines of the Fokontanies (Communes) of the region (polygons) -#' * [madDypsis]: Records of plants of the genus *Dypsis* (points) -#' * [madRivers]: Major rivers (lines) +#' * [madCoast0][madCoast0]: Outline of the region (polygon) +#' * [madCoast4][madCoast4]: Outlines of the Fokontanies (Communes) of the region (polygons) +#' * [madDypsis][madDypsis]: Records of plants of the genus *Dypsis* (points) +#' * [madRivers][madRivers]: Major rivers (lines) #' #' Rasters (objects of class `SpatRaster` from the **terra** package, saved as GeoTIFF files): -#' * [madChelsa]: Bioclimatic variables -#' * [madCover]: Land cover -#' * [`madElev`][madElev]: Elevation -#' * [madForest2000]: Forest cover in year 2000 -#' * [madForest2014]: Forest cover in year 2014 -#' * [madLANDSAT]: Surface reflectance in 2023 -#' * [madPpt], [madTmin], [madTmax]: Rasters of mean monthly precipitation, and minimum and maximum temperature. +#' * [madChelsa][madChelsa]: Bioclimatic variables +#' * [madCover][madCover]: Land cover +#' * [madElev][madElev]: Elevation +#' * [madForest2000][madForest2000]: Forest cover in year 2000 +#' * [madForest2014][madForest2014]: Forest cover in year 2014 +#' * [madLANDSAT][madLANDSAT]: Surface reflectance in 2023 +#' * [madPpt][madPpt], [madTmin][madTmin], [madTmax][madTmax]: Rasters of mean monthly precipitation, and minimum and maximum temperature. #' #' Data frames -#' * [appFunsTable]: Table of functions usable by [app()] -#' * [madCoverCats]: Land cover values and categories for [madCover] -#' * [vegIndices]: Vegetation indices that can be calculated with [vegIndex()] +#' * [appFunsTable][appFunsTable]: Table of functions usable by [app()] +#' * [madCoverCats][madCoverCats]: Land cover values and categories for `madCover` +#' * [vegIndices][vegIndices]: Vegetation indices that can be calculated with [vegIndex()] #' #' @return A `SpatRaster`, `sf` spatial vector, or a `data.frame`. #' diff --git a/R/fasterRaster.r b/R/fasterRaster.r index f8754c95..a3df5d1d 100644 --- a/R/fasterRaster.r +++ b/R/fasterRaster.r @@ -322,11 +322,11 @@ #' * Vignette on **fasterRaster** hidden functions: `vignette("hidden_functions", package = "fasterRaster")` #' #' ## Classes -#' * [`GLocation`]: Fundamental class; points to a "location/project" in **GRASS** -#' * [`GSpatial`]: Basic class of any spatial object -#' * [`GRegion`]: Points to a "region" of a "location/project" in **GRASS** -#' * [`GRaster`]: Raster class -#' * [`GVector`]: Spatial vector class +#' * [GLocation][GLocation]: Fundamental class; points to a "location/project" in **GRASS** +#' * [GSpatial][GSpatial]: Basic class of any spatial object +#' * [GRegion][GRegion]: Points to a "region" of a "location/project" in **GRASS** +#' * [GRaster][GRaster]: Raster class +#' * [GVector][GVector]: Spatial vector class #' #' @author Adam B. Smith #' @name fasterRaster diff --git a/R/fillNAs.r b/R/fillNAs.r index 71269423..5881ad76 100644 --- a/R/fillNAs.r +++ b/R/fillNAs.r @@ -121,4 +121,4 @@ methods::setMethod( .makeGRaster(srcs, names(x)) } # EOF -) \ No newline at end of file +) diff --git a/R/fragmentation.r b/R/fragmentation.r index e351784f..91e8c365 100644 --- a/R/fragmentation.r +++ b/R/fragmentation.r @@ -13,9 +13,9 @@ #' * 5: Undetermined (not possible to obtain when `w = 3`) #' * 6: Interior #' -#' @param x A `SpatRaster` or `GRaster`. +#' @param x A `SpatRaster` or `GRaster` wherein the habitat type has cell values of 1, and non-focal habitat type has cell values of 0 or `NA` #' -#' @param w An odd, positive integer: Size of the window across which fragmentation is calculated (in units of "rows" and "columns"). The default is 3, meaning the function uses a 3x3 moving window to calculate fragmentation. For large rasters, compute time is ~*O*(`N`) + *O*(`N * w^2`), where `N` is the number of cells in the raster. So, even a small increase in `w` can increase compute time by a lot. +#' @param w An odd integer between 3 and 25: Size of the window across which fragmentation is calculated (in units of "rows" and "columns"). The default is 3, meaning the function uses a 3x3 moving window to calculate fragmentation. For large rasters, compute time is ~*O*(`N`) + *O*(`N * w^2`), where `N` is the number of cells in the raster. So, even a small increase in `w` can increase compute time by a lot. Values >25 will fail. A workaround might be to use [aggregate()] with `fun = 'mode'`, apply `fragmentation()`, then `resample()` to the original resolution. #' #' @param undet Character: How to assign the "undetermined" case. Valid values are `"perforated"` (default), `"edge"`, and `"undetermined"`. Partial matching is used. If `Pf` is the proportional density raster cell value and `Pff` the proportional connectivity raster cell value, the undetermined case occurs when `Pf` > 0.6 and `Pf == Pff`. #' @@ -42,7 +42,7 @@ methods::setMethod( function(x, w = 3, undet = "undetermined", none = NA, na.rm = TRUE, cores = faster("cores"), verbose = TRUE) { # errors? - if (!omnibus::is.wholeNumber(w) || w < 1L || w %% 2 == 0) stop("Argument `w` must be an odd, positive integer.") + if (w > 25 | !omnibus::is.wholeNumber(w) | w < 1L | w %% 2 == 0) stop("Argument `w` must be an odd, positive integer in the range from 3 to 25.") undet <- omnibus::pmatchSafe(undet, c("undetermined", "perforated", "edge"), nmax = 1L) # setup @@ -307,7 +307,7 @@ methods::setMethod( } # for each layer - nSteps <- nLayers * (1 + length(offsets) + 6) + nSteps <- nLayers * (2 + length(offsets) + 6) + 1 if (verbose | faster("verbose")) { pb <- utils::txtProgressBar(min = 0, max = nSteps, initial = 0, style = 3, width = 30) } @@ -330,13 +330,14 @@ methods::setMethod( created <- c(created, srcXZeros) ### connectivities (for Pff) - srcConnectivities <- .makeSourceName("fragmentation_connectivities", "raster", n = length(offsets)) - srcEither <- .makeSourceName("fragmentation_either", "raster", n = length(offsets)) - srcBoth <- .makeSourceName("fragmentation_both", n = length(offsets)) + n <- length(offsets) + srcConnectivities <- paste0('c', 1:n) # .makeSourceName("frag_connect", "raster", n = length(offsets)) + srcEither <- paste0('e', 1:n) # .makeSourceName("frag_either", "raster", n = length(offsets)) + srcBoth <- paste0('b', 1:n) # .makeSourceName("frag_both", n = length(offsets)) for (j in seq_along(offsets)) { steps <- steps + 1 - if (verbose | faster("verbose")) utils::setTxtProgressBar(pb, steps) + if (verbose || faster("verbose")) utils::setTxtProgressBar(pb, steps) # tally occurrences of habitat in each pair of cells y1 <- offsets[[j]][1L] @@ -346,26 +347,25 @@ methods::setMethod( x2 <- offsets[[j]][4L] # sum of pairs of cells - ex <- paste0(srcConnectivities[j], " = ", srcXZeros, "[", y1, ",", x1, "] + ", srcXZeros, "[", y2, ",", x2, "]") + ex <- paste0(srcConnectivities[j], " = ", srcXZeros, "[", y1, ",", x1, "]+", srcXZeros, "[", y2, ",", x2, "]") rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) # does either cell in a pair have habitat? - ex <- paste0(srcEither[j], " = if(", srcConnectivities[j], " >= 1)") + ex <- paste0(srcEither[j], " = if(", srcConnectivities[j], ">= 1)") rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) # do both cells have habitat? - ex <- paste0(srcBoth[j], " = if(", srcConnectivities[j], " == 2)") + ex <- paste0(srcBoth[j], " = if(", srcConnectivities[j], "==2)") rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) } created <- c(created, srcConnectivities, srcEither, srcBoth) - # number of neighbor cell pairs with at least one with habitat steps <- steps + 1 if (verbose | faster("verbose")) utils::setTxtProgressBar(pb, steps) srcNumEithers <- .makeSourceName("fragmentation_sum_eithers", "raster") - ex <- paste0(srcNumEithers, " = ", paste(srcEither, collapse = " + ")) + ex <- paste0(srcNumEithers, " = ", paste(srcEither, collapse = "+")) rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) created <- c(created, srcNumEithers) @@ -486,13 +486,17 @@ methods::setMethod( # } - srcs[i] <- src - + steps <- steps + 1 + if (verbose | faster("verbose")) utils::setTxtProgressBar(pb, steps) .rm(created, type = "raster", warn = FALSE) created <- character() + srcs[i] <- src + } # next layer + steps <- steps + 1 + if (verbose | faster("verbose")) utils::setTxtProgressBar(pb, steps) if (verbose | faster("verbose")) close(pb) levs <- .fragmentationLevels(undet = undet, none = none) diff --git a/R/geomorphons.r b/R/geomorphons.r index 56a68e39..4707d203 100644 --- a/R/geomorphons.r +++ b/R/geomorphons.r @@ -50,13 +50,15 @@ methods::setMethod( ) { unit <- omnibus::pmatchSafe(unit, c("cells", "meters"), nmax = 1L) - comparison <- omnibus::pmatchSafe(mode, c("1", "2", "2d"), nmax = 1L) - comparison <-if (comparison == "1") { - "anglev1" - } else if (comparison == "2") { - "anglev2" - } else if (comparison == "2d") { - "anglev2_distance" + mode <- as.character(mode) + if (mode == "1") { + comparison <- "anglev1" + } else if (mode == "2") { + comparison <- "anglev2" + } else if (mode == "2d") { + comparison <- "anglev2_distance" + } else { + stop("Argument `mode` can only be `1`, `2`, or `2d`.") } if (nlyr(x) > 1L) stop("This function requires `x` to have just one layer.") diff --git a/R/ls.r b/R/ls.r index be21dd6b..7c60c929 100644 --- a/R/ls.r +++ b/R/ls.r @@ -14,54 +14,41 @@ .ls <- function(type = c("rasters", "vectors", "rasters3d", "groups")) { rov <- c("rasters", "vectors", "rasters3d", "groups") - type <- omnibus::pmatchSafe(type, rov, useFirst = TRUE, n = length(rov)) if (is.null(type)) type <- rov + type <- omnibus::pmatchSafe(type, rov, useFirst = TRUE, n = length(rov)) # find rasters and vector rasts <- vects <- rasts3d <- vects3d <- groups <- character() - match <- pmatch(type, rov) - # failure to match - if (anyNA(match)) { - - warning("No matches.", immediate. = TRUE) - out <- NULL - - # specific match - } else { - - # rasters - if (any(match == 1L)) { - rasts <- rgrass::execGRASS("g.list", flags = .quiet(), type = "raster", intern = TRUE, echoCmd = FALSE) - if (length(rasts) > 0L) names(rasts) <- rep("raster", length(rasts)) - rasts <- sort(rasts) - } - - # 3D rasters - if (any(match == 3L)) { - rasts3d <- rgrass::execGRASS("g.list", flags=.quiet(), type="raster_3d", intern=TRUE, echoCmd=FALSE) - if (length(rasts3d) > 0L) names(rasts3d) <- rep("raster3d", length(rasts3d)) - rasts3d <- sort(rasts3d) - } - - # vectors - if (any(match == 2L)) { - vects <- rgrass::execGRASS("g.list", flags=.quiet(), type="vector", intern=TRUE, echoCmd=FALSE) - if (length(vects) > 0L) names(vects) <- rep("vector", length(vects)) - vects <- sort(vects) - } - - # groups - if (any(match == 4L)) { - groups <- rgrass::execGRASS("g.list", flags=.quiet(), type="group", intern=TRUE, echoCmd=FALSE) - if (length(groups) > 0L) names(groups) <- rep("group", length(groups)) - groups <- sort(groups) - } - - out <- c(rasts, rasts3d, vects, groups) - + # rasters + if (any(type == "rasters")) { + rasts <- rgrass::execGRASS("g.list", flags = .quiet(), type = "raster", intern = TRUE, echoCmd = FALSE) + if (length(rasts) > 0L) names(rasts) <- rep("raster", length(rasts)) + rasts <- sort(rasts) + } + + # 3D rasters + if (any(type == "rasters3d")) { + rasts3d <- rgrass::execGRASS("g.list", flags=.quiet(), type="raster_3d", intern=TRUE, echoCmd=FALSE) + if (length(rasts3d) > 0L) names(rasts3d) <- rep("raster3d", length(rasts3d)) + rasts3d <- sort(rasts3d) } - + + # vectors + if (any(type == "vectors")) { + vects <- rgrass::execGRASS("g.list", flags=.quiet(), type="vector", intern=TRUE, echoCmd=FALSE) + if (length(vects) > 0L) names(vects) <- rep("vector", length(vects)) + vects <- sort(vects) + } + + # groups + if (any(type == "groups")) { + groups <- rgrass::execGRASS("g.list", flags=.quiet(), type="group", intern=TRUE, echoCmd=FALSE) + if (length(groups) > 0L) names(groups) <- rep("group", length(groups)) + groups <- sort(groups) + } + + out <- c(rasts, rasts3d, vects, groups) out } diff --git a/R/mask.r b/R/mask.r index 5fcccab4..67e0b78b 100644 --- a/R/mask.r +++ b/R/mask.r @@ -125,12 +125,15 @@ methods::setMethod( args$raster <- maskSrc ### mask is a vector + #################### } else if (maskType == "vector") { + args <- list( cmd = "r.mask", flags = c(.quiet(), "overwrite") ) args$vector <- mask + } if (inverse) args$flags <- c(args$flags, "i") diff --git a/R/rbind.r b/R/rbind.r index de1539b0..a002dc5c 100644 --- a/R/rbind.r +++ b/R/rbind.r @@ -80,11 +80,20 @@ methods::setMethod( # } - tables <- lapply(dots, as.data.table) - table <- tryCatch( - do.call(rbind, tables), - error = function(cond) FALSE - ) + ### combine data tables if possible + tables <- lapply(dots, as.data.frame) + tables <- lapply(tables, data.table::as.data.table) + + nRows <- sapply(tables, nrow) + nGeoms <- sapply(dots, ngeom) + if (all(nRows == nGeoms)) { + table <- tryCatch( + do.call(rbind, tables), + error = function(cond) FALSE + ) + } else { + table <- FALSE + } if (is.logical(table)) { diff --git a/R/ruggedness.r b/R/ruggedness.r index ae901491..dfd67b7c 100644 --- a/R/ruggedness.r +++ b/R/ruggedness.r @@ -6,7 +6,7 @@ #' #' @param x A `GRaster`. #' @param size Integer (default is 3): Size of the moving window. Must be an odd integer >= 3. -#' @param exponent Numeric >= 0 and <= 4. Used to reduce the influence of cells farther from the focal cell (larger areas can yield noisier results if the exponent small). All cells are weighted equally when `exponent = 0`. +#' @param exponent Numeric >= 0 and <= 4. Used to reduce the influence of cells farther from the focal cell (larger areas can yield noisier results if the exponent is small). All cells are weighted equally when `exponent = 0`. #' #' @returns A `GRaster`. #' diff --git a/R/scale.r b/R/scale.r index 6a8d5c01..97f338ae 100644 --- a/R/scale.r +++ b/R/scale.r @@ -7,11 +7,15 @@ #' @param x A `GRaster`. #' #' @param center Value depends on the function: -#' * `scale()`: Logical: If `TRUE` (default), subtract from each raster layer its mean. +#' * `scale()`: +#' * Logical: If `TRUE` (default), subtract from each raster layer its mean. If `FALSE`, do not. +#' * Numeric: A single value, in which case the same value will be used across all layers of `x`, or one value per layer in `x`. #' * `unscale()`: Numeric vector or `NULL` (default): This can be a single value, which will be recycled if there is more than one layer in the raster, or one value per raster layer. If a value is `NA`, then no un-centering will be performed on the relevant raster layer. If `NULL`, then no un-centering is done. #' #' @param scale Value depends on the function: -#' * `scale()`: Logical: If `TRUE` (default), divide each layer by its standard deviation. +#' * `scale()`: +#' * Logical: If `TRUE` (default), divide each raster layer by its standard deviation. If `FALSE`, do not. +#' * Numeric: A single value, in which case the same value will be used across all layers of `x`, or one value per layer in `x`. #' * `unscale()`: Numeric vector or `NULL` (default): This can be a single value, which will be recycled if there is more than one layer in the raster, or one value per raster layer. If a value is `NA`, then no unscaling will be done on the relevant raster layer. If `NULL`, then no un-scaling is done. #' #' @returns All functions return a `GRaster`. The output of `scale()` and `scalepop()` will have two attributes, "center" and "scale", which have the means and standard deviations of the original rasters (if `center` and `scale` are `TRUE`, otherwise, they will be `NA`). These can be obtained using `attributes(output_raster)$center` and `attributes(output_raster)$scale`. @@ -53,67 +57,214 @@ methods::setMethod( #' @noRd .scale <- function(x, center, scale, sample) { - if (!center & !scale) { - - warning("No scaling performed because neither `center` nor `scale` are TRUE.") - return(x) - - } - .locationRestore(x) .region(x) - if (center) { - fx <- "mean" - } else { - fx <- NULL - } + nLayers <- nlyr(x) + + #### calculate centers and scales then center and scale + ####################################################### + if (is.logical(center) & is.logical(scale)) { - if (scale) { - if (sample) { - sdfx <- "sd" - } else{ - sdfx <- "sdpop" + if (!center & !scale) { + warning("No scaling performed because neither `center` nor `scale` are TRUE.") + return(x) + } + + if (center) { + fx <- "mean" + } else { + fx <- NULL } - } else { - sdfx <- NULL - } - fx <- c(fx, sdfx) - stats <- global(x, fx) + if (scale) { + if (sample) { + sdfx <- "sd" + } else{ + sdfx <- "sdpop" + } + } else { + sdfx <- NULL + } + fx <- c(fx, sdfx) - nLayers <- nlyr(x) - srcs <- .makeSourceName("scale_r_mapcalc", "raster", n = nLayers) - for (i in seq_len(nLayers)) { + stats <- global(x, fx) + + srcs <- .makeSourceName("scale_r_mapcalc", "raster", n = nLayers) + for (i in seq_len(nLayers)) { + + if (center) mu <- stats[i, "mean"] + if (scale) sigma <- stats[i, sdfx] + + if (center & scale) { + ex <- paste0(srcs[i], " = (", sources(x)[i], " - ", mu, ") / ", sigma) + } else if (center & !scale) { + ex <- paste0(srcs[i], " = ", sources(x)[i], " - ", mu) + } else if (!center & scale) { + ex <- paste0(srcs[i], " = ", sources(x)[i], " / ", sigma) + } + + rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) + + } # next layer + out <- .makeGRaster(srcs, names(x)) + + if (center) { + vals <- stats[ , "mean"] + names(vals) <- names(x) + attr(out, "center") <- vals + } else { + attr(out, "center") <- NA_real_ + } + + if (scale) { + vals <- stats[ , sdfx] + names(vals) <- names(x) + attr(out, "scale") <- vals + } else { + attr(out, "scale") <- NA_real_ + } + + ### user provides centers but not scales + ######################################## + } else if (!is.logical(center) & is.logical(scale)) { - if (center) mu <- stats[i, "mean"] - if (scale) sigma <- stats[i, sdfx] + len <- length(center) + if (len != nLayers) { + if (len == 1) { + center <- rep(center, nLayers) + warning("Using the same center for all rasters.") + } else { + stop("Argument `center` must be TRUE, FALSE, a single numeric value, or have the same number of numeric values as `x` has layers.") + } + } + + if (scale) { + if (sample) { + sdfx <- "sd" + } else{ + sdfx <- "sdpop" + } + stats <- global(x, sdfx) + } - if (center & scale) { - ex <- paste0(srcs[i], " = (", sources(x)[i], " - ", mu, ") / ", sigma) - } else if (center & !scale) { - ex <- paste0(srcs[i], " = ", sources(x)[i], " - ", mu) - } else if (!center & scale) { - ex <- paste0(srcs[i], " = ", sources(x)[i], " / ", sigma) + nLayers <- nlyr(x) + srcs <- .makeSourceName("scale_r_mapcalc", "raster", n = nLayers) + for (i in seq_len(nLayers)) { + + mu <- center[i] + if (scale) sigma <- stats[i, sdfx] + + if (scale) { + ex <- paste0(srcs[i], " = (", sources(x)[i], " - ", mu, ") / ", sigma) + } else if (!scale) { + ex <- paste0(srcs[i], " = ", sources(x)[i], " - ", mu) + } + + rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) + + } # next layer + out <- .makeGRaster(srcs, names(x)) + + names(center) <- names(x) + attr(out, "center") <- center + + if (scale) { + vals <- stats[ , sdfx] + names(vals) <- names(x) + attr(out, "scale") <- vals + } else { + attr(out, "scale") <- NA_real_ + } + + ### calculate centers but user supplied scales + ############################################## + } else if (is.logical(center) & !is.logical(scale)) { + + len <- length(scale) + if (len != nLayers) { + if (len == 1) { + scale <- rep(scale, len) + warning("Using the same scale for all rasters.") + } else { + stop("Argument `scale` must be TRUE, FALSE, a single numeric value, or have the same number of numeric values as `x` has layers.") + } + } + + if (center) { + fx <- "mean" + stats <- global(x, fx) + } + + srcs <- .makeSourceName("scale_r_mapcalc", "raster", n = nLayers) + for (i in seq_len(nLayers)) { + + if (center) mu <- stats[i, fx] + sigma <- scale[i] + + if (center) { + ex <- paste0(srcs[i], " = (", sources(x)[i], " - ", mu, ") / ", sigma) + } else if (!center) { + ex <- paste0(srcs[i], " = ", sources(x)[i], " / ", sigma) + } + + rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) + + } # next layer + + if (center) { + vals <- stats[ , "mean"] + names(vals) <- names(x) + attr(out, "center") <- vals + } else { + attr(out, "center") <- NA_real_ } - rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) + names(scale) <- names(x) + attr(out, "scale") <- scale + + ### user supplied centers and scales + #################################### + } else if (!is.logical(center) & !is.logical(scale)) { - } # next layer - out <- .makeGRaster(srcs, names(x)) + len <- length(center) + if (len != nLayers) { + if (len == 1) { + center <- rep(center, nLayers) + warning("Using the same center for all rasters.") + } else { + stop("Argument `center` must be TRUE, FALSE, a single numeric value, or have the same number of numeric values as `x` has layers.") + } + } - if (center) { - attr(out, "center") <- stats[ , "mean"] - } else { - attr(out, "center") <- NA_real_ - } + len <- length(scale) + if (len != nLayers) { + if (len == 1) { + scale <- rep(scale, nLayers) + warning("Using the same scale for all rasters.") + } else { + stop("Argument `scale` must be TRUE, FALSE, a single numeric value, or have the same number of numeric values as `x` has layers.") + } + } - if (scale) { - attr(out, "scale") <- stats[ , sdfx] - } else { - attr(out, "scale") <- NA_real_ - } + srcs <- .makeSourceName("scale_r_mapcalc", "raster", n = nLayers) + for (i in seq_len(nLayers)) { + + mu <- center[i] + sigma <- scale[i] + ex <- paste0(srcs[i], " = (", sources(x)[i], " - ", mu, ") / ", sigma) + rgrass::execGRASS("r.mapcalc", expression = ex, flags = c(.quiet(), "overwrite")) + + } # next layer + out <- .makeGRaster(srcs, names(x)) + + names(center) <- names(x) + names(scale) <- names(x) + attr(out, "center") <- center + attr(out, "scale") <- scale + + } out } # EOF diff --git a/R/sun.r b/R/sun.r index 32381211..6d060840 100644 --- a/R/sun.r +++ b/R/sun.r @@ -8,13 +8,13 @@ #' #' @param coeff_dh A `GRaster`: A raster with values of the real-sky diffuse radiation coefficient. Valid values are between 0 and 1. #' -#' @param slope A `GRaster`: This is a raster representing topographic slope in radians. It can be generated using [terrain()]. +#' @param slope A `GRaster`: This is a raster representing topographic slope in degrees. It can be generated using [terrain()]. #' #' @param aspect A `GRaster`: This is a raster representing topographic aspect in degrees. It can be generated using [terrain()]. If generated with that function, "east orientation" *must* be used (i.e., argument `northIs0` must be `FALSE`). #' #' @param hh A "stack" of `GRaster`s: This represents height of the horizon in radians in particular directions. Horizon height can be calculated using [horizonHeight()]. The directions *must* be in "east orientation" (i.e., argument `northIs0` in `horzionHeight()` must be `FALSE`). The directions must correspond with the sequence given by `horizon_step` (see next argument). For example, if `horizon_step` is 90, then `hh` must contain rasters representing horizon height at 0 (east), 90 (north), 180 (west), and 270 (south) aspects. -#'#' -#' @param horizon_step Numeric >0: Difference between angular steps in which horizon height is measured. One horizon height raster will be made per value from 0 to 360 - `horizon_step` degrees. +#' +#' @param horizon_step Numeric >0: Difference between angular steps in which horizon height is measured. One horizon height raster will be made per value from 0 to (360 - `horizon_step`) degrees. #' #' @param albedo A `GRaster` or a numeric value: This is either a raster with values of ground albedo or a numeric value (in which case albedo is assumed to be the same everywhere). Albedo is unit-less, and the default value is 0.2. #' diff --git a/README.md b/README.md index 0b6d4d95..bb8b2609 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,12 @@ # fasterRaster -[![cran version](https://www.r-pkg.org/badges/version/fasterRaster)](https://cran.r-project.org/package=fasterRaster) [![R-CMD-check](https://github.com/adamlilith/fasterRaster/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/adamlilith/fasterRaster/actions/workflows/R-CMD-check.yaml) [![GPLv3 license](https://img.shields.io/badge/License-GPLv3-blue.svg)](http://perso.crans.org/besson/LICENSE.html) +[![cran version](https://www.r-pkg.org/badges/version/fasterRaster)](https://cran.r-project.org/package=fasterRaster) [![R-CMD-check](https://github.com/adamlilith/fasterRaster/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/adamlilith/fasterRaster/actions/workflows/R-CMD-check.yaml) [![GPLv3 license](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://perso.crans.org/besson/LICENSE.html) Faster raster processing in `R` using `GRASS` -fasterRaster website +fasterRaster website `fasterRaster` is an **R** package designed specifically to handle large-in-memory/large-on-disk spatial rasters and vectors. `fasterRaster` does this using Open Source Geospatial's `GRASS` @@ -189,8 +189,8 @@ Note that the `M1.M2` and `S1.S2` increment independently. For example, if the v * The Wiki on [how to run `GRASS` in `R` or `R` in `GRASS`](https://grasswiki.osgeo.org/wiki/R_statistics/rgrass) will help you to become a power-user of `GRASS` in `R`. # Citation -A publication is forthcoming! In the meantime, please cite: +A publication is forthcoming! In the meantime, please see and cite: -Smith, A.B. 2024. `fasterRaster`: Faster raster processing in `R` using `GRASS`. URL: https://cran.r-project.org/package=fasterRaster. +Smith, A.B. 2025. `fasterRaster`: GIS in `R` using `GRASS` for large vectors and rasters. EarthArXiv preprint doi: 10.31223/X52R0M ~ Adam diff --git a/data/madCoast0.rda b/data/madCoast0.rda index 609089b0..9c06390b 100644 Binary files a/data/madCoast0.rda and b/data/madCoast0.rda differ diff --git a/data/madCoast4.rda b/data/madCoast4.rda index f53bfc72..a6824d13 100644 Binary files a/data/madCoast4.rda and b/data/madCoast4.rda differ diff --git a/data/madDypsis.rda b/data/madDypsis.rda index b352405e..052ec798 100644 Binary files a/data/madDypsis.rda and b/data/madDypsis.rda differ diff --git a/data/madRivers.rda b/data/madRivers.rda index 9e1cdad5..ad437a81 100644 Binary files a/data/madRivers.rda and b/data/madRivers.rda differ diff --git a/inst/pkgdown.yml b/inst/pkgdown.yml index a0dc858b..5bbdd6d7 100644 --- a/inst/pkgdown.yml +++ b/inst/pkgdown.yml @@ -1,5 +1,5 @@ pandoc: 3.6.2 -pkgdown: 2.1.3 +pkgdown: 2.2.0 pkgdown_sha: ~ articles: faster_fasterRaster: faster_fasterRaster.html @@ -9,7 +9,7 @@ articles: projects_mapsets: projects_mapsets.html regions: regions.html three_d_objects: three_d_objects.html -last_built: 2025-06-19T21:54Z +last_built: 2025-11-18T15:29Z urls: reference: https://github.com/adamlilith/fasterRaster/reference article: https://github.com/adamlilith/fasterRaster/articles diff --git a/man/Arithmetic.Rd b/man/Arithmetic.Rd index 988d3d5d..e052519e 100644 --- a/man/Arithmetic.Rd +++ b/man/Arithmetic.Rd @@ -29,7 +29,7 @@ \S4method{Arith}{GVector,GVector}(e1, e2) } \arguments{ -\item{e1, e2}{\code{GRaster}s, \code{numeric}s, \code{integer}s, or \code{logical}s.} +\item{e1, e2}{\code{GRaster}s, \code{GVector}s, \code{numeric}s, \code{integer}s, or \code{logical}s.} } \value{ A \code{GRaster}. diff --git a/man/dot-backdoor.Rd b/man/dot-backdoor.Rd index ad34c679..9e66417e 100644 --- a/man/dot-backdoor.Rd +++ b/man/dot-backdoor.Rd @@ -4,9 +4,11 @@ \alias{.backdoor} \title{Setup fasterRaster for ABS} \usage{ -.backdoor(ver = "84") +.backdoor(start = FALSE, ver = "84") } \arguments{ +\item{start}{Logical: If \code{TRUE}, start the \strong{GRASS} session by creating the \code{madElev} \code{GRaster}.} + \item{ver}{Character: \strong{GRASS}: e.g., "83" or "84".} } \value{ diff --git a/man/fastData.Rd b/man/fastData.Rd index f41dbdbc..a0900e7b 100644 --- a/man/fastData.Rd +++ b/man/fastData.Rd @@ -21,7 +21,7 @@ Rasters (objects of class \code{SpatRaster} from the \strong{terra} package, sav \itemize{ \item \link{madChelsa}: Bioclimatic variables \item \link{madCover}: Land cover -\item \code{\link{madElev}}: Elevation +\item \link{madElev}: Elevation \item \link{madForest2000}: Forest cover in year 2000 \item \link{madForest2014}: Forest cover in year 2014 \item \link{madLANDSAT}: Surface reflectance in 2023 @@ -31,7 +31,7 @@ Rasters (objects of class \code{SpatRaster} from the \strong{terra} package, sav Data frames \itemize{ \item \link{appFunsTable}: Table of functions usable by \code{\link[=app]{app()}} -\item \link{madCoverCats}: Land cover values and categories for \link{madCover} +\item \link{madCoverCats}: Land cover values and categories for \code{madCover} \item \link{vegIndices}: Vegetation indices that can be calculated with \code{\link[=vegIndex]{vegIndex()}} }} } diff --git a/man/fasterRaster.Rd b/man/fasterRaster.Rd index 13bee210..fd5ffa70 100644 --- a/man/fasterRaster.Rd +++ b/man/fasterRaster.Rd @@ -388,11 +388,11 @@ Operations on \code{GRaster}s \subsection{Classes}{ \itemize{ -\item \code{\link{GLocation}}: Fundamental class; points to a "location/project" in \strong{GRASS} -\item \code{\link{GSpatial}}: Basic class of any spatial object -\item \code{\link{GRegion}}: Points to a "region" of a "location/project" in \strong{GRASS} -\item \code{\link{GRaster}}: Raster class -\item \code{\link{GVector}}: Spatial vector class +\item \link{GLocation}: Fundamental class; points to a "location/project" in \strong{GRASS} +\item \link{GSpatial}: Basic class of any spatial object +\item \link{GRegion}: Points to a "region" of a "location/project" in \strong{GRASS} +\item \link{GRaster}: Raster class +\item \link{GVector}: Spatial vector class } } } diff --git a/man/figures/logo.png b/man/figures/logo.png index b9e42114..bf0f14c0 100644 Binary files a/man/figures/logo.png and b/man/figures/logo.png differ diff --git a/man/fragmentation.Rd b/man/fragmentation.Rd index ebc26a02..147a01ef 100644 --- a/man/fragmentation.Rd +++ b/man/fragmentation.Rd @@ -19,9 +19,9 @@ \S4method{fragmentation}{GRaster}(x, w = 3, undet = "undetermined", none = NA, verbose = TRUE) } \arguments{ -\item{x}{A \code{SpatRaster} or \code{GRaster}.} +\item{x}{A \code{SpatRaster} or \code{GRaster} wherein the habitat type has cell values of 1, and non-focal habitat type has cell values of 0 or \code{NA}} -\item{w}{An odd, positive integer: Size of the window across which fragmentation is calculated (in units of "rows" and "columns"). The default is 3, meaning the function uses a 3x3 moving window to calculate fragmentation. For large rasters, compute time is ~\emph{O}(\code{N}) + \emph{O}(\code{N * w^2}), where \code{N} is the number of cells in the raster. So, even a small increase in \code{w} can increase compute time by a lot.} +\item{w}{An odd integer between 3 and 25: Size of the window across which fragmentation is calculated (in units of "rows" and "columns"). The default is 3, meaning the function uses a 3x3 moving window to calculate fragmentation. For large rasters, compute time is ~\emph{O}(\code{N}) + \emph{O}(\code{N * w^2}), where \code{N} is the number of cells in the raster. So, even a small increase in \code{w} can increase compute time by a lot. Values >25 will fail. A workaround might be to use \code{\link[=aggregate]{aggregate()}} with \code{fun = 'mode'}, apply \code{fragmentation()}, then \code{resample()} to the original resolution.} \item{undet}{Character: How to assign the "undetermined" case. Valid values are \code{"perforated"} (default), \code{"edge"}, and \code{"undetermined"}. Partial matching is used. If \code{Pf} is the proportional density raster cell value and \code{Pff} the proportional connectivity raster cell value, the undetermined case occurs when \code{Pf} > 0.6 and \code{Pf == Pff}.} diff --git a/man/hist.Rd b/man/hist.Rd index 9c5c25a9..4c95648a 100644 --- a/man/hist.Rd +++ b/man/hist.Rd @@ -10,7 +10,7 @@ \arguments{ \item{x}{A \code{GRaster}.} -\item{layer}{Character, numeric, or integer: Indicates which layer of a multi-layer \code{GRaster} for which to plot a histogram. The layer can be identified using its \code{\link[=name]{name()}} (character) or index (numeric or integer). If this is missing, then up to \code{maxnl} layers are plotted.} +\item{layer}{Character, numeric, or integer: Indicates which layer of a multi-layer \code{GRaster} for which to plot a histogram. The layer can be identified using its \code{\link[terra:names]{terra::name()}} (character) or index (numeric or integer). If this is missing, then up to \code{maxnl} layers are plotted.} \item{maxnl}{Maximum number of layers for which to create histograms. This is 16 by default, but ignored if \code{layer} is defined.} diff --git a/man/ruggedness.Rd b/man/ruggedness.Rd index 6969e063..d49a0f83 100644 --- a/man/ruggedness.Rd +++ b/man/ruggedness.Rd @@ -12,7 +12,7 @@ \item{size}{Integer (default is 3): Size of the moving window. Must be an odd integer >= 3.} -\item{exponent}{Numeric >= 0 and <= 4. Used to reduce the influence of cells farther from the focal cell (larger areas can yield noisier results if the exponent small). All cells are weighted equally when \code{exponent = 0}.} +\item{exponent}{Numeric >= 0 and <= 4. Used to reduce the influence of cells farther from the focal cell (larger areas can yield noisier results if the exponent is small). All cells are weighted equally when \code{exponent = 0}.} } \value{ A \code{GRaster}. diff --git a/man/scale.Rd b/man/scale.Rd index 17b0fc77..82b7aa08 100644 --- a/man/scale.Rd +++ b/man/scale.Rd @@ -20,13 +20,21 @@ \item{center}{Value depends on the function: \itemize{ -\item \code{scale()}: Logical: If \code{TRUE} (default), subtract from each raster layer its mean. +\item \code{scale()}: +\itemize{ +\item Logical: If \code{TRUE} (default), subtract from each raster layer its mean. If \code{FALSE}, do not. +\item Numeric: A single value, in which case the same value will be used across all layers of \code{x}, or one value per layer in \code{x}. +} \item \code{unscale()}: Numeric vector or \code{NULL} (default): This can be a single value, which will be recycled if there is more than one layer in the raster, or one value per raster layer. If a value is \code{NA}, then no un-centering will be performed on the relevant raster layer. If \code{NULL}, then no un-centering is done. }} \item{scale}{Value depends on the function: \itemize{ -\item \code{scale()}: Logical: If \code{TRUE} (default), divide each layer by its standard deviation. +\item \code{scale()}: +\itemize{ +\item Logical: If \code{TRUE} (default), divide each raster layer by its standard deviation. If \code{FALSE}, do not. +\item Numeric: A single value, in which case the same value will be used across all layers of \code{x}, or one value per layer in \code{x}. +} \item \code{unscale()}: Numeric vector or \code{NULL} (default): This can be a single value, which will be recycled if there is more than one layer in the raster, or one value per raster layer. If a value is \code{NA}, then no unscaling will be done on the relevant raster layer. If \code{NULL}, then no un-scaling is done. }} } diff --git a/man/sun.Rd b/man/sun.Rd index 223ea8f6..a10d7dff 100644 --- a/man/sun.Rd +++ b/man/sun.Rd @@ -35,14 +35,13 @@ sun( \item{coeff_dh}{A \code{GRaster}: A raster with values of the real-sky diffuse radiation coefficient. Valid values are between 0 and 1.} -\item{slope}{A \code{GRaster}: This is a raster representing topographic slope in radians. It can be generated using \code{\link[=terrain]{terrain()}}.} +\item{slope}{A \code{GRaster}: This is a raster representing topographic slope in degrees. It can be generated using \code{\link[=terrain]{terrain()}}.} \item{aspect}{A \code{GRaster}: This is a raster representing topographic aspect in degrees. It can be generated using \code{\link[=terrain]{terrain()}}. If generated with that function, "east orientation" \emph{must} be used (i.e., argument \code{northIs0} must be \code{FALSE}).} -\item{hh}{A "stack" of \code{GRaster}s: This represents height of the horizon in radians in particular directions. Horizon height can be calculated using \code{\link[=horizonHeight]{horizonHeight()}}. The directions \emph{must} be in "east orientation" (i.e., argument \code{northIs0} in \code{horzionHeight()} must be \code{FALSE}). The directions must correspond with the sequence given by \code{horizon_step} (see next argument). For example, if \code{horizon_step} is 90, then \code{hh} must contain rasters representing horizon height at 0 (east), 90 (north), 180 (west), and 270 (south) aspects. -#'} +\item{hh}{A "stack" of \code{GRaster}s: This represents height of the horizon in radians in particular directions. Horizon height can be calculated using \code{\link[=horizonHeight]{horizonHeight()}}. The directions \emph{must} be in "east orientation" (i.e., argument \code{northIs0} in \code{horzionHeight()} must be \code{FALSE}). The directions must correspond with the sequence given by \code{horizon_step} (see next argument). For example, if \code{horizon_step} is 90, then \code{hh} must contain rasters representing horizon height at 0 (east), 90 (north), 180 (west), and 270 (south) aspects.} -\item{horizon_step}{Numeric >0: Difference between angular steps in which horizon height is measured. One horizon height raster will be made per value from 0 to 360 - \code{horizon_step} degrees.} +\item{horizon_step}{Numeric >0: Difference between angular steps in which horizon height is measured. One horizon height raster will be made per value from 0 to (360 - \code{horizon_step}) degrees.} \item{albedo}{A \code{GRaster} or a numeric value: This is either a raster with values of ground albedo or a numeric value (in which case albedo is assumed to be the same everywhere). Albedo is unit-less, and the default value is 0.2.} diff --git a/fasterRaster_workspace.code-workspace b/workspace_fasterRaster.code-workspace similarity index 100% rename from fasterRaster_workspace.code-workspace rename to workspace_fasterRaster.code-workspace