# How to replace every element in a matrix or raster with the maximum of the values within a variable moving widow of size = n?

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.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.

• 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. Oct 31 '20 at 7:42

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
``````