R/project.r
In adamlilith/fasterRaster: Faster Raster and Spatial Vector Processing Using 'GRASS GIS'
Defines functions
.regionResCenterMethod
.regionNumericRes
.projectRaster
#' Change the coordinate reference system of a GRaster or GVector
#'
#' @description `project()` changes the coordinate reference system (CRS) of a `GRaster` or `GVector`. It has three use cases:
#' * `x` is a `GRaster` and `y` is a `GRaster`: `x` will be projected to the CRS of `y` and resampled to have the same resolution as `y`. If argument `align` is `FALSE`, then it will also be cropped to the extent of `y`.
#' * `x` is a `GRaster` and `y` is a `GVector` or a CRS string (typically in Well-Known Text format): `x` will be projected to the CRS specified by `y` and resampled but not cropped.
#' * `x` is a `GVector` and `y` is a `GRaster`, `GVector`, or CRS string: The vector will be projected to the CRS of `y`.
#'
#' @param x A `GRaster` or `GVector` to be projected.
#'
#' @param y A character or `GLocation` object (i.e., typically a `GRaster` or `GVector`): Used to set the focal `GRaster` or `GVector`'s new CRS (and resolution and possibly extent, for `GRaster`s).
#'
#' @param align Logical: If `FALSE` (default), and `x` and `y` are `GRaster`s, then the extent of `x` will be cropped to the extent of `y`. If `TRUE`, no cropping is performed.
#'
#' @param method Character or `NULL` (for `GRaster`s only): Method to use to conduct the transformation (rasters only). Partial matching is used.
#' * `NULL` (default): Automatically choose based on raster properties (`near` for categorical data, `bilinear` for continuous data).
#' * `"near"`: Nearest neighbor. Best for categorical data, and often a poor choice for continuous data. If [datatype()] is `integer`, this method will be used by default.
#' * `"bilinear"`: Bilinear interpolation (default for non-categorical data; uses weighted values from 4 cells).
#' * `"bicubic"`: Bicubic interpolation (uses weighted values from 16 cells).
#' * `"lanczos"`: Lanczos interpolation (uses weighted values from 25 cells).
#'
#' *Note #1*: If `x` and `y` are `GRaster`s, and `res = "terra"`, then the same `method` is used to resample `x` to the resolution of `y` before projecting `x`.
#'
#' *Note #2*: Methods that use multiple cells will cause the focal cell to become `NA` if there is at least one cell with an `NA` in the cells it draws from. These `NA` cells can often be filled using the `fallback` argument.
#'
#' @param fallback Logical (for projecting `GRaster`s only): If `TRUE` (default), then use "lower" resampling methods to fill in `NA` cells when a "higher" resampling method is used. For example, if `method = "bicubic"`, `NA` cells will be filled in using the `bilinear` method, except when that results in `NA`s, in which case the `near` method will be used. Fallback causes fewer cells to revert to `NA` values, so may be better at capturing complex "edges" (e.g., coastlines). Fallback does increase processing time because each "lower" method must be applied, then results merged. Fallback is not used if `method = "near"`.
#'
#' @param res Character (for projecting `GRaster`s only): Method used to set the resolution of a `GRaster` in the new CRS. This can be one of three options. Partial matching is used and case ignored:
#' * `"terra"`: This method creates an output raster that is as close as possible in values and resolution to the one that [terra::project()] would create. However, for large rasters (i.e., many cells), this can fail because `terra::project()` encounters memory limits (it is used internally to create a template). This method resamples the focal raster in its starting CRS, then projects it to the destination CRS.
#' * `"template"`: This method can only be used if `y` is a `GRaster`. The output will have the same resolution as `y` and possibly the same extent (depending on the value of `align`). However, unlike the `"terra"` method, cell values will not necessarily be as close as possible to what [terra::project()] would generate (unless `method = "near"`). Unlike the `"terra"` method, this method does not resample the focal raster in its starting CRS before projecting. For large rasters it will be faster than the `"terra"` method (especially if `"method = "near"`), and it should be less likely to fail because of memory limits.
#' * Two numeric values: Values for the new resolution (x- and y-dimensions).
#' * `"center"`: This method locates the centroid of the raster to be projected (in the same CRS as the original raster). It then creates four points north, south, east, and west of the centroid, each spaced one cell's width from the centroid. This set of points is then projected to the new CRS. The new cell size in the x-dimension will be the average of the distance between the east and west points from the centroid, and in the y-dimension the average from the centroid to the north and south points.
#' * `"fallback"` (default): This applies the `terra` method first, but if that fails, then tries `template`, then `center`. This process can take a long time for large rasters.
#'
#' @param wrap Logical:
#' * `GRaster`s: When projecting rasters that "wrap around" (i.e., whole-world rasters or rasters that have edges that actually circle around to meet on the globe), `wrap` should be `TRUE` to avoid removing rows and columns from the "edge" of the map. The default is `FALSE`.
#' * `GVector`s: When projecting vectors that span the international date line at 180E/W, `wrap` should be `TRUE` to avoid an issue where the coordinates are incorrectly mapped to the range -180 to 180.
#'
#' @param verbose Logical (for projecting `GRaster`s only): If `TRUE`, display progress. Default is `FALSE`.
#'
#' @details When projecting a raster, the "fallback" methods in **GRASS** module `r.import` are actually used, even though the `method` argument takes the strings specifying non-fallback methods. See the manual page for the `r.import` **GRASS** module.
#'
#' @returns A `GRaster` or `GVector`.
#'
#' @seealso [terra::project()], [sf::st_transform()], **GRASS** manual pages for modules `r.proj` and `v.proj` (see `grassHelp("r.proj")` and `grassHelp("v.proj")`)
#'
#' @example man/examples/ex_project.r
#'
#' @aliases project
#' @rdname project
#' @exportMethod project
methods::setMethod(
f = "project",
signature = c(x = "GRaster"),
definition = function(
x,
y,
align = FALSE,
method = NULL,
fallback = TRUE,
res = "fallback",
wrap = FALSE,
verbose = FALSE
) {
if (inherits(res, "character")) {
res <- omnibus::pmatchSafe(res, c("fallback", "terra", "template", "center"), nmax = 1L)
if (!inherits(y, "GRaster") && res == "template") stop("The `template` method for setting resolution can only be used if the argument `y` is a GRaster.")
} else if (inherits(res, "numeric")) {
if (length(res) == 1L) res <- c(res, res)
if (length(res) != 2L) stop("Argument `res` must have one or two numeric values, or be a string.")
} else {
stop("Argument `res` must be a numeric vector or a string.")
}
if (inherits(res, "numeric")) {
out <- .projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = res, wrap = wrap)
} else if (res != "fallback") {
out <- .projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = res, wrap = wrap)
} else if (res == "fallback") {
out <- FALSE
j <- 0L
n <- length(res)
try <- c("terra", "template", "center")
while (inherits(out, "logical") & j <= n) {
j <- j + 1L
thisRes <- try[j]
out <- tryCatch(
.projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = thisRes, wrap = wrap, verbose = verbose),
error = function(cond) FALSE
)
}
}
out
} # EOF
)
#' @param x A `GRaster`.
#' @param y A `GRaster` or `GVector`.
#' @param align T/F
#' @param method Character
#' @param fallback T/F
#' @param res Character.
#' @param wrap T/F
#' @param verbose T/F
#' @noRd
.projectRaster <- function(x, y, align, method, fallback, res, wrap, verbose = FALSE) {
nLayers <- nlyr(x)
# .message(msg = "project_raster", message = "This function can produce erroneous results if the raster crosses a pole or the international date line.")
# imperfect catch for cases where `wrap` should be `TRUE` but is not
unproj <- .projection(x) %in%
c("Latitude-Longitude", "WGS84", "WGS 84", "NAD83", "NAD 83", "NAD27", "NAD 27")
if (unproj & !wrap) {
extent <- ext(x, vector = TRUE)
if ((extent[1L] == -180 & extent[2L] == 180) | (extent[3L] == -90 & extent[4L] == 90)) warning("This GRaster seems to wrap around the globe to meet at the international\n date line and/or the poles. Should `wrap` be `TRUE`?")
}
xLocation <- .location(x)
yLocation <- .locationFind(y, match = "crs", return = "name")
if (!is.null(yLocation)) {
if (yLocation == xLocation) {
warning("Object is already in the desired coordinate reference system.")
return()
}
} else if (is.null(yLocation)) {
yLocation <- .locationCreate(y)
yLocation <- .location(yLocation)
}
# method
if (!is.null(method)) method <- omnibus::pmatchSafe(method, c("nearest", "bilinear", "bicubic", "lanczos"))
if (is.null(method)) {
dt <- datatype(x)
method <- if (all(dt %in% c("integer", "factor"))) {
"nearest"
} else if (all(dt %in% c("float", "double"))) {
"bilinear"
} else {
stop("Rasters are a mix of datatype integer and non-integer (categorical and continuous). \nProject each type of raster separately.")
}
}
nSteps <- nLayers + (res == "terra")
if (verbose | faster("verbose")) {
pb <- utils::txtProgressBar(min = 0, max = nSteps, initial = 0, style = 3, width = 30)
steps <- 0
}
### If y is a GRaster, reshape region in target location using y's extent/resolution.
if (inherits(y, "GRaster")) {
# user sets target cell size
if (is.numeric(res)) {
.regionNumericRes(x = x, y = y, align = align, res = res)
# resample as per terra::project()
} else if (res == "terra") {
if (verbose | faster("verbose")) {
utils::setTxtProgressBar(pb, steps)
steps <- steps + 1
}
# Use a SpatRaster as template for resampling region. We do this so we can set the "region" resolution and extent correctly.
extent <- ext(x, vector = TRUE)
xRast <- matrix(NA_real_, nrow = nrow(x), ncol = ncol(x))
xRast <- terra::rast(xRast, crs = crs(x), extent = extent)
extent <- ext(y, vector = TRUE)
yRast <- matrix(NA_real_, nrow = nrow(y), ncol = ncol(y))
yRast <- terra::rast(yRast, crs = crs(y), extent = extent)
xRast <- terra::project(xRast, yRast, align = TRUE)
xRast <- terra::project(xRast, crs(x))
xSR <- xRast
xRast[] <- 1L
xRast <- fast(xRast)
# resample x in its native location to the resolution it will have in the target location
x <- resample(
x = x,
y = xRast,
method = method,
fallback = fallback
)
xRast <- terra::project(xSR, yRast, method = "near", align = align)
# reshape region
.locationRestore(yLocation)
.region(xRast)
# } else if (res == "dimensions") {
# .regionRespectsDims(x = x, y = y, align = align)
# same extent and resolution as terra, but without the resampling before projection
} else if (res == "template") {
if (!align) {
# reshape region
.locationRestore(y)
.region(y)
} else {
# extent from x projected to y
n <- N(x)
s <- S(x)
e <- E(x)
w <- W(x)
corners <- rbind(
c(w, n),
c(w, s),
c(e, s),
c(e, n)
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- terra::project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
ewres <- xres(y)
nsres <- yres(y)
yextent <- ext(y, vector = TRUE)
if (n > yextent[4L]) {
diff <- n - yextent[4L]
nc <- ceiling(diff / nsres)
n <- yextent[4L] + nc * nsres
}
if (s < yextent[3L]) {
diff <- yextent[3L] - s
nc <- ceiling(diff / nsres)
s <- yextent[3L] - nc * nsres
}
if (e > yextent[2L]) {
diff <- e - yextent[2L]
nc <- ceiling(diff / ewres)
e <- yextent[2L] + nc * ewres
}
if (w < yextent[1L]) {
diff <- yextent[1L] - w
nc <- ceiling(diff / ewres)
w <- yextent[1L] - nc * ewres
}
rows <- round((y - s) / nsres)
cols <- round((e - w) / ewres)
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
.locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
# "center" method of determining cell resolution
} else if (res == "center") {
.regionResCenterMethod(x = x, y = y, align = align)
}
### "y" is not a raster (i.e., WKT string or GVector)
} else {
# user-supplied resolution
if (inherits(res, "numeric")) {
.regionNumericRes(x = x, y = y, align = TRUE, res = res)
# project as per terra::project()
} else if (res == "terra") {
# make template raster to match raster to be projected
extent <- ext(x, vector = TRUE)
xRast <- matrix(NA_integer_, nrow = nrow(x), ncol = ncol(x))
xSR <- terra::rast(xRast, crs = crs(x), extent = extent)
xRast <- terra::project(xSR, crs(y), method = "near", align = align)
# reshape region
.locationRestore(yLocation)
.region(xRast)
# # respect dimensions
# } else if (res == "dimensions") {
# .regionRespectsDims(x = x, y = y, align = align)
# "center" method of determining cell resolution
} else if (res == "center") {
.regionResCenterMethod(x = x, y = y, align = align)
}
}
### project raster
##################
if (method != "nearest" & fallback) method <- paste0(method, "_f")
srcs <- .makeSourceName("project", "raster", nlyr(x))
for (i in seq_len(nlyr(x))) {
if (verbose | faster("verbose")) {
utils::setTxtProgressBar(pb, steps)
steps <- steps + 1
}
args <- list(
cmd = "r.proj",
mapset = .mapset(x),
input = sources(x)[i],
output = srcs[i],
method = method,
memory = faster("memory"),
flags = c(.quiet(), "overwrite")
)
args <- .addLocationProject(args, .location(x))
if (wrap) args$flags <- c(args$flags, "n")
do.call(rgrass::execGRASS, args=args)
if (is.factor(x)[i] & method == "nearest") {
levels <- cats(x)[[i]]
} else {
levels <- NULL
}
thisOut <- .makeGRaster(srcs[i], names(x)[i], levels = levels, ac = activeCat(x, layer = i))
if (i == 1L) {
out <- thisOut
} else {
out <- c(out, thisOut)
}
} # project next raster
.locationRestore(y)
# if using y as extent to which to crop
if (!align & inherits(y, "GRaster")) out <- crop(out, y)
if (verbose | faster("verbose")) close(pb)
out
}
#' @aliases project
#' @rdname project
#' @exportMethod project
methods::setMethod(
f = "project",
signature = c(x = "GVector"),
definition = function(
x,
y,
wrap = FALSE
) {
# imperfect catch for cases where `wrap` should be `TRUE` but is not
unproj <- .projection(x) %in%
c("Latitude-Longitude", "WGS84", "WGS 84", "NAD83", "NAD 83", "NAD27", "NAD 27")
if (unproj & !wrap) {
extent <- ext(x, vector = TRUE)
if (
(extent[1L] == -180 & extent[2L] == 180) |
(extent[3L] == -90 & extent[4L] == 90) |
(extent[1L] > 0 & extent[2L] < 0) |
(extent[3L] > extent[4L]) |
(extent[4L] < extent[3L])
) warning("This GVector seems to span the international date line and/or poles. Should `wrap` be `TRUE`?")
}
xLocation <- .location(x)
yLocation <- .locationFind(y, match = "crs", return = "name")
if (!is.null(yLocation)) {
if (yLocation == xLocation) {
warning("Object is already in the desired coordinate reference system.")
return(x)
}
} else if (is.null(yLocation)) {
yLocation <- .locationCreate(y)
yLocation <- .location(yLocation)
}
.locationRestore(yLocation)
src <- .makeSourceName("project_v_proj", "vector")
args <- list(
cmd = "v.proj",
dbase = faster("workDir"),
mapset = .mapset(x),
input = sources(x),
output = src,
flags = c(.quiet(), "overwrite")
)
args <- .addLocationProject(args, .location(x))
if (wrap) args$flags <- c(args$flags, "w") # if crosses international date line
# if (geomtype(x, grass = TRUE) == "point") args$flags <- c(args$flags, "b") # disable topology build for points
do.call(rgrass::execGRASS, args = args)
out <- .makeGVector(src, table = x@table)
out
} # EOF
)
#' Set region in target location so it has user-defined resolution
#'
#' @param x A `GRaster`.
#' @param y `GRaster` or `GVector`
#' @param align Logical
#'
#' @noRd
.regionNumericRes <- function(x, y, align, res) {
if (!align) {
# extent from y
n <- N(y)
s <- S(y)
e <- E(y)
w <- W(y)
} else {
# extent from x projected to y
n <- N(x)
s <- S(x)
e <- E(x)
w <- W(x)
corners <- rbind(
c(w, n),
c(w, s),
c(e, s),
c(e, n)
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- terra::project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
}
cols <- ceiling((e - w) / res[1L])
rows <- ceiling((n - s) / res[2L])
e <- w + cols * res[1L]
s <- n - rows * res[2L]
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
.locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
# #' Set region in target location so it has same dimensions (rows, columns) as focal raster
# #'
# #' @param x `GRaster` to project
# #' @param y `GRaster` or `GVector`.
# #' @param align Logical.
# #'
# #' @noRD
# .regionRespectsDims <- function(x, y, align) {
# # get new extent in y location
# if (!align) {
# extent <- ext(y, vector = TRUE)
# } else {
# # get new extent from projected corners of x raster
# n <- N(x)
# s <- S(x)
# e <- E(x)
# w <- W(x)
# corners <- rbind(
# c(w, n),
# c(w, s),
# c(e, s),
# c(e, n)
# )
# if (inherits(y, "GLocation")) {
# yCrs <- crs(y)
# } else {
# yCrs <- y
# }
# corners <- terra::vect(corners, type = "points", crs = crs(x))
# corners <- terra::project(corners, yCrs)
# corners <- terra::crds(corners)
# extent <- c(
# xmin = min(corners[1L:2L, 1L]),
# xmax = max(corners[2L:3L, 1L]),
# ymin = min(corners[2L:3L, 2L]),
# ymax = max(corners[c(1L, 4L), 2L])
# )
# }
# xdim <- dim(x)
# rows <- xdim[1L]
# cols <- xdim[2L]
# n <- as.character(extent[4L])
# s <- as.character(extent[3L])
# e <- as.character(extent[2L])
# w <- as.character(extent[1L])
# # reshape region
# .locationRestore(y)
# rgrass::execGRASS(
# cmd = "g.region",
# n = n, s = s, e = e, w = w,
# rows = rows, cols = cols,
# flags = c("o", .quiet())
# )
# }
#' Set target region resolution using "center" method
#'
#' Set target region resolution using "center" method. Extent will be from y or projected corners of x
#'
#' @param x `GRaster`.
#' @param x `GRaster` or `GVector`.
#' @param align Logical.
#'
#' @noRd
.regionResCenterMethod <- function(x, y, align) {
extent <- ext(x, vector = TRUE)
# get resolution in location of y
center <- c(mean(extent[1L:2L]), mean(extent[3L:4L]))
xres <- res(x)
wp <- c(center[1L] - xres[1L], center[2L])
ep <- c(center[1L] + xres[1L], center[2L])
np <- c(center[1L], center[2L] + xres[2L])
sp <- c(center[1L], center[2L] - xres[2L])
pts <- rbind(center, wp, ep, np, sp)
pts <- terra::vect(pts, type = "points", crs = crs(x))
pts <- terra::project(pts, crs(y))
pts <- terra::crds(pts)
centerx <- pts[1L, 1L]
centery <- pts[1L, 2L]
ewres <- mean(c(centerx - pts[2L, 1L], pts[3L, 1L] - centerx))
nsres <- mean(c(pts[4L, 2L] - centery, centery - pts[5L, 2L]))
# get extent in y location
if (!align & inherits(y, "GRegion")) {
# extent from y
yextent <- ext(y, vector = TRUE)
n <- N(y)
s <- S(y)
e <- E(y)
w <- W(y)
} else {
# extent by projecting corners of x to y location
corners <- rbind(
extent[c(1L, 4L)],
extent[c(1L, 3L)],
extent[c(2L, 3L)],
extent[c(2L, 4L)]
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
}
cols <- round((e - w) / ewres)
rows <- round((n - s) / nsres)
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
if (inherits(y, "GSpatial")) .locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
adamlilith/fasterRaster documentation built on Oct. 26, 2024, 4:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .regionResCenterMethod .regionNumericRes .projectRaster
#' Change the coordinate reference system of a GRaster or GVector
#'
#' @description `project()` changes the coordinate reference system (CRS) of a `GRaster` or `GVector`. It has three use cases:
#' * `x` is a `GRaster` and `y` is a `GRaster`: `x` will be projected to the CRS of `y` and resampled to have the same resolution as `y`. If argument `align` is `FALSE`, then it will also be cropped to the extent of `y`.
#' * `x` is a `GRaster` and `y` is a `GVector` or a CRS string (typically in Well-Known Text format): `x` will be projected to the CRS specified by `y` and resampled but not cropped.
#' * `x` is a `GVector` and `y` is a `GRaster`, `GVector`, or CRS string: The vector will be projected to the CRS of `y`.
#'
#' @param x A `GRaster` or `GVector` to be projected.
#'
#' @param y A character or `GLocation` object (i.e., typically a `GRaster` or `GVector`): Used to set the focal `GRaster` or `GVector`'s new CRS (and resolution and possibly extent, for `GRaster`s).
#'
#' @param align Logical: If `FALSE` (default), and `x` and `y` are `GRaster`s, then the extent of `x` will be cropped to the extent of `y`. If `TRUE`, no cropping is performed.
#'
#' @param method Character or `NULL` (for `GRaster`s only): Method to use to conduct the transformation (rasters only). Partial matching is used.
#' * `NULL` (default): Automatically choose based on raster properties (`near` for categorical data, `bilinear` for continuous data).
#' * `"near"`: Nearest neighbor. Best for categorical data, and often a poor choice for continuous data. If [datatype()] is `integer`, this method will be used by default.
#' * `"bilinear"`: Bilinear interpolation (default for non-categorical data; uses weighted values from 4 cells).
#' * `"bicubic"`: Bicubic interpolation (uses weighted values from 16 cells).
#' * `"lanczos"`: Lanczos interpolation (uses weighted values from 25 cells).
#'
#' *Note #1*: If `x` and `y` are `GRaster`s, and `res = "terra"`, then the same `method` is used to resample `x` to the resolution of `y` before projecting `x`.
#'
#' *Note #2*: Methods that use multiple cells will cause the focal cell to become `NA` if there is at least one cell with an `NA` in the cells it draws from. These `NA` cells can often be filled using the `fallback` argument.
#'
#' @param fallback Logical (for projecting `GRaster`s only): If `TRUE` (default), then use "lower" resampling methods to fill in `NA` cells when a "higher" resampling method is used. For example, if `method = "bicubic"`, `NA` cells will be filled in using the `bilinear` method, except when that results in `NA`s, in which case the `near` method will be used. Fallback causes fewer cells to revert to `NA` values, so may be better at capturing complex "edges" (e.g., coastlines). Fallback does increase processing time because each "lower" method must be applied, then results merged. Fallback is not used if `method = "near"`.
#'
#' @param res Character (for projecting `GRaster`s only): Method used to set the resolution of a `GRaster` in the new CRS. This can be one of three options. Partial matching is used and case ignored:
#' * `"terra"`: This method creates an output raster that is as close as possible in values and resolution to the one that [terra::project()] would create. However, for large rasters (i.e., many cells), this can fail because `terra::project()` encounters memory limits (it is used internally to create a template). This method resamples the focal raster in its starting CRS, then projects it to the destination CRS.
#' * `"template"`: This method can only be used if `y` is a `GRaster`. The output will have the same resolution as `y` and possibly the same extent (depending on the value of `align`). However, unlike the `"terra"` method, cell values will not necessarily be as close as possible to what [terra::project()] would generate (unless `method = "near"`). Unlike the `"terra"` method, this method does not resample the focal raster in its starting CRS before projecting. For large rasters it will be faster than the `"terra"` method (especially if `"method = "near"`), and it should be less likely to fail because of memory limits.
#' * Two numeric values: Values for the new resolution (x- and y-dimensions).
#' * `"center"`: This method locates the centroid of the raster to be projected (in the same CRS as the original raster). It then creates four points north, south, east, and west of the centroid, each spaced one cell's width from the centroid. This set of points is then projected to the new CRS. The new cell size in the x-dimension will be the average of the distance between the east and west points from the centroid, and in the y-dimension the average from the centroid to the north and south points.
#' * `"fallback"` (default): This applies the `terra` method first, but if that fails, then tries `template`, then `center`. This process can take a long time for large rasters.
#'
#' @param wrap Logical:
#' * `GRaster`s: When projecting rasters that "wrap around" (i.e., whole-world rasters or rasters that have edges that actually circle around to meet on the globe), `wrap` should be `TRUE` to avoid removing rows and columns from the "edge" of the map. The default is `FALSE`.
#' * `GVector`s: When projecting vectors that span the international date line at 180E/W, `wrap` should be `TRUE` to avoid an issue where the coordinates are incorrectly mapped to the range -180 to 180.
#'
#' @param verbose Logical (for projecting `GRaster`s only): If `TRUE`, display progress. Default is `FALSE`.
#'
#' @details When projecting a raster, the "fallback" methods in **GRASS** module `r.import` are actually used, even though the `method` argument takes the strings specifying non-fallback methods. See the manual page for the `r.import` **GRASS** module.
#'
#' @returns A `GRaster` or `GVector`.
#'
#' @seealso [terra::project()], [sf::st_transform()], **GRASS** manual pages for modules `r.proj` and `v.proj` (see `grassHelp("r.proj")` and `grassHelp("v.proj")`)
#'
#' @example man/examples/ex_project.r
#'
#' @aliases project
#' @rdname project
#' @exportMethod project
methods::setMethod(
f = "project",
signature = c(x = "GRaster"),
definition = function(
x,
y,
align = FALSE,
method = NULL,
fallback = TRUE,
res = "fallback",
wrap = FALSE,
verbose = FALSE
) {
if (inherits(res, "character")) {
res <- omnibus::pmatchSafe(res, c("fallback", "terra", "template", "center"), nmax = 1L)
if (!inherits(y, "GRaster") && res == "template") stop("The `template` method for setting resolution can only be used if the argument `y` is a GRaster.")
} else if (inherits(res, "numeric")) {
if (length(res) == 1L) res <- c(res, res)
if (length(res) != 2L) stop("Argument `res` must have one or two numeric values, or be a string.")
} else {
stop("Argument `res` must be a numeric vector or a string.")
}
if (inherits(res, "numeric")) {
out <- .projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = res, wrap = wrap)
} else if (res != "fallback") {
out <- .projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = res, wrap = wrap)
} else if (res == "fallback") {
out <- FALSE
j <- 0L
n <- length(res)
try <- c("terra", "template", "center")
while (inherits(out, "logical") & j <= n) {
j <- j + 1L
thisRes <- try[j]
out <- tryCatch(
.projectRaster(x = x, y = y, align = align, method = method, fallback = fallback, res = thisRes, wrap = wrap, verbose = verbose),
error = function(cond) FALSE
)
}
}
out
} # EOF
)
#' @param x A `GRaster`.
#' @param y A `GRaster` or `GVector`.
#' @param align T/F
#' @param method Character
#' @param fallback T/F
#' @param res Character.
#' @param wrap T/F
#' @param verbose T/F
#' @noRd
.projectRaster <- function(x, y, align, method, fallback, res, wrap, verbose = FALSE) {
nLayers <- nlyr(x)
# .message(msg = "project_raster", message = "This function can produce erroneous results if the raster crosses a pole or the international date line.")
# imperfect catch for cases where `wrap` should be `TRUE` but is not
unproj <- .projection(x) %in%
c("Latitude-Longitude", "WGS84", "WGS 84", "NAD83", "NAD 83", "NAD27", "NAD 27")
if (unproj & !wrap) {
extent <- ext(x, vector = TRUE)
if ((extent[1L] == -180 & extent[2L] == 180) | (extent[3L] == -90 & extent[4L] == 90)) warning("This GRaster seems to wrap around the globe to meet at the international\n date line and/or the poles. Should `wrap` be `TRUE`?")
}
xLocation <- .location(x)
yLocation <- .locationFind(y, match = "crs", return = "name")
if (!is.null(yLocation)) {
if (yLocation == xLocation) {
warning("Object is already in the desired coordinate reference system.")
return()
}
} else if (is.null(yLocation)) {
yLocation <- .locationCreate(y)
yLocation <- .location(yLocation)
}
# method
if (!is.null(method)) method <- omnibus::pmatchSafe(method, c("nearest", "bilinear", "bicubic", "lanczos"))
if (is.null(method)) {
dt <- datatype(x)
method <- if (all(dt %in% c("integer", "factor"))) {
"nearest"
} else if (all(dt %in% c("float", "double"))) {
"bilinear"
} else {
stop("Rasters are a mix of datatype integer and non-integer (categorical and continuous). \nProject each type of raster separately.")
}
}
nSteps <- nLayers + (res == "terra")
if (verbose | faster("verbose")) {
pb <- utils::txtProgressBar(min = 0, max = nSteps, initial = 0, style = 3, width = 30)
steps <- 0
}
### If y is a GRaster, reshape region in target location using y's extent/resolution.
if (inherits(y, "GRaster")) {
# user sets target cell size
if (is.numeric(res)) {
.regionNumericRes(x = x, y = y, align = align, res = res)
# resample as per terra::project()
} else if (res == "terra") {
if (verbose | faster("verbose")) {
utils::setTxtProgressBar(pb, steps)
steps <- steps + 1
}
# Use a SpatRaster as template for resampling region. We do this so we can set the "region" resolution and extent correctly.
extent <- ext(x, vector = TRUE)
xRast <- matrix(NA_real_, nrow = nrow(x), ncol = ncol(x))
xRast <- terra::rast(xRast, crs = crs(x), extent = extent)
extent <- ext(y, vector = TRUE)
yRast <- matrix(NA_real_, nrow = nrow(y), ncol = ncol(y))
yRast <- terra::rast(yRast, crs = crs(y), extent = extent)
xRast <- terra::project(xRast, yRast, align = TRUE)
xRast <- terra::project(xRast, crs(x))
xSR <- xRast
xRast[] <- 1L
xRast <- fast(xRast)
# resample x in its native location to the resolution it will have in the target location
x <- resample(
x = x,
y = xRast,
method = method,
fallback = fallback
)
xRast <- terra::project(xSR, yRast, method = "near", align = align)
# reshape region
.locationRestore(yLocation)
.region(xRast)
# } else if (res == "dimensions") {
# .regionRespectsDims(x = x, y = y, align = align)
# same extent and resolution as terra, but without the resampling before projection
} else if (res == "template") {
if (!align) {
# reshape region
.locationRestore(y)
.region(y)
} else {
# extent from x projected to y
n <- N(x)
s <- S(x)
e <- E(x)
w <- W(x)
corners <- rbind(
c(w, n),
c(w, s),
c(e, s),
c(e, n)
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- terra::project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
ewres <- xres(y)
nsres <- yres(y)
yextent <- ext(y, vector = TRUE)
if (n > yextent[4L]) {
diff <- n - yextent[4L]
nc <- ceiling(diff / nsres)
n <- yextent[4L] + nc * nsres
}
if (s < yextent[3L]) {
diff <- yextent[3L] - s
nc <- ceiling(diff / nsres)
s <- yextent[3L] - nc * nsres
}
if (e > yextent[2L]) {
diff <- e - yextent[2L]
nc <- ceiling(diff / ewres)
e <- yextent[2L] + nc * ewres
}
if (w < yextent[1L]) {
diff <- yextent[1L] - w
nc <- ceiling(diff / ewres)
w <- yextent[1L] - nc * ewres
}
rows <- round((y - s) / nsres)
cols <- round((e - w) / ewres)
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
.locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
# "center" method of determining cell resolution
} else if (res == "center") {
.regionResCenterMethod(x = x, y = y, align = align)
}
### "y" is not a raster (i.e., WKT string or GVector)
} else {
# user-supplied resolution
if (inherits(res, "numeric")) {
.regionNumericRes(x = x, y = y, align = TRUE, res = res)
# project as per terra::project()
} else if (res == "terra") {
# make template raster to match raster to be projected
extent <- ext(x, vector = TRUE)
xRast <- matrix(NA_integer_, nrow = nrow(x), ncol = ncol(x))
xSR <- terra::rast(xRast, crs = crs(x), extent = extent)
xRast <- terra::project(xSR, crs(y), method = "near", align = align)
# reshape region
.locationRestore(yLocation)
.region(xRast)
# # respect dimensions
# } else if (res == "dimensions") {
# .regionRespectsDims(x = x, y = y, align = align)
# "center" method of determining cell resolution
} else if (res == "center") {
.regionResCenterMethod(x = x, y = y, align = align)
}
}
### project raster
##################
if (method != "nearest" & fallback) method <- paste0(method, "_f")
srcs <- .makeSourceName("project", "raster", nlyr(x))
for (i in seq_len(nlyr(x))) {
if (verbose | faster("verbose")) {
utils::setTxtProgressBar(pb, steps)
steps <- steps + 1
}
args <- list(
cmd = "r.proj",
mapset = .mapset(x),
input = sources(x)[i],
output = srcs[i],
method = method,
memory = faster("memory"),
flags = c(.quiet(), "overwrite")
)
args <- .addLocationProject(args, .location(x))
if (wrap) args$flags <- c(args$flags, "n")
do.call(rgrass::execGRASS, args=args)
if (is.factor(x)[i] & method == "nearest") {
levels <- cats(x)[[i]]
} else {
levels <- NULL
}
thisOut <- .makeGRaster(srcs[i], names(x)[i], levels = levels, ac = activeCat(x, layer = i))
if (i == 1L) {
out <- thisOut
} else {
out <- c(out, thisOut)
}
} # project next raster
.locationRestore(y)
# if using y as extent to which to crop
if (!align & inherits(y, "GRaster")) out <- crop(out, y)
if (verbose | faster("verbose")) close(pb)
out
}
#' @aliases project
#' @rdname project
#' @exportMethod project
methods::setMethod(
f = "project",
signature = c(x = "GVector"),
definition = function(
x,
y,
wrap = FALSE
) {
# imperfect catch for cases where `wrap` should be `TRUE` but is not
unproj <- .projection(x) %in%
c("Latitude-Longitude", "WGS84", "WGS 84", "NAD83", "NAD 83", "NAD27", "NAD 27")
if (unproj & !wrap) {
extent <- ext(x, vector = TRUE)
if (
(extent[1L] == -180 & extent[2L] == 180) |
(extent[3L] == -90 & extent[4L] == 90) |
(extent[1L] > 0 & extent[2L] < 0) |
(extent[3L] > extent[4L]) |
(extent[4L] < extent[3L])
) warning("This GVector seems to span the international date line and/or poles. Should `wrap` be `TRUE`?")
}
xLocation <- .location(x)
yLocation <- .locationFind(y, match = "crs", return = "name")
if (!is.null(yLocation)) {
if (yLocation == xLocation) {
warning("Object is already in the desired coordinate reference system.")
return(x)
}
} else if (is.null(yLocation)) {
yLocation <- .locationCreate(y)
yLocation <- .location(yLocation)
}
.locationRestore(yLocation)
src <- .makeSourceName("project_v_proj", "vector")
args <- list(
cmd = "v.proj",
dbase = faster("workDir"),
mapset = .mapset(x),
input = sources(x),
output = src,
flags = c(.quiet(), "overwrite")
)
args <- .addLocationProject(args, .location(x))
if (wrap) args$flags <- c(args$flags, "w") # if crosses international date line
# if (geomtype(x, grass = TRUE) == "point") args$flags <- c(args$flags, "b") # disable topology build for points
do.call(rgrass::execGRASS, args = args)
out <- .makeGVector(src, table = x@table)
out
} # EOF
)
#' Set region in target location so it has user-defined resolution
#'
#' @param x A `GRaster`.
#' @param y `GRaster` or `GVector`
#' @param align Logical
#'
#' @noRd
.regionNumericRes <- function(x, y, align, res) {
if (!align) {
# extent from y
n <- N(y)
s <- S(y)
e <- E(y)
w <- W(y)
} else {
# extent from x projected to y
n <- N(x)
s <- S(x)
e <- E(x)
w <- W(x)
corners <- rbind(
c(w, n),
c(w, s),
c(e, s),
c(e, n)
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- terra::project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
}
cols <- ceiling((e - w) / res[1L])
rows <- ceiling((n - s) / res[2L])
e <- w + cols * res[1L]
s <- n - rows * res[2L]
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
.locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
# #' Set region in target location so it has same dimensions (rows, columns) as focal raster
# #'
# #' @param x `GRaster` to project
# #' @param y `GRaster` or `GVector`.
# #' @param align Logical.
# #'
# #' @noRD
# .regionRespectsDims <- function(x, y, align) {
# # get new extent in y location
# if (!align) {
# extent <- ext(y, vector = TRUE)
# } else {
# # get new extent from projected corners of x raster
# n <- N(x)
# s <- S(x)
# e <- E(x)
# w <- W(x)
# corners <- rbind(
# c(w, n),
# c(w, s),
# c(e, s),
# c(e, n)
# )
# if (inherits(y, "GLocation")) {
# yCrs <- crs(y)
# } else {
# yCrs <- y
# }
# corners <- terra::vect(corners, type = "points", crs = crs(x))
# corners <- terra::project(corners, yCrs)
# corners <- terra::crds(corners)
# extent <- c(
# xmin = min(corners[1L:2L, 1L]),
# xmax = max(corners[2L:3L, 1L]),
# ymin = min(corners[2L:3L, 2L]),
# ymax = max(corners[c(1L, 4L), 2L])
# )
# }
# xdim <- dim(x)
# rows <- xdim[1L]
# cols <- xdim[2L]
# n <- as.character(extent[4L])
# s <- as.character(extent[3L])
# e <- as.character(extent[2L])
# w <- as.character(extent[1L])
# # reshape region
# .locationRestore(y)
# rgrass::execGRASS(
# cmd = "g.region",
# n = n, s = s, e = e, w = w,
# rows = rows, cols = cols,
# flags = c("o", .quiet())
# )
# }
#' Set target region resolution using "center" method
#'
#' Set target region resolution using "center" method. Extent will be from y or projected corners of x
#'
#' @param x `GRaster`.
#' @param x `GRaster` or `GVector`.
#' @param align Logical.
#'
#' @noRd
.regionResCenterMethod <- function(x, y, align) {
extent <- ext(x, vector = TRUE)
# get resolution in location of y
center <- c(mean(extent[1L:2L]), mean(extent[3L:4L]))
xres <- res(x)
wp <- c(center[1L] - xres[1L], center[2L])
ep <- c(center[1L] + xres[1L], center[2L])
np <- c(center[1L], center[2L] + xres[2L])
sp <- c(center[1L], center[2L] - xres[2L])
pts <- rbind(center, wp, ep, np, sp)
pts <- terra::vect(pts, type = "points", crs = crs(x))
pts <- terra::project(pts, crs(y))
pts <- terra::crds(pts)
centerx <- pts[1L, 1L]
centery <- pts[1L, 2L]
ewres <- mean(c(centerx - pts[2L, 1L], pts[3L, 1L] - centerx))
nsres <- mean(c(pts[4L, 2L] - centery, centery - pts[5L, 2L]))
# get extent in y location
if (!align & inherits(y, "GRegion")) {
# extent from y
yextent <- ext(y, vector = TRUE)
n <- N(y)
s <- S(y)
e <- E(y)
w <- W(y)
} else {
# extent by projecting corners of x to y location
corners <- rbind(
extent[c(1L, 4L)],
extent[c(1L, 3L)],
extent[c(2L, 3L)],
extent[c(2L, 4L)]
)
corners <- terra::vect(corners, type = "points", crs = crs(x))
corners <- project(corners, crs(y))
corners <- terra::crds(corners)
n <- max(corners[c(1L, 4L), 2L])
s <- min(corners[c(2L, 3L), 2L])
e <- max(corners[c(3L, 4L), 1L])
w <- min(corners[c(1L, 2L), 1L])
}
cols <- round((e - w) / ewres)
rows <- round((n - s) / nsres)
n <- as.character(n)
s <- as.character(s)
e <- as.character(e)
w <- as.character(w)
# reshape region
if (inherits(y, "GSpatial")) .locationRestore(y)
rgrass::execGRASS(
cmd = "g.region",
n = n, s = s, e = e, w = w,
rows = rows, cols = cols,
flags = c("o", .quiet())
)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.