1

I have a training dataset containing around 350 polygons and their land cover classes which I use to train a random forest classifier to classify my raster image. I would like to use the values of all pixels covered by the training polygons in order to have a big enough training dataset. I have tried extracting pixel values to my raster but the classification result turned out to look very wrong. I have attached the codes below. Any suggestions on how I could approach this problem?

    training<- shapefile("2022_try_train.shp")
    training$class <- as.factor(training$class)
    raster <- raster("22_notfull_clip.tif")
    # Extract pixel values from training polygons
    training_pixels <- extract(raster, training, df=TRUE)
    training_pixels$ID <- as.factor(training_pixels$ID)
    training_pixels <- as.data.frame(training_pixels)
    # Train random forest model with response variable "landcover 
    class"
    rf_model <- randomForest(ID ~ ., data=training_pixels, 
    type="classification", ntree=500, mtry=3)

Update:I have modified my codes based on aldo_tapia's answer which is super clear and helpful. But I just ran into a new issue, all data is lost from my test data as I convert it to a data frame.

    training = vect('~/desktop/qgis/2022_try_train.shp')
    testing = vect('~/desktop/qgis/2022_try_test.shp')
    r = rast('~/desktop/qgis/22_notfull_clip.tif')


    training$ID = 1:nrow(training)
    r[r<0] = 0

    training_pixels = extract(r, training, df=TRUE)
    training_pixels %<>%
    left_join(data.frame(training[,c('ID','class')]), by = 'ID') %>% 
    select(-ID) %>% mutate(class = factor(class))

    training_pixels %>% group_by(class) %>% 
    summarise(n = n()) %>% ggplot(aes(y=n, x=class, fill=class)) +
    geom_bar(stat = "identity")

    training_pixels %>%
    group_by(class) %>% 
    summarise(total = n()) %>%
    ungroup() %>% 
    summarise(sample_size = min(total)) %>% 
    unlist() -> sample_size

    training_pixels %>%
    group_by(class) %>% 
    sample_n(sample_size) %>% 
    ungroup() -> samples

    samples$ID = 1:nrow(samples)

    train = list()
    test = list()
    classes = unique(samples$class)
    for (i in seq_along(classes)) {
    train[[i]] <- training[training$class == classes[i], ]
    test[[i]] <- testing[testing$class == classes[i], ]
    }

    train <- lapply(train, function(x) {
    x[, -which(names(x) == "ID")]
    })
    train <- do.call(rbind, train)  

    test <- lapply(test, function(x) {
    x[, -which(names(x) == "ID")]
    })
    test <- do.call(rbind, test)

    train <- as.data.frame(train)
    train %>% group_by(class) %>% 
    summarise(n = n())

    test <- as.data.frame(test)
    test %>% group_by(class) %>% 
    summarise(n = n())

    train$class <- as.factor(train$class)

    rf_model <- randomForest(class ~ ., data=train, 
    type="classification", ntree=500, mtry=3)

enter image description here

3
  • What have you tried and what problems were they? If you don't tell us, we might just hit the same problems.
    – Spacedman
    Commented Apr 12, 2023 at 11:13
  • Actually, the process is to extract values from the raster file and use it in the RF model. Isn't problematic at all, so please share you specific problem for helping you
    – aldo_tapia
    Commented Apr 12, 2023 at 13:20
  • Ok, I see the code but I don't see EDA. Maybe you have class imbalance, or you need to scale data. A machine-learning classification needs many pre-process before to train the model, maybe that's why you're failing here (for RF class imbalance is a huge problem, try bootstrap or reduce samples to the least frequent class, there are several techniques to deal with this issue)
    – aldo_tapia
    Commented Apr 12, 2023 at 15:08

2 Answers 2

0

Why are you using the "ID" column from extract? This explicitly relates to the individual polygon(s) row indices so, the assumptions would be that you only have a single AOI for each landcover class. Even in a nonparametric model this would equate to virtually no power and you cannot count on having enough variability within a given AOI to act as, what would equate to, independent training observations. However, with small sample sizes, I understand the motivation here.

You just have to ask, do I have enough power to discriminate the classes? With no insight on your data the first thing that comes to mind would be the amount of variability within each polygon, particularly given the possibility that you only have one polygon per landcover type. It is critical that you perform some exploratory data analysis (EDA) to see if there is any power. Ideally, you want low variability for an AOI so that the mean adequately represents the spectral characteristics of the training sample. You are banking on the opposite but, I do not see how high variability in an AOI would ever be a positive thing. It is very likley that you need more training samples and that is why your model is performing poorly.

Here is a quick worked example that results in the summary statistics for each band intersecting each polygon. I have moved the workflow to terra as the raster library is being depreciated in Oct.

First, create some example data that emulates polygon AOI regions and a multiband raster acting as spectral data. I am adding a "class" column to the polygons as well as an "ID" column representing row indices.

library(terra)
library(sf)

p <- st_read(system.file("ex/lux.shp", package="terra"))
  p$class <- sample(c("grass","forest","built"), nrow(p), replace=TRUE)
  p$ID <- 1:nrow(p)
    plot(p["class"])

r <- rast(ext(p), resolution=0.001, nlyrs=7, names=paste0("B",1:7))
  values(r) <- runif(ncell(r))

Here we extract the raster bands and relate the class column to the extracted data.frame.

v <- extract(r,vect(p))
  v <- data.frame(class=merge(st_drop_geometry(p), v, by="ID")$class, v)
    head(v)

Now, we can generate some summary statistics for each polygon. This is how you can evaluate the variability in the AOI pixels as training samples. Look at things like the deviation of the mean from the min and max and the variance. This will at least given some degree of insight to the expected behavior of your data. If too little variability you do not have power to predict to the entire scene and if too much there may be a poor signal to noise ratio.

by(v[,3:ncol(v)], v$class, FUN=summary)
by(v[,3:ncol(v)], v$class, FUN=var)




 
0

As Mr. Evans says, you can't use ID as your RF label. You need to add the labels from your polygons. Also, apply a EDA as first step.

Using a simple example with few classes (prepared really fast so isn't the best output)

library(randomForest)
library(terra)
library(dplyr)
library(ggplot2)

training = vect('~/path/samples.shp')
r = rast('~/path/spline.tif')
training$ID = 1:nrow(training)

r[r<0] = 0

plotRGB(r, r=1, g=6, b=12, scale = 300)
plot(training, 'class', add = TRUE)

enter image description here

training_pixels = extract(r, training, df=TRUE)

training_pixels %<>%
  left_join(data.frame(training[,c('ID','class')]), by = 'ID') %>% 
  select(-ID) %>% mutate(class = factor(class))

training_pixels %>% group_by(class) %>% 
  summarise(n = n()) %>% ggplot(aes(y=n, x=class, fill=class)) +
  geom_bar(stat = "identity")

enter image description here

Since there are more observations for class 5, the RF classifier will be biased to class 5.

You can either bootstrap the data or limit to the least frequent class (here I present option 2, which is easier to apply)

training_pixels %>%
  group_by(class) %>% 
  summarise(total = n()) %>%
  ungroup() %>% 
  summarise(sample_size = min(total)) %>% 
  unlist() -> sample_size

training_pixels %>%
  group_by(class) %>% 
  sample_n(sample_size) %>% 
  ungroup() -> samples

Here you get a balanced sample dataset for train/test the model. Then you can split data for training and testing:

samples$ID = 1:nrow(samples)

train = list()
test = list()

classes = unique(samples$class)

# 70% training, 30% testing
for(i in seq_along(classes)){
  temp = samples[samples$class == classes[i],]
  train[[i]] = sample_frac(temp, 0.7)
  test[[i]] = temp[!(temp$ID %in% train[[i]]$ID),]
}

train = do.call(rbind.data.frame, train) %>% select(-ID)
test = do.call(rbind.data.frame, test) %>% select(-ID)

train %>% group_by(class) %>% 
  summarise(n = n())
# # A tibble: 4 × 2
#   class     n
#   <fct> <int>
# 1 1       181
# 2 2       181
# 3 3       181
# 4 4       181
# 5 5       181

test %>% group_by(class) %>% 
  summarise(n = n())
# # A tibble: 4 × 2
#   class     n
#   <fct> <int>
# 1 1        77
# 2 2        77
# 3 3        77
# 4 4        77
# 5 5        77

Finally, train the model and compute accuracy:

rf_model <- randomForest(class ~ ., data=train, type="classification", ntree=500, mtry=3)
rf_model
# 
# Call:
#   randomForest(formula = class ~ ., data = train, type = "classification",      ntree = 500, mtry = 3) 
#                 Type of random forest: classification
#                       Number of trees: 500
# No. of variables tried at each split: 3
# 
#           OOB estimate of  error rate: 0%
# Confusion matrix:
#     1   2   3   4   5 class.error
# 1 181   0   0   0   0           0
# 2   0 181   0   0   0           0
# 3   0   0 181   0   0           0
# 4   0   0   0 181   0           0
# 5   0   0   0   0 181           0

modelled = predict(rf_model, test)
caret::confusionMatrix(test$class, modelled)
# Confusion Matrix and Statistics
# 
#             Reference
# Prediction  1  2  3  4  5
#          1 77  0  0  0  0
#          2  0 77  0  0  0
#          3  0  0 77  0  0
#          4  0  0  0 77  0
#          5  0  0  0  0 77
#
# more printed stuff...

As final step, predict using raster data:

class_values = predict(rf_model, r)
rclass = setValues(r[[1]], class_values)
plot(rclass)

enter image description here

9
  • Hi aldo_tapia, your answer is super helpful! One thing that wasn't clear to me was that when I convert my train/test dataset to dataframe somehow all data contained is lost. Would you mind clarifying that for me a bit? Thanks a lot again! Commented Apr 20, 2023 at 16:39
  • Can you edit the question and add more details of this issue? It's weird to lose all data
    – aldo_tapia
    Commented Apr 20, 2023 at 19:20
  • Hi I have edited the question and added the codes I am using. Hopefully it is clear enough and it would it great if you could have a look into the issue. Please let me know if there are any details needed. Thank you in advance for your help! Commented Apr 21, 2023 at 10:17
  • I see you have a training and testing vector. The code I shared involves splitting train and test data after balancing classes. I suggest you to merge both vectors. The most important step is the EDA, are your samples balanced?
    – aldo_tapia
    Commented Apr 21, 2023 at 13:51
  • In model specification, I strongly advocate against data withhold. This can be done efficiently and in a more informative way by bootstrapping a fit model. In a withhold you have no idea on the distributional characteristics that you are removing from the data. Also, in class balance, I believe that it is better dealing with it on the model side and not down or up sampling data (eg., SMOTE). Within random forests you can weight the bootstrap sample across classes or run multiple balanced ensembles with fixed parameter space, converging on sample covariance, then combine the results. Commented Apr 21, 2023 at 15:07

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.