#' Binarise the values in a raster
#'
#' Transform a raster to binary form.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param thresh [\code{integerish(1)}]\cr Value above which the cell will be set
#' to 1, below which it will be set to 0.
#' @param match [\code{integerish(.)}]\cr One or more values which will be set to
#' 1, while the remaining values will be set to 0.
#' @return A binary \code{RasterLayer} of the same dimension as \code{obj}.
#' @details A binary object contains only values 1 (foreground) and 0
#' (background).
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#' visualise(rBinarise(input, thresh = 30))
#'
#' visualise(rBinarise(input, match = 1), new = TRUE)
#' @importFrom checkmate assertClass assertNumber assertNumeric
#' @importFrom raster as.matrix raster extent extent<- crs crs<-
#' @export
rBinarise <- function(obj, thresh = NULL, match = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumber(thresh, null.ok = TRUE)
assertNumeric(match, null.ok = TRUE)
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
if(is.null(thresh)){
thresh <- max(vals)
} else{
if(!min(vals, na.rm = TRUE)<=thresh | !thresh<=max(vals, na.rm = TRUE)){
stop("please provide a value for 'thresh' within the range of the values of 'obj'.")
}
}
# in case 'match' is defined, remove these values from 'values'
if(!is.null(match)){
values <- vals[!vals %in% match]
} else{
values <- vals
}
# select only values that are above 'thresh'
temp <- morphC(mat = mat, kernel = matrix(thresh, 1, 1), value = values, blend = 4,
merge = 3, rotateKernel = FALSE, strictKernel = FALSE)
# in case there are several values in 'match', the matrix needs to be flattened
if(length(match) > 1){
temp <- morphC(mat = temp, kernel = matrix(0, 1, 1), value = vals, blend = 4,
merge = 3, rotateKernel = FALSE, strictKernel = FALSE)
}
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history) == 0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("values have been binarised")))
names(out) <- "values_binarised"
return(out)
}
#' Blend two rasters with each other
#'
#' Blend two rasters simply by adding (and weighing) their values per cell.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param overlay [\code{RasterLayer(1)} | \code{matrix(1)}]\cr the raster that
#' should be blended with the primary raster (obj); has to have the same
#' dimension as the primary raster.
#' @param patches [\code{logical(1)}]\cr should the blend be carried out while
#' considering foreground patches in the primary raster (\code{TRUE}), or
#' should the blend be carried out without this consideration (\code{FALSE},
#' default).
#' @param weight [\code{numeric(1)}]\cr the weight the overlay raster should
#' have as multiple (or fraction) of the primary raster.
#' @return A \code{RasterLayer} of the same dimension as \code{primary}, in
#' which an overlay has been blended with the primary raster.
#' @details In case \code{patches = TRUE}, the cells in the overlay are grouped
#' according to the foreground patches in the primary raster and all cells of
#' this group are assigned their common average, then the blend is carried
#' out.
#'
#' A blend operation currently is defined as 'overlay*weight + obj'. To
#' aggregate more than two \code{RasterLayers} with arbitrary functions, see
#' \code{\link{rReduce}}.
#' @family operators to modify a raster
#' @examples
#' # define primary raster ...
#' input <- rtRasters$continuous
#' patches <- rPatches(rBinarise(input, thresh = 30))
#'
#' # ... and an overlay from a matrix
#' m <- matrix(nrow = 56, ncol = 60, data = 0)
#' m[c(5:25), c(5:50)] <- 10
#' mask <- raster::raster(m, xmn=0, xmx=60, ymn=0, ymx=56, crs=NA)
#'
#' # blend while considering patches in the primary raster
#' blended <- rBlend(patches, overlay = mask, patches = TRUE)
#' visualise(raster::stack(patches, blended))
#'
#' # blend while not considering patches
#' visualise(rBlend(patches, overlay = mask, weight = 10))
#' @importFrom checkmate assertClass assert testClass assertLogical assertNumber
#' assertTRUE
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rBlend <- function(obj, overlay, patches = FALSE, weight = 1){
# check arguments
assertClass(obj, "RasterLayer")
assert(
isRaster <- testClass(overlay, "RasterLayer"),
isMatrix <- testClass(overlay, "matrix")
)
assertLogical(patches)
assertNumber(weight)
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
if(isRaster){
mat_overlay <- as.matrix(overlay)
} else{
mat_overlay <- overlay
}
assertTRUE(all(dim(mat)==dim(mat_overlay)))
if(patches){
clusterIds <- vals[!is.na(vals)]
for(i in clusterIds){
mat_overlay[mat%in%i] <- mean(overlay[mat%in%i], na.rm = TRUE)
}
}
# I keep this here because at some point I might implement '#include <boost/...>' in reduceMatrixC
# matList <- list(mat, mat_overlay*weight)
# temp <- reduceMatrixC(matList, f = sum)
temp <- mat_overlay*weight + mat
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the raster has been blended")))
names(out) <- "blended"
return(out)
}
#' Select cells with values between an upper and lower threshold in a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param range [\code{numeric(2)}]\cr minimum and maximum value between which
#' the values will be selected.
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, where the
#' cells with values covered by \code{range} retain their value and other
#' cells are set to \code{NA}.
#' @family operators to select a subset of cells
#' @examples
#' input <- rtRasters$continuous
#' visualise(rBounded(input, range = c(40, 80)))
#' @importFrom checkmate assertClass assertNumeric
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rBounded <- function(obj, range = NULL, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumeric(range, len = 2)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
lower <- min(range, na.rm = TRUE)
upper <- max(range, na.rm = TRUE)
if(!min(vals)<lower | !lower<max(vals) | !min(vals)<upper | !upper<max(vals)){
stop("please provide values for 'lower' and 'upper' within the range of the values of 'obj'.")
}
values <- vals[!duplicated(vals)]
temp <- morphC(mat = mat, kernel = matrix(lower, 1, 1), value = values, blend = 4,
merge = 12, rotateKernel = FALSE, strictKernel = FALSE)
temp <- morphC(mat = temp, kernel = matrix(upper, 1, 1), value = values, blend = 3,
merge = 12, rotateKernel = FALSE, strictKernel = FALSE)
# set na
temp[is.na(temp)] <- background
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("values greater than ", lower, " and less than ", upper, " have been selected")))
names(out) <- paste0("between_", lower, "_and_", upper)
return(out)
}
#' Assign categories to the values in a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param breaks [\code{integerish(.)}]\cr the values, where categories should
#' be delimited.
#' @param n [\code{integerish(1)}]\cr number of categories.
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, in which
#' the cells have the category number into which their values fall.
#' @details Using \code{n} will determine \code{breaks} based on the value-range
#' of \code{obj} so that the values are assigned to n categories.
#'
#' Assigning \code{breaks} is mostly usefull when values are to be non-linear,
#' such as \code{log(min:max)*max/log(max)}, but could also be \code{seq(min,
#' max, length.out = 21)}, which corresponds to \code{n = 20}.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#' visualise(rCategorise(input, n = 5))
#'
#' # use as algorithm in modify, to combine two iterations in one run
#' algorithm <- list(list(operator = "rCategorise", breaks = c(25, 50, 75, 90)),
#' list(operator = "rCategorise", breaks = log(1:5)*5/log(5)*20))
#'
#' obj_mod <- modify(input, by = algorithm, merge = TRUE)
#' visualise(obj_mod)
#' @importFrom checkmate assertClass assertNumeric assertNumber
#' @importFrom raster as.matrix values<- values crs crs<-
#' @export
rCategorise <- function(obj, breaks = NULL, n = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumeric(breaks, null.ok = TRUE)
assertIntegerish(n, null.ok = TRUE)
out <- obj
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
theRange <- range(vals, na.rm = T)
if(is.null(breaks)){
breaks <- seq(theRange[1], theRange[2], length.out = n+1)
}
# manage the tails
if(!any(breaks==theRange[1])){
breaks <- c(theRange[1], breaks)
}
if(!any(breaks==theRange[2])){
breaks <- c(breaks, theRange[2])
}
values(out) <- findInterval(values(out), breaks, rightmost.closed = TRUE)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("categories based on the breaks ", paste0(breaks, collapse = "/"), " have been defined")))
names(out) <- paste0("new_breaks")
return(out)
}
#' Determine the centroid of patches in a raster
#'
#' The centroid is the average location of all cells of a foreground patch.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param output [\code{character(1)}]\cr the type in which the output should be
#' returned. Either \code{"raster"} (default) or \code{"geom"}.
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return Depending on \code{output}, either \enumerate{ \item a
#' \code{RasterLayer} of the same dimension as \code{obj}, in which the
#' centroid of each foreground patch has the value of the patch-number
#' identified with \code{\link{rPatches}} and where the remaining cells of
#' each foreground patch have the value NA or \item a geometry object
#' of the centroids.}
#' @family operators to determine objects
#' @examples
#' input <- rtRasters$continuous
#' patches <- rPatches(rBinarise(input, thresh = 30))
#' visualise(rCentroid(patches))
#'
#' # get centroid coordinates
#' centroid_coords <- rCentroid(patches, output = "geom")
#' visualise(raster = patches, geom = centroid_coords, linecol = "red")
#' @importFrom checkmate assertClass
#' @importFrom raster values rasterToPoints rasterize crs
#' @importFrom methods new
#' @importFrom tibble tibble
#' @export
rCentroid <- function(obj, output = "raster", background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
output <- match.arg(output, c("raster", "geom"))
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
if(isBinaryC(mat)){
mat <- components(mat, shapeKernel(c(3, 3), type = "diamond"))
}
containsNA <- ifelse(anyNA(mat), TRUE, FALSE)
if(containsNA){
mat[is.na(mat)] <- 0
}
vals <- unique(as.vector(mat))
# determine centroids by averaging all cell coordinates per patch
thePoints <- cellToPointsC(mat)
theMeans <- NULL
for(i in vals){
temp <- thePoints[thePoints$value == i,]
theMeans <- rbind(theMeans, round(colMeans(temp)/0.5)*0.5)
}
theMeans <- as_tibble(theMeans)
theMeans <- theMeans[order(theMeans$value),]
if(containsNA){
theMeans <- theMeans[-which(theMeans$value == 0),]
}
if(output == "geom"){
theMeans <- tibble(fid = theMeans$value, vid = rep(1, length(theMeans$value)), x = theMeans$x, y = theMeans$y)
out <- new(Class = "geom",
type = "point",
coords = theMeans,
attr = tibble(fid = unique(theMeans$fid), n = 1),
window = tibble(x = c(0, dim(obj)[2]), y = c(0, dim(obj)[1])),
scale = "absolute",
crs = as.character(NA),
history = list(paste0("the centroids of patches have been determined")))
return(out)
} else{
out <- rasterize(x = theMeans[c(1, 2)], y = obj, background = background)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the centroids of patches have been determined")))
names(out) <- "centroids"
return(out)
}
}
#' Morphologically dilate foreground patches in a raster
#'
#' The morphological operation 'dilate' increases the value of cells along a
#' gradient of values.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param kernel [\code{matrix(1)}]\cr scan the raster with this kernel (default
#' is a 3 by 3 cells diamond kernel).
#' @details A morphological dilation changes the values in a raster after
#' comparing each value to a kernel that moves systematically across the
#' raster. A kernel is any 2D array with an odd number of cells in each
#' dimension and a focal cell in the middle. The cells of a kernel can have
#' four categories of values: 0, 1, NA and any value > 1. When the kernel
#' values match the values of the input raster, the input raster cell covered
#' by the focal kernel cell is modified according to a specific operation, for
#' instance the smallest value covered by the kernel is assigned.
#'
#' Rasters are either binary or non-binary (greyscale) rasters. \itemize{
#' \item For binary rasters the kernel cells with values 0 and 1 are matched
#' accurately in the input raster. Cells with value NA will be ignored, i.e.
#' "it does not matter what value is there". \item For a non-binary raster the
#' value 0 does not have any meaning and every 0 is turned into an NA. In case
#' a non-binary kernel is defined to dilate a non-binary (greyscale) raster,
#' each value of the input raster that is currently covered by the kernel is
#' summarised with the focal kernel value and the maximum of this set of
#' values is assigned to the input raster cell covered by the focal kernel
#' cell.}
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, where a
#' dilation has been performed.
#' @family operators to morphologically modify a raster
#' @examples
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#' visualise(rDilate(binarised))
#'
#' # use another kernel
#' (myKernel <- matrix(1, 3, 3))
#' visualise(rDilate(binarised, kernel = myKernel), new = TRUE)
#'
#' # dilate also non-binary rasters
#' visualise(raster::stack(input, rDilate(obj = input)))
#' @importFrom checkmate assertClass assertMatrix
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rDilate <- function(obj, kernel = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertMatrix(kernel, mode = "integerish", all.missing = FALSE, null.ok = TRUE)
if(is.null(kernel)) kernel <- matrix(c(0, 1, 0, 1, 1, 1, 0, 1, 0), 3, 3)
mat <- as.matrix(obj)
blend <- 1 # morphC::blendIdentity
if(!isBinaryC(mat)){
vals <- values(obj)
values <- vals[!duplicated(vals) & vals!=0]
if(!isBinaryC(kernel)){
blend <- 5 # morphC::blendPlus
}
} else{
values <- c(0)
}
kernel[kernel==0] <- NA
temp <- morphC(mat = mat, kernel = kernel,
value = values, blend = blend,
merge = 2, # morphC::mergeMax
rotateKernel = FALSE, strictKernel = FALSE)
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the raster has been morphologically dilated")))
names(out) <- paste0("dilated")
return(out)
}
#' Calculate the distance map for a raster
#'
#' The distance map of a binarised raster contains the distance of each
#' background cell to the nearest foreground cell.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param method [\code{character(1)}]\cr the distance measure to be calculated.
#' Either \code{"euclidean"} (default), \code{"manhatten"} or
#' \code{"chessboard"} distance.
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, where the
#' value of the background cells has been replaced with the distance to the
#' nearest foreground cell.
#' @details In contrast to \code{\link[raster]{distance}}, the distance
#' values here do not warp around the boundaries of the map.
#'
#' @references Meijster, A., Roerdink, J.B.T.M., Hesselink, W.H., 2000. A
#' general algorithm for computing distance transforms in linear time, in:
#' Goutsias, J., Vincent, L., Bloomberg, D.S. (Eds.), Mathematical Morphology
#' and Its Applications to Image and Signal Processing. Springer, pp. 331–340.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#'
#' # the different distance metrics
#' binarised <- rBinarise(input, thresh = 40)
#' disEuc <- rDistance(binarised)
#' disMan <- rDistance(binarised, method = "manhattan")
#' disChb <- rDistance(binarised, method = "chessboard")
#'
#' distances <- raster::stack(binarised, disEuc, disMan, disChb)
#' visualise(distances)
#'
#' # calculate distance from edge to patch interior
#' inverted <- rPermute(binarised)
#' visualise(rDistance(inverted))
#' @importFrom checkmate assertClass assertIntegerish
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rDistance <- function(obj, method = "euclidean"){
# check arguments
assertClass(obj, "RasterLayer")
method <- match.arg(method, c("euclidean", "manhattan", "chessboard"))
mat <- as.matrix(obj)
if(!isBinaryC(mat)){
stop("spatial object is not binary, please provide a binarised RasterLayer.")
}
temp <- meijsterDistanceC(mat, method = method)
if(method=="euclidean"){
temp <- sqrt(temp)
}
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("distance values have been calculated according to the ", method, " distance")))
names(out) <- paste0(method, "_distance")
# manage the bibliography entry
bib <- bibentry(bibtype = "incollection",
title = "A general algorithm for computing distance transforms in linear time",
volume = "18",
isbn = "978-0-306-47025-7",
booktitle = "Mathematical Morphology and its Applications to Image and Signal Processing",
publisher = "Springer",
author = c(
person("A", "Meijster"),
person(c("J", "B", "T", "M"), "Roerdink"),
person(c("W", "H"), "Hesselink")
),
editor = c(
person("John", "Goutsias"),
person("Luc", "Vincent"),
person(c("Dan", "S"), "Bloomberg")
),
year = "2000",
pages = "331--340"
)
if(is.null(getOption("bibliography"))){
options(bibliography = bib)
} else{
currentBib <- getOption("bibliography")
if(!bib%in%currentBib){
options(bibliography = c(currentBib, bib))
}
}
return(out)
}
#' Morphologically erode foreground patches in a raster
#'
#' The morphological operation 'erode' removes cells at the boundary of a
#' foreground patch in a binarised raster.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param kernel [\code{matrix(1)}]\cr scan the raster with this kernel (default
#' is a 3 by 3 cells diamond kernel).
#' @details A morphological dilation changes the values in a raster after
#' comparing each value to a kernel that moves systematically across the
#' raster. A kernel is any 2D array with an odd number of cells in each
#' dimension and a focal cell in the middle. The cells of a kernel can have
#' four categories of values: 0, 1, NA and any value > 1. When the kernel
#' values match the values of the input raster, the input raster cell covered
#' by the focal kernel cell is modified according to a specific operation, for
#' instance the smallest value covered by the kernel is assigned.
#'
#' Rasters are either binary or non-binary (greyscale) rasters. \itemize{
#' \item For binary rasters the kernel cells with values 0 and 1 are matched
#' accurately in the input raster. Cells with value NA will be ignored, i.e.
#' "it does not matter what value is there". \item For a non-binary raster the
#' value 0 does not have any meaning and every 0 is turned into an NA. In case
#' a non-binary kernel is defined to erode a non-binary (greyscale) raster,
#' from each value of the input raster that is currently covered by the kernel
#' is the focal kernel value is subtracted and the minimum of this set of
#' values is assigned to the input raster cell covered by the focal kernel
#' cell.}
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, where an
#' erosion has been performed.
#' @family operators to morphologically modify a raster
#' @examples
#' # use as standalone
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#' visualise(rErode(obj = binarised))
#'
#' # use another kernel
#' (myKernel <- matrix(1, 3, 3))
#' visualise(rErode(binarised, kernel = myKernel), new = TRUE)
#'
#' # dilate also non-binary rasters
#' visualise(raster::stack(input, rErode(obj = input)))
#' @importFrom checkmate assertClass assertMatrix
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rErode <- function(obj, kernel = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertMatrix(kernel, mode = "integerish", all.missing = FALSE, null.ok = TRUE)
if(is.null(kernel)) kernel <- matrix(c(0, 1, 0, 1, 1, 1, 0, 1, 0), 3, 3)
mat <- as.matrix(obj)
blend <- 1 # morphC::blendIdentity
if(!isBinaryC(mat)){
vals <- values(obj)
values <- vals[!duplicated(vals) & vals!=0]
if(!isBinaryC(kernel)){
blend <- 6 # morphC::blendMinus
}
} else{
values <- c(1)
}
kernel[kernel==0] <- NA
temp <- morphC(mat = mat, kernel = kernel,
value = values, blend = blend,
merge = 1, # morphC::mergeMin
rotateKernel = FALSE, strictKernel = FALSE)
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the raster has been morphologically eroded")))
names(out) <- paste0("eroded")
return(out)
}
#' Fill NA values in a raster
#'
#' Assign a value to cells with NA
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param with [\code{integerish(1)}]\cr value with which NA should be replaced.
#' @return a \code{RasterLayer} of the same dimension as \code{obj}, where value
#' \code{NA} has been replaced by value \code{with}.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#' patches <- rPatches(rBinarise(input, thresh = 30))
#'
#' visualise(raster::stack(patches, rFillNA(patches)))
#' @importFrom checkmate assertClass assertIntegerish
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rFillNA <- function(obj, with = 0){
# check arguments
assertClass(obj, "RasterLayer")
assertIntegerish(with)
mat <- as.matrix(obj)
mat[is.na(mat)] <- with
out <- raster(mat)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("NA has been replaced with ", with)))
names(out) <- paste0("NA_replaced")
return(out)
}
#' Select cells with values above a threshold in a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param thresh [\code{numeric(1)}]\cr values above this threshold should be
#' retained.
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' all values lower than \code{thresh} have been set to \code{background}.
#' @family operators to select a subset of cells
#' @examples
#' input <- rtRasters$continuous
#' visualise(rGreater(input, thresh = 40))
#' @importFrom checkmate assertClass assertIntegerish
#' @importFrom raster as.matrix values raster extent crs crs<-
#' @importFrom methods as
#' @export
rGreater <- function(obj, thresh, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumeric(thresh)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
if(!min(vals, na.rm = TRUE)<=thresh | !thresh<=max(vals, na.rm = TRUE)){
stop("please provide a value for 'thresh' within the range of the values of 'obj'.")
}
temp <- morphC(mat = mat, kernel = matrix(thresh, 1, 1), value = vals, blend = 4,
merge = 12, rotateKernel = FALSE, strictKernel = FALSE)
# set na
temp[is.na(temp)] <- background
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("values greater than ", thresh, " have been selected")))
names(out) <- paste0("greater_than_", thresh)
return(out)
}
#' Select cells with values below a threshold in a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param thresh [\code{numeric(1)}]\cr values below this threshold should be
#' retained.
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' all values greater than \code{thresh} have been set to \code{NA}.
#' @family operators to select a subset of cells
#' @examples
#' input <- rtRasters$continuous
#' visualise(rLess(input, thresh = 80))
#' @importFrom checkmate assertClass assertIntegerish
#' @importFrom raster as.matrix values raster extent crs crs<-
#' @export
rLess <- function(obj, thresh, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumeric(thresh)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
if(!min(vals, na.rm = TRUE)<=thresh | !thresh<=max(vals, na.rm = TRUE)){
stop("please provide a value for 'thresh' within the range of the values of 'obj'.")
}
temp <- morphC(mat = mat, kernel = matrix(thresh, 1, 1), value = vals, blend = 3,
merge = 12, rotateKernel = FALSE, strictKernel = FALSE)
# set na
temp[is.na(temp)] <- background
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("values less than ", thresh, " have been selected")))
names(out) <- paste0("less_than_", thresh)
return(out)
}
#' Select cells of a raster based on a mask
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param mask [\code{RasterLayer(1)} | \code{matrix(1)} | \code{geom(1)}]\cr
#' Either binary object of the same dimension as \code{obj} where the cells
#' that should be retained have the value 1 and all other cells the value 0 or\cr
#' \code{geom} where values inside the geom should be retained.
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' all cells with value 0 in the mask have been set to \code{NA} and all other
#' values are retained.
#' @details If used in an algorithm, \code{mask} can also contain the name of a
#' sub-algorithm to use the final output thereof as mask. Moreover, \code{mask
#' = "input"} would select the original raster as mask.
#' @family operators to select a subset of cells
#' @examples
#' input <- rtRasters$continuous
#' m <- matrix(nrow = 56, ncol = 60, data = 0)
#' m[c(5:25), c(5:50)] <- 1
#'
#' visualise(rMask(input, mask = m))
#'
#' \dontrun{
#'
#' # determine mask interactively
#' mask <- geomPolygon(template = input, vertices = 5, show = T, col = "deeppink")
#' mask <- gToRaster(mask)
#'
#' visualise(raster = rMask(obj = input, mask = mask))
#' }
#' @importFrom checkmate assertClass
#' @importFrom geometr gc_sp
#' @importFrom raster as.matrix raster extent crs crs<- xres yres colortable<- mask
#' @export
rMask <- function(obj, mask = NULL, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
isRaster <- testClass(mask, "RasterLayer")
isMatrix <- testClass(mask, "matrix")
isGeom <- testClass(mask, "geom")
isSp <- testClass(mask, "Spatial")
isSf <- testClass(mask, "sf")
if(!isRaster & !isMatrix & !isGeom & !isSp & !isSf){
stop("please provide a mask of the supported type.")
}
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
temp <- as.matrix(obj)
ext <- extent(obj)
target_crs <- crs(obj)
rat <- levels(obj)[[1]]
cltbl <- colortable(obj)
if(isGeom){
mask <- gc_sp(geom = mask)
isSp <- TRUE
isGeom <- FALSE
}
if(isSp | isSf){
out <- mask(x = obj, mask = mask)
} else{
if(isRaster){
mask <- as.matrix(mask)
}
if(!isBinaryC(mask)){
stop("please provide a binary mask.")
}
if(!all(dim(temp)==dim(mask))){
stop("please provide a mask of the same dimension as 'obj'.")
}
# modify obj
temp[mask==0] <- background
out <- raster(temp,
xmn = ext[1], xmx = ext[2], ymn = ext[3], ymx = ext[4],
crs = target_crs)
}
# assign attribute table
if(!is.null(rat)){
tempRat <- rat[rat$id %in% unique(values(out)),]
out@data@isfactor <- TRUE
out@data@attributes <- list(tempRat)
}
# assign colourtable
if(!length(cltbl) == 0){
colortable(out) <- cltbl
}
# assign history
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the raster has been masked")))
names(out) <- "masked"
return(out)
}
#' Match cells of a raster with a kernel
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param kernel [\code{matrix(.)} | \code{list(.)} thereof]\cr scan the
#' raster with this kernel (default is a 3 by 3 cells diamond kernel).
#' @param rotate [\code{logical(1)}]\cr should the kernel be applied for all
#' possible rotations (\code{TRUE}, default) or should the kernel be used as
#' is (\code{FALSE})?
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} or \code{RasterStack} of the same dimension as
#' \code{obj} in which all cells that match with the kernel(s) have the kernel
#' value and all other cells have the value \code{background}.
#' @details This is also known as the 'hit-or-miss'-transform.\cr\cr Wrapper of
#' \code{\link{rMorph}} with blend = 2 and merge = 12.
#'
#' The cells of a kernel can have the values 0, 1 and NA. The kernel
#' cells with values 0 and 1 are matched accurately in the input raster. The
#' kernel cells with value NA will be ignored.
#' @family operators to select a subset of cells
#' @examples
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#'
#' # match isolated cells
#' (kIso <- matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0), 3, 3))
#' matched <- rMatch(binarised, kernel = kIso)
#' visualise(matched)
#'
#' # match all right angled corners
#' (kCorner <- matrix(c(NA, 1, NA, 0, 1, 1, 0, 0, NA), 3, 3))
#' matched <- rMatch(binarised, kernel = kCorner, background = 0)
#' visualise(matched, new = TRUE)
#'
#' # match north-east facing right angled corners and isolated cells
#' matched <- rMatch(binarised, kernel = list(kIso, kCorner), rotate = FALSE)
#' visualise(matched, new = TRUE)
#'
#' # match endpoints of a skeleton
#' skeletonised <- rSkeletonise(binarised, background = 0)
#' (kEnd <- matrix(c(NA, 0, 0, NA, 1, 0, NA, 0, 0), 3, 3))
#' endpoints <- rMatch(skeletonised, kernel = kEnd, background = 0)
#'
#' # match triple points (conjunctions) of a skeleton
#' kConj <- list(matrix(c(NA, 0, 1, 1, 1, NA, NA, 0, 1), 3, 3),
#' matrix(c(1, NA, 1, NA, 1, NA, NA, NA, 1), 3, 3),
#' matrix(c(NA, 1, NA, 0, 1, 1, 1, 0, NA), 3, 3))
#' conjunctions <- rMatch(skeletonised, kernel = kConj, background = 0)
#' out <- raster::stack(
#' rBlend(skeletonised, overlay = endpoints),
#' rBlend(skeletonised, overlay = conjunctions[[1]]),
#' rBlend(skeletonised, overlay = conjunctions[[2]]),
#' rBlend(skeletonised, overlay = conjunctions[[3]]))
#' names(out) <- c("endpoints", "conj1", "conj2", "conj3")
#' visualise(out)
#' @export
rMatch <- function(obj, kernel = NULL, rotate = TRUE, background = NULL){
out <- rMorph(obj, kernel = kernel, blend = "equal", merge = "sumNa",
rotate = rotate, background = background)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("cells have been matched with a ", dim(kernel)[1], "x", dim(kernel)[2], " kernel with values ", paste0(as.vector(kernel), collapse = " "))))
return(out)
}
#' Morphologically modify a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param kernel [\code{matrix(.)} | \code{list(.)} thereof]\cr scan the raster
#' with this kernel (default is a 3 by 3 cells diamond kernel).
#' @param blend [\code{character(1)}]\cr \code{identity}, \code{equal},
#' \code{lower}, \code{greater}, \code{plus}, \code{minus}, \code{product};
#' see Details.
#' @param merge [\code{character(1)}]\cr \code{min}, \code{max}, \code{all},
#' \code{any}, \code{sum}, \code{mean}, \code{median}, \code{sd}, \code{cv},
#' \code{one}, \code{zero}, \code{na}; see Details.
#' @param strictKernel [\code{logical(1)}]\cr
#' @param rotate [\code{logical(1)}]\cr should the kernel be applied for all
#' possible rotations (\code{TRUE}, default) or should the kernel be used as
#' is (\code{FALSE})?
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @details The morphC function (internal) is the basis of many modify
#' operations in \code{rasterTools} and \code{rMorph} exposes a fully
#' functional interface of this C++ function. The morphC function iteratively
#' goes through each pixel of a raster (\code{obj}) and compares the kernel
#' with the raster at that location. The result of this comparison depends on
#' the arguments \code{blend} and \code{merge}: \itemize{ \item First, all
#' values in raster that are covered by the kernel (which may be larger than
#' one pixel) are "cut out" and summarised pairwise with \code{kernel} (i.e.
#' they are \emph{blended}) by the function defined in \code{blend} (e.g.
#' "+"), leaving as many values as cells in \code{kernel}. \item Then these
#' values are summarised into a single value (i.e. they are \emph{merged}) by
#' the function defined in \code{merge} (e.g. "mean") and the resulting value
#' is assigned in the current location in the raster.}
#'
#' The following functions are defined for \code{blend}: \enumerate{ \item
#' identity: the value of \code{obj} where \code{kernel} is not NA. \item
#' equal: the value 1 where \code{obj} and \code{kernel} are equal, otherwise
#' the value 0. \item lower: the values of \code{obj} that are lower than
#' \code{kernel}, otherwise 0. \item greater: the values of \code{obj} that
#' are greater than \code{kernel}, otherwise 0. \item plus: the values of
#' \code{obj} added to the values of \code{kernel}. \item minus: the values of
#' \code{kernel} subtracted from the values of \code{obj}. \item product: the
#' product of the values of \code{obj} and \code{kernel}.} The following
#' functions are defined for \code{merge}: \enumerate{ \item min: the minimum
#' value. \item max: the maximum value. \item all: the value 1 if all non-NA
#' values are not 0, otherwise 0. \item any: the value 1 if any of the non-NA
#' values are not 0, otherwise 0. \item sum: the sum of all non-NA values.
#' \item mean: the mean of all non-NA values. \item median: the median of all
#' non-NA values. \item sd: the standard deviation of all non-NA values. \item
#' cv: the coefficient of variation of all non-NA values. \item sumNa: if the
#' sum of all values is greater than 0 than this sum, otherwise NA.}
#' @references Credit for the original idea/architecture of the C++ part of this
#' function goes to Jon Clayden
#' (\href{https://github.com/jonclayden/mmand}{R::mmand}). The functionality
#' has been slightly extended here.
#' @family operators to morphologically modify a raster
#' @importFrom raster as.matrix
#' @importFrom checkmate assertClass testClass testList assertList assertMatrix
#' assertIntegerish assertLogical
#' @export
rMorph <- function(obj, kernel = NULL, blend = NULL, merge = NULL, rotate = TRUE,
strictKernel = TRUE, background = NULL){
assertClass(obj, "RasterLayer")
if(is.null(kernel)){
kernel <- matrix(c(0, 1, 0, 1, 1, 1, 0, 1, 0), 3, 3)
isList <- FALSE
} else{
isMatrix <- testClass(kernel, "matrix")
isList <- testList(kernel)
}
if(isList){
assertList(kernel, types = "matrix")
} else{
assertMatrix(kernel, null.ok = TRUE)
}
assertSubset(blend, choices = c("identity", "equal", "lower", "greater", "plus", "minus", "product"))
blendID <- which(c("identity", "equal", "lower", "greater", "plus", "minus", "product") %in% blend)
assertSubset(merge, choices = c("min", "max", "all", "any", "!all", "!any", "sum", "mean", "median", "sd", "cv", "sumNa"))
mergeID <- which(c("min", "max", "all", "any", "!all", "!any", "sum", "mean", "median", "sd", "cv", "sumNa") %in% merge)
assertLogical(rotate, any.missing = FALSE)
assertLogical(strictKernel, any.missing = FALSE)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
if(isList){
temp <- lapply(seq_along(kernel), function(i){
temp2 <- morphC(mat = mat, kernel = kernel[[i]], value = vals, blend = blendID,
merge = mergeID, rotateKernel = rotate, strictKernel = TRUE)
temp2[is.na(temp2)] <- background
temp2 <- raster(temp2)
names(temp2) <- paste0("morphed")
return(temp2)
})
out <- stack(temp)
} else{
temp <- morphC(mat = mat, kernel = kernel, value = vals, blend = blendID,
merge = mergeID, rotateKernel = rotate, strictKernel = TRUE)
temp[is.na(temp)] <- background
out <- raster(temp)
names(out) <- paste0("morphed")
}
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, paste0("the object has been morphologically modified (blend:", blend, ",merge:", merge, ")"))
return(out)
}
#' Offset the values in a raster.
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param value [\code{integerish(1)}]\cr the value by which to offset all
#' values in the raster (can also be negative).
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, in which
#' all values have been offsetted by \code{value}.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#' visualise(rOffset(binarised, value = 2))
#' @importFrom checkmate assertClass assertIntegerish
#' @importFrom raster as.matrix extent crs crs<-
#' @export
rOffset <- function(obj, value = 1){
assertClass(obj, "RasterLayer")
assertIntegerish(value, any.missing = FALSE, len = 1)
mat <- as.matrix(obj)
temp <- mat + value
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, paste0("the raster values have been offset by ", value, "."))
names(out) <- paste0("offset_by_", value)
return(out)
}
#' Determine foreground patches in a raster
#'
#' (Foreground) Patches are sets of cells which are connected either directly or
#' via other cells in a binarised raster and which should hence be treated as
#' distinct objects.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param kernel [\code{matrix(1)}]\cr scan the raster with this kernel (default
#' is a 3 by 3 cells diamond kernel).
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, in which
#' neighbouring cells of the foreground have been assigned the same value,
#' forming patches.
#' @family operators to determine objects
#' @examples
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#' visualise(rPatches(binarised))
#' @importFrom checkmate assertClass assertMatrix assertIntegerish
#' @importFrom mmand components
#' @importFrom raster as.matrix raster extent crs crs<-
#' @export
rPatches <- function(obj, kernel = NULL, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertMatrix(kernel, mode = "integerish", all.missing = FALSE, null.ok = TRUE)
if(is.null(kernel)) kernel <- matrix(c(0, 1, 0, 1, 1, 1, 0, 1, 0), 3, 3)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
mat <- as.matrix(obj)
if(!isBinaryC(mat)){
stop("spatial object is not binary, run 'rBinarise()' first.")
}
temp <- components(mat, kernel)
# set na
temp[is.na(temp)] <- background
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("patches have been determined")))
names(out) <- "patches"
# manage the bibliography entry (mmand and/or the algorithm)
bib <- bibentry(bibtype = "Manual",
title = "mmand: Mathematical Morphology in Any Number of Dimensions",
author = person("Jon", "Clayden"),
year = "2017",
note = "R package version 1.5.3",
url = "https://CRAN.R-project.org/package=mmand"
)
if(is.null(getOption("bibliography"))){
options(bibliography = bib)
} else{
currentBib <- getOption("bibliography")
if(!bib%in%currentBib){
options(bibliography = c(currentBib, bib))
}
}
return(out)
}
#' Apply a permutation to the cell values of a raster
#'
#' The permutation of a set of cell values leads to a systematic (re)arrangement
#' of all the members of the set.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param type [\code{character(1)}]\cr the type of permutation that should be
#' applied. Either \code{"invert"}, \code{"revert"}, \code{"descending"},
#' \code{"ascending"} or \code{"cycle"}.
#' @param by [\code{integerish(1)}]\cr value by which to apply the permutation;
#' recently only for \code{type = "cycle"}.
#' @return A \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' all values have been permuted according to the chosen \code{type}.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#'
#' # invert background/foreground
#' visualise(raster::stack(input, rPermute(input)))
#'
#' # sort patch values ascending
#' patches <- rPatches(rBinarise(input, thresh = 30))
#' visualise(raster::stack(patches, rPermute(patches, type = "ascending")))
#'
#' # cycle values backwards by 30
#' visualise(raster::stack(input, rPermute(input, type = "cycle", by = -30)))
#' @importFrom checkmate assertClass
#' @importFrom raster values as.matrix crs crs<-
#' @export
rPermute <- function(obj, type = "invert", by = NULL){
# check arguments
assertClass(obj, "RasterLayer")
type <- match.arg(type, c("invert", "revert", "descending", "ascending", "cycle"))
if(type == "cycle"){
assertIntegerish(by)
}
mat <- as.matrix(obj)
vals <- unique(as.vector(mat))
vals <- vals[!is.na(vals)]
if(type == "invert"){
newVals <- max(vals, na.rm = TRUE) - vals
action <- "inverted"
} else if(type == "revert"){
newVals <- rev(vals)
action <- "reverted"
} else if(type == "descending"){
newVals <- sort(vals, decreasing = TRUE)
action <- "sorted into descending order"
} else if(type == "ascending"){
newVals <- sort(vals, decreasing = FALSE)
action <- "sorted into ascending order"
} else if(type == "cycle"){
newVals <- vals + by
newVals[newVals > max(vals, na.rm = TRUE)] <- newVals[newVals > max(vals, na.rm = TRUE)] - max(vals, na.rm = TRUE)
newVals[newVals < min(vals, na.rm = TRUE)] <- newVals[newVals < min(vals, na.rm = TRUE)] + max(vals, na.rm = TRUE)
action <- paste0("cycled by ", by)
}
temp <- subNumNumC(mat = mat, replace = vals, with = newVals)
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("values have been ", action)))
names(out) <- "values_permuted"
return(out)
}
#' Rescale the range of values of a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param range [\code{integer(2)}]\cr vector of minimum and maximum value to
#' which the values will be scaled.
#' @return A \code{RasterLayer} of the same dimension as \code{obj}, where the
#' cell values have been rescaled.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$continuous
#' scaled <- rRange(input, range = c(0, 25))
#' visualise(raster::stack(input, scaled))
#' @importFrom checkmate assertClass
#' @importFrom raster extent<- crs crs<-
#' @export
rRange <- function(obj, range = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertIntegerish(range, len = 2)
mat <- as.matrix(obj)
# scale to range
temp <- (mat - min(mat, na.rm = TRUE)) * (range[2] - range[1]) / (max(mat, na.rm = TRUE) - min(mat, na.rm = TRUE)) + range[1]
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the values have been scaled between ", range[1], " and ", range[2])))
names(out) <- "scaled"
return(out)
}
#' Combine a raster stack after segregation
#'
#' Successively combine the layers of a raster stack with each other (similar to
#' \code{\link{Reduce}}).
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param by [\code{list(.)}]\cr the cell values by which the layers shall be
#' aggregated; see Details.
#' @param fun [\code{function(1)}]\cr find the resulting value of the combined
#' layers.
#' @param weights [\code{numeric(length(obj))}]\cr weight by which the values in
#' each layer are multiplied before applying \code{fun}.
#' @param direction [\code{character(1)}]\cr the direction into which the
#' \code{RasterLayer} in \code{obj} is supposed to be combine. Either
#' \code{"right"} (default) or \code{"left"}; see Details.
#' @details The argument \code{direction} takes the direction into which the
#' layers should be combined. \code{"right"} means that layers are combined
#' from left to right. \code{rReduce} is based on the functional
#' \code{\link{Reduce}}, where this wording is handled differently.
#'
#' The number of layers in the aggregated raster depends on the length of the
#' list in \code{by}. If by is left empty, everything is written into one
#' \code{RasterLayer} object. Values to be aggregated in a \code{RasterLayer}
#' are in the same list element.
#' @return a \code{RasterLayer} or \code{RasterStack} of the same dimensions as
#' \code{obj}, in which the layers of \code{obj} are aggregated into a smaller
#' set of layers.
#' @family operators to modify a raster
#' @examples
#' input <- rtRasters$continuous
#' patches <- rPatches(rBinarise(input, thresh = 30))
#' myPatches <- rSegregate(patches)
#'
#' # revert the segregation
#' visualise(rReduce(myPatches))
#'
#' # group patches
#' twoGroups <- list(c(1:14), c(15:28))
#' visualise(rReduce(myPatches, by = twoGroups))
#'
#' # select a subset of patches
#' someLayers <- list(c(1, 3, 5, 7, 9))
#' visualise(rReduce(myPatches, by = someLayers))
#' @importFrom checkmate assertClass assertList assertFunction assertIntegerish
#' @importFrom raster as.matrix addLayer overlay brick extent<- crs crs<-
#' @export
rReduce <- function(obj, by = NULL, fun = function(x) sum(x, na.rm = TRUE),
weights = NULL, direction = "right"){
# check arguments
assertClass(obj, "RasterStackBrick")
assertList(by, types = "integerish", min.len = 1, null.ok = TRUE)
assertFunction(fun)
assertIntegerish(weights, null.ok = TRUE)
direction <- match.arg(direction, c("right", "left"))
if(is.null(weights)){
weights <- rep_len(1, dim(obj)[3])
} else{
if(length(weights) != dim(obj)[3]){
weights <- rep_len(weights, dim(obj)[3])
}
}
if(direction == "right"){
right <- FALSE
} else{
right <- TRUE
}
if(!is.null(by)){
out <- brick()
for(i in seq_along(by)){
tempBy <- by[i]
tempName <- names(tempBy)
tempObj <- obj[[tempBy[[1]]]]
tempColTab <- tempObj[[1]]@legend@colortable
tempWeights <- weights[by[[i]]]
theRasters <- lapply(X = 1:dim(tempObj)[3], function(x){
tempObj[[x]] * tempWeights[x]
})
tempOut <- Reduce(f = function(x, y) overlay(x, y, fun = fun),
x = theRasters,
right = right)
colortable(tempOut) <- tempColTab
names(tempOut) <- tempName
out <- addLayer(out, tempOut)
}
action <- paste0(length(by), " new layers")
} else{
theRasters <- lapply(X = 1:dim(obj)[3], function(x){
obj[[x]] * weights[x]
})
out <- Reduce(f = function(x, y) overlay(x, y, fun = fun),
x = theRasters,
right = right)
# I keep this here because at some point I might implement '#include <boost/...>' in reduceMatrixC
# theMatrixes <- lapply(1:dim(obj)[3], function(x) as.matrix(obj[[x]]))
# out <- reduceMatrixC(theMatrixes, f = fun)
action <- paste0("1 new layer")
names(out) <- "reduced"
}
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
extent(out) <- extent(obj)
crs(out) <- crs(obj)
out@history <- c(history, list(paste0("layers have been reduced into ", action)))
return(out)
}
#' Change the resolution of a raster
#'
#' Rescaling a raster changes the raster's resolution.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param factor [\code{integer(1)} | \code{numeric(1)}]\cr an integer for
#' up-scaling and a fraction for down-scaling (see Details).
#' @param kernelFunction [\code{kernelFunction}]\cr A kernel function to
#' determine new values of the rescaled raster, defaulting to
#' \code{boxKernel()}; see Details.
#' @return A \code{RasterLayer} with a dimension that has been up or down-scaled
#' from \code{obj}, where the new values have been determined by applying
#' \code{kernelFunction}.
#' @details The factor is a value by which the number of cells in x and y
#' dimension will be multiplied. If the raster should be up-scaled
#' (\code{factor > 1}), only integer values are allowed (will be truncated if
#' it is not an integer). If the raster should be down-scaled (\code{factor <
#' 1}) any inverse of an integer value (1/n) is allowed (denominator will be
#' rounded to the next integer, such as 1/2 or 1/5).
#'
#' The \code{kernelFunction} used here deviates from the other
#' (morpholological) kernels in that it is a function, which interpolates the
#' values of new cells. See \code{\link[mmand]{kernelFunction}} for an
#' in-depth description.
#' @family operators to modify a raster
#' @examples
#' input <- rtRasters$continuous
#' visualise(input)
#' visualise(rRescale(input, factor = 0.5))
#'
#' # up-scaling can be useful for follow-up morphological operations
#' binarised <- rBinarise(input, thresh = 30)
#' visualise(rSkeletonise(binarised))
#'
#' rescaled <- rRescale(binarised, factor = 2)
#' visualise(rSkeletonise(rescaled), new = TRUE)
#' @importFrom checkmate assertClass
#' @importFrom raster as.matrix raster extent extent<- crs crs<-
#' @importFrom mmand boxKernel rescale
#' @export
rRescale <- function(obj, factor = NULL, kernelFunction = NULL){
# check arguments
assertClass(obj, "RasterLayer")
assertNumeric(factor, any.missing = FALSE, len = 1)
if(is.null(kernelFunction)){
kernelFunction <- boxKernel()
} else{
assertClass(kernelFunction, "kernelFunction")
}
mat <- as.matrix(obj)
if(factor > 1){
factor <- trunc(factor)
} else{
denominator <- round(1/factor)
factor <- 1/denominator
}
temp <- rescale(mat, factor, kernelFunction)
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the raster has been rescaled by the factor ", factor)))
names(out) <- "rescaled"
# manage the bibliography entry (mmand)
bib <- bibentry(bibtype = "Manual",
title = "mmand: Mathematical Morphology in Any Number of Dimensions",
author = person("Jon", "Clayden"),
year = "2017",
note = "R package version 1.5.3",
url = "https://CRAN.R-project.org/package=mmand"
)
if(is.null(getOption("bibliography"))){
options(bibliography = bib)
} else{
currentBib <- getOption("bibliography")
if(!bib%in%currentBib){
options(bibliography = c(currentBib, bib))
}
}
return(out)
}
#' Segregate values of a raster into layers
#'
#' Distinct values in a raster will be assigned to layers in a raster stack.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param by [\code{RasterLayer(1)}]\cr optional object by which \code{obj}
#' should be segregated. If left empty, the distinct values of \code{obj} will
#' be taken.
#' @param flatten [\code{logical(1)}]\cr should all values be set to value 1
#' (\code{TRUE}) or should the original \code{obj} values be retained
#' (\code{FALSE}, default)?
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return a \code{RasterStack} of the same dimensions as \code{obj}, in which
#' the elements specified in \code{by} or the distinct values of \code{obj}
#' have each been assigned to a seperate layer.
#' @family operators to modify a raster
#' @examples
#' input <- rtRasters$continuous
#' patches <- rPatches(rBinarise(input, thresh = 30), background = 0)
#' myPatches <- rSegregate(patches)
#' visualise(myPatches[[c(2, 3, 12, 16)]])
#'
#' # when flattening, all values are set to 1
#' myPatches2 <- rSegregate(patches, flatten = TRUE)
#' visualise(myPatches2[[c(2, 3, 12, 16)]])
#'
#' # cut out by 'patches'
#' patchValues <- rSegregate(input, by = patches)
#' visualise(patchValues[[c(2, 3, 12, 16)]])
#' @importFrom checkmate assertClass testNull assert testClass assertLogical
#' assertIntegerish
#' @importFrom raster as.matrix values extent<- stack crs crs<-
#' @importFrom mmand shapeKernel components
#' @export
rSegregate <- function(obj, by = NULL, flatten = FALSE, background = NULL){
# check arguments
objMeta <- checkGDAL(x = obj)
byMeta <- checkGDAL(x = by)
assertLogical(flatten)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
# get provenance
history <- getProv(x = obj)
if(!objMeta$isGDAL & !byMeta$isGDAL){
mat <- as.matrix(obj)
if(!byMeta$isRaster){
subMat <- mat
} else {
subMat <- as.matrix(by)
}
vals <- unique(as.vector(subMat))
vals <- vals[!is.na(vals)]
theRasters <- list()
for(i in seq_along(vals)){
tempMat <- mat
tempMat[subMat != vals[i] | is.na(subMat)] <- background
if(flatten){
tempMat[subMat == vals[i]] <- 1
}
tempRaster <- raster(tempMat)
extent(tempRaster) <- extent(obj)
crs(tempRaster) <- crs(obj)
tempRaster@history <- c(history, list(paste0("element ", vals[i], " has been segregated from the main raster")))
theRasters <- c(theRasters, list(tempRaster))
}
theRasters <- stack(theRasters)
theRasters@history <- c(history, list(paste0("the raster has been segregated")))
names(theRasters) <- paste0("element_", vals)
} else {
theRasters <- gdal_segregate(src = objMeta$path, by = byMeta$path,
output_Raster = TRUE)
if(objMeta$isGDAL){
}
if(byMeta$isGDAL){
# save obj as tif and proceed with gdal_segregate
}
}
# manage the bibliography entry (mmand)
setBib(bib = bibentry(bibtype = "Manual",
title = "mmand: Mathematical Morphology in Any Number of Dimensions",
author = person("Jon", "Clayden"),
year = "2017",
note = "R package version 1.5.3",
url = "https://CRAN.R-project.org/package=mmand"
))
return(theRasters)
}
#' Determine the skeleton of foreground patches in a raster
#'
#' The morphological skeleton of a patch is a binary skeletal remnant that
#' preserves the extent and connectivity (i.e. the topology) of the patch.
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param method [\code{character(1)}]\cr the method to determine the skeleton.
#' Either \code{"hitormiss"} (default), \code{"lantuejoul"} or
#' \code{"beucher"}.
#' @param kernel [\code{matrix(1)}]\cr scan the raster with this kernel (default
#' is a 3 by 3 cells diamond kernel).
#' @param background [\code{integerish(1)}]\cr the value any cell with value NA
#' should have.
#' @return a \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' foreground patches have been transformed into their morphological
#' skeletons.
#' @seealso \code{\link{rDistance}}
#' @details For details, refer to \code{\link[mmand]{skeletonise}} of which
#' \code{rSkeletonise} is a wrapper.
#' @family operators to determine objects
#' @examples
#' input <- rtRasters$continuous
#' binarised <- rBinarise(input, thresh = 30)
#'
#' visualise(raster::stack(binarised, rSkeletonise(binarised)))
#' @importFrom checkmate assertClass
#' @importFrom mmand shapeKernel skeletonise
#' @importFrom raster as.matrix raster extent extent<- crs crs<-
#' @export
rSkeletonise <- function(obj, method = "hitormiss", kernel = NULL, background = NULL){
# check arguments
assertClass(obj, "RasterLayer")
method <- match.arg(arg = method, choices = c("lantuejoul", "beucher", "hitormiss"))
assertMatrix(kernel, mode = "integerish", all.missing = FALSE, null.ok = TRUE)
if(is.null(kernel)) kernel <- matrix(c(0, 1, 0, 1, 1, 1, 0, 1, 0), 3, 3)
assertIntegerish(background, null.ok = TRUE)
if(is.null(background)){
background <- NA
}
if(is.null(kernel)){
kernel <- shapeKernel(c(3, 3), type="diamond")
}
mat <- as.matrix(obj)
if(!isBinaryC(mat)){
stop("spatial object is not binary, please run 'rBinarise()' initially.")
}
temp <- skeletonise(mat, kernel, method = method)
# set na
temp[temp == 0] <- background
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("the morphological skeleton has been determined")))
names(out) <- "skeletonised"
# manage the bibliography entry (mmand)
bib <- bibentry(bibtype = "Manual",
title = "mmand: Mathematical Morphology in Any Number of Dimensions",
author = person("Jon", "Clayden"),
year = "2017",
note = "R package version 1.5.3",
url = "https://CRAN.R-project.org/package=mmand"
)
if(is.null(getOption("bibliography"))){
options(bibliography = bib)
} else{
currentBib <- getOption("bibliography")
if(!bib%in%currentBib){
options(bibliography = c(currentBib, bib))
}
}
return(out)
}
#' Substitute values in a raster
#'
#' @param obj [\code{RasterLayer(1)}]\cr The object to modify.
#' @param old [\code{integerish(.)}]\cr values to be substituted.
#' @param new [\code{integerish(1)}]\cr value to substitute with.
#' @param col [\code{character(1)}]\cr six-digit hex colour code.
#' @return A \code{RasterLayer} of the same dimensions as \code{obj}, in which
#' \code{old} values have been replaced by the \code{new} value.
#' @family operators to modify cell values
#' @examples
#' input <- rtRasters$categorical
#' substituted <- rSubstitute(input, old = c(41, 44, 47), new = list(id = 40))
#' visualise(raster::stack(input, substituted))
#' @importFrom checkmate assertClass
#' @importFrom tibble as_tibble
#' @importFrom raster values as.matrix extent<- crs crs<- ratify
#' @export
rSubstitute <- function(obj, old = NULL, new = NULL, col = "#FF55FF"){
# check arguments
assertClass(obj, "RasterLayer")
isFactor <- obj@data@isfactor
assertIntegerish(x = old, any.missing = FALSE, min.len = 1)
assertList(x = new, min.len = 1)
assertNames(x = names(new), must.include = "id")
assertCharacter(x = col, pattern = "^\\#[a-zA-Z0-9]{6}$", null.ok = TRUE)
if(length(old)!=length(new$id)){
newValues <- rep(new$id, length.out = length(old))
} else{
newValues <- new$id
}
mat <- as.matrix(obj)
rat <- levels(obj)[[1]]
cltbl <- colortable(obj)
if(!is.null(col)){
cltbl[new$id] <- col
}
temp <- subNumNumC(mat = mat, replace = old, with = newValues)
out <- raster(temp)
extent(out) <- extent(obj)
crs(out) <- crs(obj)
if(!is.null(rat)){
outRat <- rat[-which(rat$id %in% old),]
outRat <- rbind(outRat, data.frame(new))
out <- ratify(out)
out@data@attributes[[1]] <- outRat
}
# assign colourtable
if(!length(cltbl) == 0){
colortable(out) <- cltbl
}
# assign history
if(length(obj@history)==0){
history <- list(paste0("the object was loaded from memory"))
} else{
history <- obj@history
}
out@history <- c(history, list(paste0("value of ", old, " replaced with ", new)))
names(out) <- "values_substituted"
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.