1

So essentially, what I want to do is replace every element in a matrix with the maximum of neighboring cells within a window that is determined by the value in that cell.

The window size would be determined by this function (fitlwr), where Tree_Height calls a linear model that was fit to a dataset of Tree Height and Crown Diameter data:

RoundOdd <- function(x) {2*floor(x/2)+1} #makes sure window size is an odd number

fitlwr <- function(x){for(i in x){
  if(i > 13){
    m <- RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == i, "fit.lwr"]) 
  return(matrix(1, nrow = m, ncol = m))
    }
  else {
    return(matrix(1, 3, 3))
    }
}}

I then want to replace every value in that matrix with the maximum of the values within that window.The matrix was derived from a raster layer and the values represent the height above ground for a given cell. The dimensions are 6,571 x 5,764. A section of the data might look like this:

    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    9   47  103   58   80   55   72   56   14    52
 [2,]   68   49   49   43   62   80   62   23   55    82
 [3,]   58   10   79   70   75   49   68   60   74    79
 [4,]   78   19   51   26   61   77   57   70   51    43
 [5,]   47   88   57   80   25   33   24   30   56    63
 [6,]   73   36   53   25   63   30   19   59   17    63
 [7,]   95    9   49   95    6   13   21   75   60    34
 [8,]   36   65   47   64   22   66   52    9   71    20
 [9,]   45   53   31   47  114   55   44   42   44    44
[10,]   47   23  102   34   67   60    5   23   61    32

The raster focal functions were my go-to, but they don't let you use a variable window size (see below).

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

I am open to using another language or software tools to accomplish this task, including GRASS, Python, QGIS, or ArcGIS if necessary.

1
  • The GRASS module r.neighbors calculates pixels values by applying some function to a user defined neighborhood, such as mean, max, min, etc. No need to reinvent the wheel.
    – Micha
    Oct 31 '20 at 7:42
0

Okay! So I think I figured it out:

First I convert my raster to a matrix, remove 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.