Classifying raster that takes neighbours into account?

Question

I am working on classifying forest types using Sentinel-2 satellite imagery (10-m resolution; band-2, band-3, band-4 & band-8). I have ground-truthed field data as a training set. I want to set a parameter that constrains classification so that nearby pixels must be classified into the same class or that there must be a minimum of 10 neighbouring pixels within the same class. Ideally, I want the classification to take into account the heterogeneity of an area and not just classify each pixel based on the average in the training set. I would prefer a solution in R but I could also use QGIS or ArcGIS.

Details:

I have been using this helpful tutorial for the classification in R: http://ceholden.github.io/open-geo-tutorial/R/chapter_5_classification.html.

While I have successfully run the functions and algorithms to produce a classification using the randomForest package and the raster package, I am not satisfied with it because the result is lots of isolated pixels of different classes. (For example, five of the six classes are adjacent neighbours with no dominant class to use in post-processing to cluster together). Visually, I noticed that some forest types are distinguished by being heterogeneous (pixels of very different colours side-by-side) and other forest types are homogeneous (pixels of similar colours side-by-side). I want the classification algorithm to pick it on the pattern within the training set (polygons of circles 60m in diameter) and/or to take into account neighbours in the classification so that nearby pixels have the same value. I also have polygons of delineated forest types from photo-interpretation, which I could use to constrain to have only single value per polygon (the values of the photo interpretation are not accurate but the limits are good enough).

I've taken screenshots of the RGB true-colour image of three (out of six) of forest types and their respective classifications. I've chosen two for which it works well and I could use post-classification majority filter to smooth groups. And one for which there is a high heterogeneity and the classification algorithm does not work.

What function can I use in R to classify raster that takes neighbours into account?

Answer 1

Thank to all those who made comments on my question. Based on those comments, I was able to write code that took neighbours into account in my classification, which was successful. I am sharing my procedure here for others.

I created additional raster layers to use in the classification that used the values of neighbours with the focal function based on Matifou's suggestion. I created one raster where the value of each cell was the weighted mean of the neighbouring cell values (excluding itself) to include neighbouring values in calculations. I created another raster where the value of each cell was the standard deviation of a 3x3 window around the cell to capture the heterogeneity of a location.

#run loop to create RasterStack with weighted mean of 
#neighbour values for each pixel and each band
#img_2017_2dates_4band is my raster file with several dates and bands for the same location
library(raster)
m=matrix(c(1, 2, 1, 2, 0, 2, 1, 2, 1), ncol=3, nrow=3)
rs_focal=focal(img_2017_2dates_4band[[1]], fun=mean, w=m)
for (i in 2:nlayers(img_2017_2dates_4band))
{
rl=focal(img_2017_2dates_4band[[i]], fun=mean, w=m)
rs_focal=stack(rs_focal, rl)
}
names(rs_focal)=paste(names(img_2017_2dates_4band), "f")
#end of loop

#run loop to create RasterStack with standard deviation within band for
#set of pixels for each band
m=matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1), ncol=3, nrow=3)
rs_focal_sd=focal(img_2017_2dates_4band[[1]], fun=sd, w=m)
for (i in 2:nlayers(img_2017_2dates_4band))
{
rl=focal(img_2017_2dates_4band[[i]], fun=sd, w=m)
rs_focal_sd=stack(rs_focal_sd, rl)
}
names(rs_focal_sd)=paste(names(img_2017_2dates_4band), "f")
#end of loop

Instead of just using the raw raster values as predictors, I added the neighbour-value rasters to the rasterStack into one big raster.

img_focal=stack(img, rs_focal_n, rs_focal_sd)
roi_data=extract(x=img_layers, y=trainData, df=TRUE) #where trainData is  training polygon for classification
roi_data$class <- as.factor(trainData$Class_num[roi_data$ID])#assign class to each row as a column
rf <- randomForest(class ~ ., data=roi_data[,c(2:ncol(roi_data)], importance=TRUE)