3

I'm trying to normalize a DEM raster with focal and a function to set values between 0 and 1 (iteratingly from the neighboring cells). Thus, for each cell, a 3x3 matrix is to be created and the function applied.

While it works well with the built-in functions (min, max, mean...), I can't get it to work with a normalizing function, ie (x-min(x))/max(x)-min(x)).

Example:

# create a raster
x <- raster(ncol=3,nrow=3)
x[] <- c(13,14,23,46,25,18,32,39,33)
extent(x) <- extent(c(-1,1,-1,1))

#applying the function
res <- focal(x, w=matrix(1,3,3), fun = function(x){(x-min(x))/(max(x)-min(x))})

> res
class      : RasterLayer 
dimensions : 3, 3, 9  (nrow, ncol, ncell)
resolution : 0.6666667, 0.6666667  (x, y)
extent     : -1, 1, -1, 1  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0, 0  (min, max)

I would expect the central cell to return a value of 0,3636... But it just gives me 0. Any idea what I'm doing wrong?

2
  • If I omit setting the extent, it works for me. I am adding pad=TRUE and padValue=1 as additional argument, just to keep the original dimensions. The formula that you are using is in accounting for potential negative numbers. It is very likely that you can just row standardize by [x / max(x)] which would simplify your focal function. However, this seems very dubious. You should be standardizing to a global max and not local. Feb 7 '20 at 19:18
  • Thanks for the quick answer... I tried with {x/max(x)} as you suggested and got 0.283... Which would be the value of first cell divided by the max, 13/46... Any idea how can I get the result for the value of the central cell instead, 25/46? Or am I doing something wrong?
    – gvanhavre
    Feb 7 '20 at 22:11
2

If the focal value is not specifically indexed in the function it looks like the focal function is assigning the incorrect values and positions back to the matrix. However, this is expected behavior. For these types of operations, you have to specifically index the focal value to get the correct results. One would think that the focal value would automatically be referenced as the value that is operated on so, simply doing x / max(x) would work. However, keep in mind that a vector is returned in within the focal function so your function needs to operate on the vector of values and not just assume a single value in relation to the other focal values. That is to say that a 3x3 window will return a vector of 9 values with the 5th value being the focal value. So an operator would look something like this x[5] / max(x).

In the following example, on a 3x3 matrix/raster and window the center matrix value should be 0.54. If I pass a function that uses x / max(x) the results are miss-assigned back to the raster. This is because the results are a vector with 9 values, the first being 0.2826087, which is what is assigned to the resulting focal value. In this case if it took the 5 value, (x / max(x))[5], the result would be correct.

library(raster)

r <- raster(matrix(c(13,14,23,46,25,18,32,39,33),3,3))
as.matrix(r[2,2])/max(as.matrix(r))
(r[] / max(r[]))[5]

f1 <- function(x,na.rm) { return(x/max(x)) }
as.matrix(focal(r, w=matrix(1,3,3), pad=TRUE, padValue=1, fun=f1))

However, if you explicitly index the focal value (which is the center value of the vector) the the correct results are returned.

f3 <- function(x,y=46,na.rm) { return(x[ceiling(length(x)/2)] / max(x)) }
as.matrix(focal(r, w=matrix(1,3,3), pad=TRUE, padValue=1, fun=f4))
3
  • I just found an interesting hint from another post. Adding the function > test <- function(x) {center <- x[ceiling(length(x)/2)] center/max(x)} gives me the correct values on the initial 3x3 raster.
    – gvanhavre
    Feb 8 '20 at 0:02
  • @gvanhavre, As I was looking at code last week I thought that somehow my logic was not tracking. Please see my updated answer. The focal function was behaving exactly as expected the problem is that you have to operate on the vector that is returned within the function and not assume a single value. Hopefully I am clear in how this is explained in my answer. If you have any questions please feel free to contact me. Feb 10 '20 at 17:20
  • Of course, now it makes perfect sense. I thought the focal was returning 9 separate values indeed. Thank you so much for all the care.
    – gvanhavre
    Feb 11 '20 at 10:28
0

Looks like your raster datatype may be set to integer.

Try: rasterOptions(datatype='FLT4S')

1
  • This doesn't work for the example in the question. Indeed, inspecting the results shows the datatype was double all along.
    – whuber
    Feb 7 '20 at 20:46

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.