R/vwf.R

Defines functions vwf

Documented in vwf

#' Variable Window Filter
#'
#' Implements the variable window filter algorithm (Popescu & Wynne, 2004) for detecting treetops from a canopy height model.
#'
#' This function uses the resolution of the raster to figure out how many cells the window needs to cover.
#' This means that the raster value (representing height above ground) and the map unit (represented by the raster's resolution),
#' need to be in the _same unit_. This can cause issues if the raster is in lat/long, whereby its resolution is in decimal degrees.
#'
#' @param CHM Canopy height model. Either in \link[raster]{raster} format, or a path directing to a raster file. A character vector of multiple paths directing to a
#' tiled raster dataset can also be used.
#' @param winFun function. The function that determines the size of the window at any given location on the
#' canopy. It should take the value of a given \code{CHM} pixel as its only argument, and return the desired *radius* of
#' the circular search window when centered on that pixel. Size of the window is in map units.
#' @param minHeight numeric. The minimum height value for a \code{CHM} pixel to be considered as a potential treetop. All \code{CHM} pixels beneath
#' this value will be masked out.
#' @param maxWinDiameter numeric. Sets a cap on the maximum window diameter (in cells). If an
#' improperly calibrated function is set for \code{winFun}, it may produce overly large windows that would perform poorly
#' and significantly slow processing time. This setting can be disabled by setting to \code{NULL}.
#' @param minWinNeib character. Define whether the smallest possible search window (3x3) should use a \code{queen} or
#' a \code{rook} neighborhood.
#' @param verbose logical. Print progress to console if set to \code{TRUE}.
#'
#' @references Popescu, S. C., & Wynne, R. H. (2004). Seeing the trees in the forest. \emph{Photogrammetric Engineering & Remote Sensing, 70}(5), 589-604.
#'
#' @return \link[sp:SpatialPoints]{SpatialPointsDataFrame}. The point locations of detected treetops. The object contains two fields in its
#' data table: \emph{height} is the height of the tree, as extracted from the \code{CHM}, and \emph{winRadius} is the radius
#' of the search window when the treetop was detected. Note that \emph{winRadius} does not necessarily correspond to the radius
#' of the tree's crown.
#'
#' @examples
#' # Set function for determining variable window radius
#' winFunction <- function(x){x * 0.06 + 0.5}
#'
#' # Set minimum tree height (treetops below this height will not be detected)
#' minHgt <- 2
#'
#' # Detect treetops in demo canopy height model
#' ttops <- vwf(CHMdemo, winFunction, minHgt)
#'
#' @seealso \code{\link{mcws}} \code{\link{sp_summarise}}
#'
#' @export

vwf <- function(CHM, winFun, minHeight = NULL, maxWinDiameter = 99, minWinNeib = "queen", verbose = FALSE){

  ### CHECK INPUTS

    if(verbose) cat("Checking inputs", "\n")

    # Check for valid inputs for 'minWinNeib'
    if(!minWinNeib %in% c("queen", "rook")) stop("Invalid input for 'minWinNeib'. Set to 'queen' or 'rook'")

    # Check for unprojected rasters
    CHM.crs <- as.character(raster::crs(CHM))
    CHM.prj <- regmatches(CHM.crs, regexpr("(?<=proj=).*?(?=\\s)", CHM.crs, perl = TRUE))
    if(length(CHM.prj) > 0 && CHM.prj %in% c(c("latlong", "latlon", "longlat", "lonlat"))){
      warning("'CHM' map units are in degrees. Ensure that 'winFun' outputs values in this unit.")
    }

    # Round out CHM resolution to fifth decimal and check that CHM has square cells.
    # Rounding is necessary since a lack of precision in CHM cell size call cause the
    # 'focalWeight' function to misbehave
    roundRes <- round(raster::res(CHM), 5)
    if(roundRes[1] != roundRes[2]) stop("Input 'CHM' does not have square cells")
    if(roundRes[1] == 0)           stop("The map units of the 'CHM' are too small")

    # Ensure that 'minHeight' argument is given a positive value
    if(!is.null(minHeight) && minHeight <= 0) stop("Minimum canopy height must be set to a positive value.")

    # Get range of CHM values
    CHM.rng <- if(CHM@data@haveminmax){
      c(CHM@data@min, CHM@data@max)
    }else{
      suppressWarnings(raster::cellStats(CHM, "range"))
    }
    names(CHM.rng) <- c("min", "max")

    # Check if CHM has usable values
    if(is.infinite(CHM.rng["max"]) | is.infinite(CHM.rng["min"])){stop("Input 'CHM' does not contain any usable values. Check input data.")}


  ### APPLY MINIMUM CANOPY HEIGHT ----

    if(!is.null(minHeight)){

      if(minHeight >= CHM.rng["max"]) stop("'minHeight' is set to a value higher than the highest cell value in 'CHM'")

      # Mask sections of CHM that are lower than 'minHeight'
      if(minHeight > CHM.rng["min"]){

        CHM[CHM < minHeight] <- NA
        CHM.rng["min"] <- minHeight

      }
    }


  ### CREATE WINDOWS ----

    # Here, the variably sized windows used for detecting trees are "pre-generated". First, a series of 'winRadii'
    # is generated, representing all the sizes the windows can take based on the range of potential values returned
    # by 'winFun'. These radii are then converted to binary matrices, where values of 1 represent the circular shape
    # of each window. These matrices are then converted to vectors and

    if(verbose) cat("Creating windows", "\n")

    # Generate a list of radii
    seqFloor   <- APfun::AProunder(winFun(CHM.rng["min"]), interval = roundRes[1], direction = "down")
    seqCeiling <- APfun::AProunder(winFun(CHM.rng["max"]), interval = roundRes[1], direction = "up")
    if(is.infinite(seqFloor)) seqFloor <- 0 # Watch out for parabola!
    winRadii <- seq(seqFloor, seqCeiling, by = roundRes[1])

    # Remove radii that are smaller than the CHM's resolution
    winRadii <- winRadii[winRadii >= roundRes[1]]
    if(length(winRadii) == 0){
      warning("The maximum window radius computed with 'winFun' is smaller than the CHM's resolution",
              "\nA 3x3 cell search window will be uniformly applied",
              "\nUse a higher resolution 'CHM' or adjust 'winFun' to produce wider dynamic windows")
      winRadii <- roundRes[1]
    }

    # Calculate the dimensions of the largest matrix to be created from the generated list of radii

    # NOTE: 'winDiameter' should be an uneven integer. This was causing a problem until Michael Koontz
    # pointed out the issue on 2018/11/15. Having the sequence of 'winRadii' created using AProunder should
    # solve this issue but just a failsafe, 'ceiling' is used to force the number into an integer and
    # an explicit check for even numbers was added
    winDiameter <- ceiling((max(winRadii) / roundRes[1]) * 2)
    if (winDiameter %% 2 == 0) winDiameter <- winDiameter + 1

    # Check if input formula will yield a window size bigger than the maximum set by 'maxWinDiameter'
    if(!is.null(maxWinDiameter) && winDiameter > maxWinDiameter){

        stop("Input function for 'winFun' yields a window diameter of ",  winDiameter, " cells, which is wider than the maximum allowable value in \'maxWinDiameter\'.",
             "\nChange the 'winFun' function or set 'maxWinDiameter' to a higher value (or to NULL).")
    }

    # Convert radii into windows
    windows <- lapply(winRadii, function(radius){

      # Based on the unit size of the input CHM and a given radius, this function will create a matrix whose non-zero
      # values will form the shape of a circular window
      win.mat <- raster::focalWeight(raster::raster(resolution = roundRes), radius, type = "circle")

      # Apply Queen's neighborhood if circle is 3x3
      if(nrow(win.mat) == 3 && minWinNeib == "queen") win.mat[] <- 1

      # Pad the window to the size of the biggest matrix created from the list of radii
      win.pad <- raster::extend(raster::raster(win.mat), (winDiameter - ncol(win.mat)) /2, value = 0)

      # The matrix values are then transformed into a vector
      win.vec <- as.vector(win.pad != 0)

      return(win.vec)

    })
    names(windows) <- winRadii


  ### CREATE VWF FUNCTION ----

    .vwMax <- function(x, ...){

      # Locate central value in the moving window.
      centralValue <- x[length(x) / 2 + 0.5]

      # If central value is NA, then return NA.
      if(is.na(centralValue)){

        return(NA)

      }else{

        # Calculate the expected crown radius.
        radius <- winFun(centralValue)

        # Retrieve windows size closest to radius
        window <- windows[[which.min(abs(as.numeric(names(windows)) -  radius))]]

        # If the central value is the highest value within the variably-sized window (i.e.: local maxima), return 1. If not, return 0.
        return(if(max(x[window], na.rm = TRUE) == centralValue) 1 else 0)

      }
    }


  ### APPLY VWF FUNCTION ----

    if(verbose) cat("Detecting local maxima", "\n")

    # Apply local maxima-finding function to raster
    localMax <- raster::rasterToPoints(

      raster::focal(CHM, matrix(1, winDiameter, winDiameter), .vwMax,
        pad = TRUE, padValue = NA),

      fun = function(x) x == 1, spatial = TRUE)

    # Add attributes
    localMax[["height"]]    <- raster::extract(CHM, localMax)
    localMax[["winRadius"]] <- winFun(localMax[["height"]])
    localMax[["treeID"]]    <- 1:length(localMax)
    localMax[["layer"]] <- NULL


  ### RETURN OUTPUT ----

    return(localMax)

  }

Try the ForestTools package in your browser

Any scripts or data that you put into this service are public.

ForestTools documentation built on Sept. 11, 2021, 9:07 a.m.