0

I have issues understanding the operating of the terra::focal function. I created the following raster:

r<-rast(ncol=5,nrow=5,vals=1:25)
r_df <- as.data.frame(r, xy=TRUE)
colnames(r_df) <- c("x", "y", "z")

ggplot(r_df, aes(x=x, y=y, fill=z)) +
geom_tile() +
geom_text(label = r_df$z, size = 4, colour = "black") +
scale_fill_gradient2(low = "blue", high = "red",
                   limits = c(min(r_df$z), max(r_df$z)),
                   midpoint = median(r_df$z)) +
scale_x_discrete(expand = c(0,0)) +
scale_y_reverse() +
theme(panel.grid.minor = element_blank()) +
theme(panel.grid.major = element_blank()) +
theme(axis.ticks = element_blank()) +
theme(panel.background = element_rect(fill = "transparent"))+
theme(legend.position = "none") +
theme()

enter image description here

Now I want to calculate for each cell the sum of the neighboring cells values:

weight_matrix <- matrix(1, nrow=3, ncol=3)
weight_matrix[2,2] <- NA
r<-terra::focal(r,weight_matrix, fun=sum, na.policy="all", na.rm=TRUE, expand=FALSE, fillvalue=NA)

Which then yields the following matrix:

enter image description here

I do not understand how the values on the outside of the raster are calculated. Looking at the bottom left cell, I would expect the value to be 15 (2 + 6 + 7), however, it is equal to 30. It seems like the values from the bottom right are included in the sum. I played around with the parameters but don't think I get it right.

1 Answer 1

1

Your ggplot function is hiding the X axis values (and reversing the Y axis for some reason?). When I used a simpler function that didn't, I saw the problem.

plr <- function(r){
    rd = as.data.frame(r, xy=TRUE)
    plot(r)
    text(rd[,1],rd[,2], as.character(rd[,3]))
}

When I plot your r, its from -180 to 180, so it wraps around the sphere, and focal is respecting that by wrapping the weight matrix across the 180:-180 boundary:

enter image description here

As a demonstration, I made a matrix of all 1 values and smoothed it with your function call:

> r1 = rast(ncols=5, nrows=5, vals=1)
> s1 = terra::focal(r1,weight_matrix, fun=sum, na.policy="all", na.rm=TRUE, expand=FALSE, fillvalue=NA)
> plr(s1)

enter image description here

But if I change the CRS to 27700, so instead of being a -180 to 180 degree raster covering the globe, its a 360 x 180 metre tiny plot of ocean off the coast of the UK, I get a different result:

> crs(r1) = "epsg:27700"
> s1 = terra::focal(r1,weight_matrix, fun=sum, na.policy="all", na.rm=TRUE, expand=FALSE, fillvalue=NA)
> plr(s1)

enter image description here

If you want to stop if wrapping round the world, give it a missing CRS and replace the CRS back on the computed raster. I'm not sure why you wouldn't want a focal operation to not look round the back of the world but here's how to do that:

> r1 = rast(ncols=5, nrows=5, vals=1)
> # save original crs and set to nothing
> c1 = crs(r1) ; crs(r1) = "" 
> # do the focal now on a raster with no CRS so focal doesn't know the
> # earth is round
> s1 = terra::focal(r1,weight_matrix, fun=sum, na.policy="all", na.rm=TRUE, expand=FALSE, fillvalue=NA)
> # restore the non-flat earth
> crs(s1) = c1 ; crs(r1) = c1
> plot(s1)
1
  • Thanks, now I do understand it! Commented Mar 21 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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