5

I am trying to use R with the focal function from the raster package to smooth raster images.

I am working with Radar images that have already been preprocessed. As part of that, some calculations introduce unrealistic values that are then deleted from the raster. So basically, I am working with images with values ranging from -2 (black) to 7 (white) like this and want to smooth over the small datagaps: enter image description here

To do so, I tried the focal function from the raster package with the following code:

raster1 <- raster("example_raster.tif")

mat <- matrix(1/25,ncol=5, nrow=5)
raster_res <- focal(raster1, mat, FUN="mean", na.rm=T) 

writeRaster(raster_res, filename="smooth_raster.tif", format="GTiff")

The problem with this is, that the NAs seem to introduce some kind of edge effect. na.rm is set to TRUE to exclude the NAs from the actual calculations, but if you look at the edge at the east or how spotty the areas around former NA values appear, it looks like R still somehow includes them in the moving window calculations, somehow influencing the sourrounding values. Here is what the output looks like: enter image description here

To check, if this is a normal focal-problem, I also used r.neighbors in the QGIS processing toolbox set to a window size of 5 to smooth my data and it does exactly what I expected. Unfortunately I can't post a third link due to being new, but the raster is very smooth and does not have the spotty appearence of the output R produces.

I already tried using the pad and padValue arguments in focal but did not manage to make it work the way I want. Unfortunately, I don't know enough about GRASS to understand, what exactly r.neighbors does differently than focal. Since this will be part of a script that automatically processes several rasters, doing everything manually with QGIS is not really an option.

  • Have you tried not having quotation marks around 'mean'? – Mikkel Lydholm Rasmussen Mar 31 '16 at 13:09
  • I am unable to see any "edge effects" at all, perhaps because the images do not appear to distinguish between NA and actual values. Could you produce versions that visually differentiate NA from everything else? The focal mean looks like it's working correctly. Note, furthermore, that you only care about the results of the smoothing at the NA cells themselves: you will then use these imputed values to patch the original NA values, rather than replacing the original image with its focal mean. – whuber Apr 1 '16 at 19:27
5

There is a couple of things wrong with this:

mat <- matrix(1/25,ncol=5, nrow=5)
raster_res <- focal(raster1, mat, FUN="mean", na.rm=T) 

The argument to supply the function is called fun, not FUN. You give each value a weight of 1/25, and then you want to use "mean". However, you should use "sum" in that case! (which is the default, and because you were using FUN instead of fun, that is what happened anyway.) However, since you have missing values and you are using na.rm=TRUE, you really need to give each value a weight of 1 and then use mean.

Also, as you want to fill in missing values, not change existing values, I would use NAonly = TRUE, together with pad=TRUE (to pad virtual rows and columns with NAs outside of the raster).

Here is an example:

library(raster)

# example data
logo <- raster(system.file("external/rlogo.grd", package="raster")) 
set.seed(0)
i <- sample(ncell(logo), 200)
logo[i] <- NA
plot(logo)

m <- matrix(1, ncol=5, nrow=5)
r <- focal(logo, m, fun="mean", na.rm=TRUE, NAonly=TRUE, pad=TRUE) 
plot(r)

To see for which cells the values were estimated

plot(as(reclassify(is.na(logo), cbind(0, NA)), 'SpatialPolygons'), add=TRUE)

Also, if you want to write the raster to a file, instead of using writeRaster you should do that in one step:

r <- focal(logo, m, fun="mean", na.rm=TRUE, NAonly=TRUE, pad=TRUE, filename="smooth.tif", overwrite=TRUE)
  • Thank you very much for your response! So it was pretty much stupid rookie mistakes on my part... Now it works nicely and produces smooth images. The data values are also really close to what QGIS-r.neighbors produces. – weis_fn Apr 4 '16 at 6:01

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.