rspatial / raster

R raster package https://rspatial.github.io/raster/reference/raster-package.html
GNU General Public License v3.0
161 stars 53 forks source link

Using Variable Window Sizes within Focal Function #176

Open azh3 opened 3 years ago

azh3 commented 3 years ago

So I am attempting to implement a maximum filter on a CHM (Canopy Height Model) with a variable window size here, 'Tree_Heights' is a data frame containing the fitted values of a linear regression model to predict tree diameter, the window size is supposed to match the lower prediction interval for the radius of a tree at a given height.

RoundOdd <- function(x) {2*floor(x/2)+1}

fitlwr <- function(x){ RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"]/2) }

m <- raster::focalWeight(x = CMM, d = fitlwr(), type = "circle")

CMM <- raster::focal(x = CMM, w = m, fun = max)

This returns the following error:

Error in [.data.frame(Tree_Heights, Tree_Heights$Tree_Height == x, "fit.lwr") : argument "x" is missing, with no default 6.[.data.frame(Tree_Heights, Tree_Heights$Tree_Height == x, "fit.lwr") 5.Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"] 4.RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"]/2) 3.fitlwr() 2..circular.weight(x, d[1]) 1.raster::focalWeight(x = CMM, d = fitlwr(), type = "circle")

If I try instead to use the function in the argument for window size, I get this error:

Error in .local(x, ...) : is.matrix(w) is not TRUE

  1. stop(simpleError(msg, call = if (p <- sys.parent(1L)) sys.call(p)))
  2. stopifnot(is.matrix(w))
  3. .local(x, ...)
  4. raster::focal(x = CMM, w = fitlwr, fun = max)
  5. raster::focal(x = CMM, w = fitlwr, fun = max)

Is there no way to use a variable window size inherently in this function? Are there any other options for doing this in the Raster package?

azh3 commented 3 years ago

I think I found a temporary

First I converted my raster to a matrix, removed NAs, and round the values (this is unique to my use case and not necessary for the algorithm to function.

    X <- raster::as.matrix(Z)
    X <- round(X, digits = 0)
    X[is.na(X)] <- 0

This is to calculate the maximum for a variable size rectangular moving window:

  Y <- X
    for (i in 1:nrow(X)){ 
       for (j in 1:ncol(X)){ 
          N <- fitlwr(X[i,j])
          Y[i,j] = max(X[max(1, i-N):min(nrow(X), i+N), max(1, j-N):min(ncol(X), j+N)]) 
      }
    }

fitlwr() is a custom function that calls a linear model that matches the value of a cell to the expected radius of the moving window.

And here is for a circular moving window:

    Y <- X for (i in 1:nrow(X)){     for (j in 1:ncol(X)){ 
           N = fitlwr(X[i,j])
           M = X[max(1, i-N):min(nrow(X), i+N), max(1, j-N):min(ncol(X), j+N)]
           W = reshape2::melt(M)
           W$d2 = sqrt((W$Var1-mean(W$Var1))^2 + (W$Var2-mean(W$Var2))^2)
           Y[i,j] = max(X[i,j], max(subset(W, d2 <= N, select = value)))}

I then write the values back to my raster, this maintains the CRS, projection, etc.

raster::values(Z) <- Y

I imagine this would be more efficient if it could be implemented within the raster package itself and think there are many potential use-cases for being able to specify a variable window size within the focal function.