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.


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

forest type eg

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?

  • 1
    As in other software, to delete salt and pepper effect you need to make a pixel-by-pixel classification and next apply a majority filter or aggregation. Another option could be to use OBIA approach
    – aldo_tapia
    Commented Sep 27, 2017 at 21:36
  • you seem to be describing a OBIA process, as aldo_tapia suggests: youtube.com/watch?v=aBw83S5gW2w
    – Elio Diaz
    Commented Sep 27, 2017 at 22:32
  • I don't think object-based will work because it is so heteregeneous that there aren't any clear distinctions based on individual pixel values. I'll try but I'm still hoping for something that aggregates before classification. Commented Sep 28, 2017 at 1:22
  • Taking account your update. Did you use random forest? How is your out of bag error? Are samples different enough? A good selection of samples is the most important part to obtain a good classification
    – aldo_tapia
    Commented Sep 28, 2017 at 13:49
  • It sounds like you want to establish a minimal mapping unit and not apply a conditional classification, which in many cases would be very problematic. Besides the above suggestions, you could look at a sieve filter (details here: gis.stackexchange.com/questions/91609/…). This is available in ArcGIS (via script), GDAL, ERDAS, QGIS and GRASS. Commented Sep 28, 2017 at 22:01

1 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
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)

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.