Nothing
#' @name aoi
#' @title Build area of interest (A.O.I.)
#'
#' @description This function creates a spatial object (raster or vector) defining the area of interest, by taking a bounding box or
#' a spatial object, or unioning the input spatial objects if more than one are provided.
#' When \code{msk} is a list of rasters, extent is set equal to their combined extent (union) and resolution
#' to the finest resolution among them.
#' @param msk a character (path to raster or vector file), or a bounding box as numeric (xmin,xmax,ymin,ymax),
#' or one or more (as list of) rasters of class "RasterLayer", or a single object of class "sf" or "SpatialPolygonsDataFrame".
#' The reference data (raster or vector) to be used as mask. All model outputs will have the same extent (outline) as this object.
#' All locations with no data (i.e. NA) cells in \code{msk} input will be ignored as well.
#' @param mskSub vector of values, for raster data only. The subset values from \code{msk} which should be considered to build the area
#' of interest. All other values will be ignored and masked out during modelling.
#' @param xy logical. Should return a two column matrix of x and y coordinates of cells centre (raster data) or
#' the IDs of features? Defaults to FALSE, returning an object of class RasterLayer or sf.
#' @param bbox numeric of four elements, the coordinates defining a rectangle (bounding box) to limit the area of interest.
#' Must be ordered as xmin, xmax, ymin, ymax. Coordinates must be in the same reference system as spatial data.
#' @details If rasters are used, all model outputs will have the same resolution and same extent as inherited from \code{msk}.
#' All locations with no data (i.e. NA) cells from \code{msk} will be ignored as well.
#' @return An object of class RasterLayer or sf, or a matrix of coordinates of mask cells (raster only). In the former case, valid cells
#' (i.e. the area of interest) will have value 1, \code{NA} otherwise.
#' @seealso \code{\link{extractByMask}}
#' @examples
#' ## Make a mask from a group of input layers:
#' list2env(ConwyData, environment())
#'
#' network <- LandUseChange
#' spatialData <- c(ConwyLU, ConwySlope, ConwyStatus)
#' m <- aoi(spatialData)
#' m
#'
#' ## Plot mask
#' library(raster)
#' m <- aoi(ConwyLU)
#' plot(m)
#'
#' ## Make mask from a subset of values and plot
#' m <- aoi(ConwyLU, mskSub=c(2,3))
#' plot(m)
#'
#' ## Return coordinates of valid mask locations
#' coord <- aoi(ConwyLU, xy=TRUE)
#' head(coord)
#'
#' ## Using a bounding box
#' aoi(ConwyLU, bbox=c(270000, 284950, 347000, 365000))
#'
#' ## For vectorial spatial data. Note xy=TRUE shall return the features IDs
#' Conwy = sf::st_read(system.file("extdata", "Conwy.shp", package = "bnspatial"))
#' aoi(Conwy, bbox=c(270000, 284950, 347000, 365000))
#'
#' @export
aoi <- function(msk, mskSub=NULL, xy=FALSE, bbox=NULL){ ## Check if aoi and extractByMask can be condensed in one or nested.
msk <- .loadSpatial(msk, checkfld=FALSE)
msk <- .bbox(bbox, msk)
if('RasterLayer' %in% class(msk[[1]])){
.aoiRaster(msk, mskSub, xy)
} else {
if(!is.null(mskSub)) warning('mskSub argument ignored for vectorial spatial data')
.aoiVector(msk, xy)
}
}
##
.bbox <- function(bbox, msk){
if(!is.null(bbox)){
if(!is.numeric(bbox) | length(bbox) != 4) stop('Bounding box must be a numeric vector of length 4 (i.e xmin, xmax, ymin, ymax)')
if('RasterLayer' %in% class(msk[[1]])){
if(is.list(msk)){
msk <- lapply(msk, function(m) { raster::crop(m, bbox) })
} else {
msk <- raster::crop(msk, bbox)
}
} else {
msk <- sf::st_crop(sf::st_buffer(msk, dist = 0),
xmin=bbox[1],xmax=bbox[2],ymin=bbox[3],ymax=bbox[4])
}
}
return(msk)
}
##
.aoiRaster <- function(msk, mskSub, xy){
# if(class(msk) == 'character'){
# msk <- lapply(msk, function(x) { raster::raster(x, RAT=FALSE) }) ## Should apply this even when an input raster is provided
# }
if(is.list(msk)){
.checkSR(msk)
r <- msk[[1]]
ext <- raster::extent(r)
cellSizeX <- raster::res(r)[1]
cellSizeY <- raster::res(r)[2]
if(length(msk) > 1){
for(i in 2:length(msk)){
r <- msk[[i]]
ext <- raster::union(ext, raster::extent(r) )
cellSizeX <- ifelse(cellSizeX > raster::res(r)[1], raster::res(r)[1], cellSizeX)
cellSizeY <- ifelse(cellSizeY > raster::res(r)[2], raster::res(r)[2], cellSizeY)
}
}
msk <- raster::raster(ext=ext, resolution=c(cellSizeX, cellSizeY), crs=raster::crs(r), vals = 1)
if(xy == TRUE){
return( raster::xyFromCell(msk, seq_along(msk)) )
} else {
return( msk )
}
# TODO cope with lower resolution raster
} else {
mskVals <- raster::getValues(msk)
if(!is.null(mskSub)){
id = mskVals %in% mskSub
} else {
id = is.finite(mskVals)
}
if(xy == TRUE) {
id = seq_along(msk)[id]
return( raster::xyFromCell(msk, id) )
} else {
msk[] <- NA
msk[id] <- 1
names(msk) <- 'layer'
return( msk )
}
}
}
##
.aoiVector <- function(msk, xy=FALSE){ # use xy to return id?
if(!'FID' %in% names(msk)) {
# warning('No features ID found in attributes of "msk" (FID column); ',
# 'a column named "FID" will be guessed from ',
# 'attributes table, with unique integers for each feature.')
msk["FID"] <- 1:nrow(msk)
}
if(xy){
return(msk[['FID']])
} else {
return(msk['FID'])
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.