Clasificación de imágenes con SVM

Librerias aplicadas

library(caret)
library(raster)

Carga de las bandas de imágenes satelitales.

band1 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B1.TIF")
band2 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B2.TIF")
band3 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B3.TIF")
band4 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B4.TIF")
band5 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B5.TIF")
band6 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B6.TIF")
band7 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B7.TIF")
band10 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B10.TIF")
band11 <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/capturas-landsat/roma/LC81910312013208LGN00/LC81910312013208LGN00_B11.TIF")

roma_rst_lan <- brick(band1,band2,band3,band4,band5,band6,band7,band10,band11)
names(roma_rst_lan) <- c(paste("Banda",1:7,sep=""),"Banda10","Banda11")

Procesado de datos

#Carga de muestras de entrenamiento
rstTrainroma <- raster("/Users/elianaalizchuquillanquijulcapari/Documents/workspace-r/demos/DataTraining/rome_lcz_GT.tif")

# Remover ceros (fondo NA)
rstTrainroma[rstTrainroma==0] <- NA

# Convertir a tipo discreto
rstTrainroma <- ratify(rstTrainroma)

# Cambiamos el nombre
names(rstTrainroma) <- "layer"

#Muestras de valores
table(values(rstTrainroma))

   2    3    5    6    8   10   11   12   14   17 
1551  104 1495  480  435   51  284  555  984  500 

# Eliminación de NAs 
rstDFlanroma.df <- na.omit(values(stack(rstTrainroma, roma_rst_lan)))

cols <- c("#fde725","#b5de2b", "#6ece58","#35b779", "#1f9e89","#26828e", "#31688e", "#3e4989", "#482878", "#440154")
plot(rstTrainroma, main="Muestras de entrenamiento de Roma", col = cols)

Fase de Clasificación

#Si es la primera vez se ejecuta, sino ya esta guardado los datos del dataframe.
#write.csv(rstDFlanroma.df,"rstDFRoma.csv", row.names = FALSE)

#Leemos los datos
samplesLanroma = read.csv("rstDFRoma.csv")

# Se divide las muestras en 70 y 30
trainx = list(0)
evalx = list(0)
for (i in 1:17){ # se itera en todas las posibles categorías
  cls = samplesLanroma[samplesLanroma$layer == i,]
  smpl <- floor(0.70 * nrow(cls))
  tt <- sample(seq_len(nrow(cls)), size = smpl)
  trainx[[i]] <- cls[tt,]
  evalx[[i]] <- cls[-tt,]
}


# Se combina en dataframe referenciados s entrenamiento y evaluación
trn = do.call(rbind, trainx) 
eva = do.call(rbind, evalx)

Algoritmo SVM

# Definiendo el control del entrenamiento
fitControl <- trainControl(method = "repeatedcv", # repeated cross-validation of the training data
                   number = 10, # number of folds
                   repeats = 5) # view the training iterations


# Entrenamiento con svm
svm_roma_model<- caret::train(x = trn[,(2:ncol(trn))], y = as.factor(trn$layer),
                      method = "svmRadial", metric="Accuracy", 
                      trControl = fitControl, 
                      preProcess = c("center", "scale"),
                      tuneLength = 10)
# Entrenamiento modelo 2
svm_grid <- expand.grid(sigma=seq(from = 0.01, to = 0.10, by = 0.02),
                        C=seq(from = 2, to = 10, by = 2))

svm_roma_model2<- caret::train(x = trn[,(2:ncol(trn))], y = as.factor(trn$layer),
                      method = "svmRadialSigma", metric="Accuracy", 
                      trControl = fitControl,
                      preProcess = c("center", "scale"),
                      tuneGrid = svm_grid)

Exploración del modelo

#Información modelo svm
print(svm_roma_model)

Support Vector Machines with Radial Basis Function Kernel 

4502 samples
   9 predictor
  10 classes: '2', '3', '5', '6', '8', '10', '11', '12', '14', '17' 

Pre-processing: centered (9), scaled (9) 
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 4050, 4053, 4053, 4052, 4051, 4054, ... 
Resampling results across tuning parameters:

  C       Accuracy   Kappa    
    0.25  0.7602406  0.7111069
    0.50  0.7725452  0.7261232
    1.00  0.7821401  0.7377471
    2.00  0.7939981  0.7521567
    4.00  0.8023949  0.7623404
    8.00  0.8109229  0.7729743
   16.00  0.8186064  0.7826110
   32.00  0.8265611  0.7923832
   64.00  0.8266494  0.7925697
  128.00  0.8223409  0.7874989

Tuning parameter 'sigma' was held constant at a value of 0.369381
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were sigma = 0.369381 and C = 64.

print(svm_roma_model2) #modelo 2

Support Vector Machines with Radial Basis Function Kernel 

4502 samples
   9 predictor
  10 classes: '2', '3', '5', '6', '8', '10', '11', '12', '14', '17' 

Pre-processing: centered (9), scaled (9) 
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 4051, 4051, 4052, 4053, 4053, 4054, ... 
Resampling results across tuning parameters:

  sigma  C   Accuracy   Kappa    
  0.01    2  0.7425628  0.6887592
  0.01    4  0.7530479  0.7018675
  0.01    6  0.7573136  0.7071743
  0.01    8  0.7608221  0.7114416
  0.01   10  0.7635315  0.7148006
  0.03    2  0.7585560  0.7087919
  0.03    4  0.7675344  0.7198353
  0.03    6  0.7753058  0.7293252
  0.03    8  0.7801483  0.7352591
  0.03   10  0.7835245  0.7394246
  0.05    2  0.7667807  0.7189443
  0.05    4  0.7795289  0.7345235
  0.05    6  0.7860566  0.7425342
  0.05    8  0.7905899  0.7480190
  0.05   10  0.7942355  0.7524134
  0.07    2  0.7748659  0.7288664
  0.07    4  0.7860116  0.7424993
  0.07    6  0.7924144  0.7502067
  0.07    8  0.7958818  0.7544567
  0.07   10  0.7985892  0.7577439
  0.09    2  0.7794400  0.7344683
  0.09    4  0.7903707  0.7477745
  0.09    6  0.7957463  0.7542856
  0.09    8  0.7986329  0.7578170
  0.09   10  0.8013002  0.7610445

Accuracy was used to select the optimal model using the largest value.
The final values used for the model were sigma = 0.09 and C = 10.

#Gráfica del modelo svm
plot(svm_roma_model)

plot(svm_roma_model2) #modelo 2

Evaluación del modelo 1.

pred_svm_roma_model<- predict(svm_roma_model,newdata = eva[,-1])

cm <- confusionMatrix(data = pred_svm_roma_model, as.factor(eva$layer))
cm

Confusion Matrix and Statistics

          Reference
Prediction   2   3   5   6   8  10  11  12  14  17
        2  394  12  89   7   4   2   0   0   0   0
        3    4  19   2   4   0   0   0   0   0   0
        5   63   1 314  17  11   8   2  13   1   0
        6    3   0  17 108   0   1   1   6   4   0
        8    2   0  12   0 116   0   0   0   0   0
        10   0   0   0   0   0   3   0   0   0   0
        11   0   0   0   0   0   0  81   2   0   0
        12   0   0  10   2   0   0   2 140   6   0
        14   0   0   5   6   0   2   0   6 285   0
        17   0   0   0   0   0   0   0   0   0 150

Overall Statistics
                                          
               Accuracy : 0.8312          
                 95% CI : (0.8137, 0.8476)
    No Information Rate : 0.2406          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.7981          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: 2 Class: 3 Class: 5 Class: 6 Class: 8
Sensitivity            0.8455 0.593750   0.6993  0.75000  0.88550
Specificity            0.9225 0.994751   0.9220  0.98215  0.99225
Pos Pred Value         0.7756 0.655172   0.7302  0.77143  0.89231
Neg Pred Value         0.9496 0.993187   0.9104  0.97997  0.99170
Prevalence             0.2406 0.016520   0.2318  0.07434  0.06763
Detection Rate         0.2034 0.009809   0.1621  0.05576  0.05989
Detection Prevalence   0.2623 0.014972   0.2220  0.07228  0.06711
Balanced Accuracy      0.8840 0.794250   0.8107  0.86608  0.93887
                     Class: 10 Class: 11 Class: 12 Class: 14 Class: 17
Sensitivity           0.187500   0.94186   0.83832    0.9628   1.00000
Specificity           1.000000   0.99892   0.98870    0.9884   1.00000
Pos Pred Value        1.000000   0.97590   0.87500    0.9375   1.00000
Neg Pred Value        0.993278   0.99730   0.98481    0.9933   1.00000
Prevalence            0.008260   0.04440   0.08622    0.1528   0.07744
Detection Rate        0.001549   0.04182   0.07228    0.1471   0.07744
Detection Prevalence  0.001549   0.04285   0.08260    0.1569   0.07744
Balanced Accuracy     0.593750   0.97039   0.91351    0.9756   1.00000

Representación de Matriz de confusión

cm_d <- as.data.frame(cm$table)

cm_st <-data.frame(cm$overall)

cm_st$cm.overall <- round(cm_st$cm.overall,2)

cm_p <- as.data.frame(prop.table(cm$table))
cm_d$Perc <- round(cm_p$Freq*100,2)

# dibujando la matriz
ggplot(data = cm_d, aes(x = Prediction , y =  Reference, fill = Freq))+
  geom_tile(color = "white", lwd = 0.5) +
  geom_text(aes(label = paste("",Freq)), color = 'white', size = 3.5)+ 
  theme_minimal() +
  ggtitle("Matriz de Confusión")

Evaluación del modelo 2.

pred_svm_roma_model2<- predict(svm_roma_model2,newdata = eva[,-1])

cm <- confusionMatrix(data = pred_svm_roma_model2, as.factor(eva$layer))
cm

Confusion Matrix and Statistics

          Reference
Prediction   2   3   5   6   8  10  11  12  14  17
        2  385  26  86   8   1   2   0   0   0   0
        3    0   0   0   0   0   0   0   0   0   0
        5   78   3 319  24  16  11   0  18   6   0
        6    2   3  18  99   0   0   0   2   3   0
        8    1   0   8   0 114   0   0   0   0   0
        10   0   0   0   0   0   2   0   0   0   0
        11   0   0   0   0   0   0  80   3   0   0
        12   0   0  10  10   0   0   5 130  13   0
        14   0   0   8   3   0   1   1  14 274   0
        17   0   0   0   0   0   0   0   0   0 150

Overall Statistics
                                          
               Accuracy : 0.8018          
                 95% CI : (0.7833, 0.8193)
    No Information Rate : 0.2406          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.7617          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: 2 Class: 3 Class: 5 Class: 6 Class: 8
Sensitivity            0.8262  0.00000   0.7105  0.68750  0.87023
Specificity            0.9164  1.00000   0.8952  0.98438  0.99502
Pos Pred Value         0.7579      NaN   0.6716  0.77953  0.92683
Neg Pred Value         0.9433  0.98348   0.9111  0.97514  0.99063
Prevalence             0.2406  0.01652   0.2318  0.07434  0.06763
Detection Rate         0.1988  0.00000   0.1647  0.05111  0.05885
Detection Prevalence   0.2623  0.00000   0.2452  0.06557  0.06350
Balanced Accuracy      0.8713  0.50000   0.8028  0.83594  0.93262
                     Class: 10 Class: 11 Class: 12 Class: 14 Class: 17
Sensitivity           0.125000   0.93023   0.77844    0.9257   1.00000
Specificity           1.000000   0.99838   0.97853    0.9835   1.00000
Pos Pred Value        1.000000   0.96386   0.77381    0.9103   1.00000
Neg Pred Value        0.992765   0.99676   0.97908    0.9866   1.00000
Prevalence            0.008260   0.04440   0.08622    0.1528   0.07744
Detection Rate        0.001033   0.04130   0.06711    0.1415   0.07744
Detection Prevalence  0.001033   0.04285   0.08673    0.1554   0.07744
Balanced Accuracy     0.562500   0.96431   0.87849    0.9546   1.00000

Representación de Matriz de confusión

cm_d <- as.data.frame(cm$table)

cm_st <-data.frame(cm$overall)

cm_st$cm.overall <- round(cm_st$cm.overall,2)

cm_p <- as.data.frame(prop.table(cm$table))
cm_d$Perc <- round(cm_p$Freq*100,2)

# dibujando la matriz
ggplot(data = cm_d, aes(x = Prediction , y =  Reference, fill = Freq))+
  geom_tile(color = "white", lwd = 0.5) +
  geom_text(aes(label = paste("",Freq)), color = 'white', size = 3.5)+ 
  theme_minimal() +
  ggtitle("Matriz de Confusión")

Predicción y visualización del modelo 1.

classi_svm_model_roma <- raster::predict(roma_rst_lan, model=svm_roma_model)
plot(classi_svm_model_roma, main="Clasificación SVM Modelo 1",col = cols)

Predicción y visualización del modelo 2.

classi_svm_model_roma2 <- raster::predict(roma_rst_lan, model=svm_roma_model2)
plot(classi_svm_model_roma2, main="Clasificación SVM Modelo 2",col = cols)