Introducción

El objetivo de este ejercicio es realizar la clasificación supervidada de un mosaico de imágenes multiespectrales de alta resolución capturadas con un vehículo aéreo no tripulado (UAV) sobre un cultivo experimental de papa ubicado en el municipio de Subachoque, Cundinamarca.

Datos

Se utilizará un mosaico de imágenes áereas multiespectrales obtenidas con el sensor Micasense RedEdge, la información de las bandas del sensor se muestra a continuación:

rededge <- read.csv("./datos/rededge.csv", sep = ",")

knitr::kable(rededge)
X..Banda Nombre Centro.de.Banda..nm. Ancho.Banda.FWHM..nm.
Banda 1 Azul 475 20
Banda 2 Verde 560 20
Banda 3 Rojo 668 10
Banda 4 Infrarrojo cercano 840 40
Banda 5 Red edge 717 10 
Def_clases <- read.csv("./datos/Clases.csv", sep = ",")

knitr::kable(Def_clases)
ID Clase Descripción
1 1 Agua
2 2 Suelo-Via
3 3 Pasto-Alto
4 4 Suelo-Cultivo
5 5 Pasto-Bajo
6 6 Papa
7 7 Papa-Tizón

A continuación se carga el mosaico multiespectral:

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.3-4, (SVN revision 766)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/jorod/OneDrive/Documentos/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/jorod/OneDrive/Documentos/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(tmap)
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
# Load the multi-band GeoTIFF file with brick function
rst <- brick("./datos/M_20180526_SPEC_3116_clip_10cm.tif")
rst_crs <- crs(rst, asText = TRUE)

# Change band names      
names(rst) <- paste("b",1:5,sep="")

El mosaico se puede ver a continuación, por conveniencia removemos los pixeles con valor 0 del mosaico. En la zona suroccidental del mosaico se encuentra un experimento de papa en el cual se utilizó una variedad de papa con diferentes tratamientos. Las imágenes aéreas fueron capturadas en un momento en el que el cultivo experimental había sufrido graves daños por causa del tizón tardío.

rst[rst==0] <- NA
plotRGB(rst, r=4, g=3, b=2, scale=10000, stretch="lin", main="Composición NIR-R-G (b4,b3,b2) Micasense RedEdge")

A continuación se cargan los poligonos de entrenamiento, se carga un archivo shapefile, que luego será rasterizado:

entrenamiento <- readOGR("./entrenamiento/entrenamiento.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jorod\OneDrive\Documentos\PRA\entrenamiento\entrenamiento.shp", layer: "entrenamiento"
## with 56 features
## It has 3 fields
## Integer64 fields read as strings:  id clase
entreCol <- "clase"

plotRGB(rst, r=3, g=2, b=1, scale=10000, stretch="lin", main="Composición NIR-R-G (b4,b3,b2) Micasense RedEdge")

plot(entrenamiento, col="yellow", border="black", lwd=1,
     main="Poligonos de entrenamiento", add=TRUE)

Método

Primero se rasteriza la capa de entrenamiento, para esto se crea un raster del mismo número de filas y columnas, además de la misma extencion del mosaico multiespectral. También se remueven los valores con 0 en la capa de entrenamiento.

rows <- nrow(rst)
cols <- ncol(rst)

Temptrain <- raster(ncol=cols, nrow=rows)
extent(Temptrain) <- extent(entrenamiento)

rtrain <- rasterize(entrenamiento,Temptrain,"clase")
rtrain[rtrain==0] <- NA

# Convert the data to factor/discrete representation
rtrain <- ratify(rtrain)

# Se cambia el nombre de la capa
names(rtrain) <- "trainClass"

plot(rtrain)

Ahora podemos ver cuantos píxeles tenemos por cada una de las 8 clases de entrenamiento.

tab_pix <- table(values(rtrain))
print(tab_pix)
## 
##     1     2     3     4     5     6     7     8 
## 18143  7315 12763   529  2682 69675  3567   288

Aún cuando se iguala la extensión del archivo vectorial con la imagen multiespectral, puede que el raster de entrenamiento generado difiera ligeramente en extensión por lo cual se debe realizar nuevamente la igualación de la extensión de la capa raster de entrenamiento con la del mosaico multiespectral.

## Extensión de la imagen multiespectral
rst
## class       : RasterBrick 
## dimensions  : 2214, 2435, 5391090, 5  (nrow, ncol, ncell, nlayers)
## resolution  : 0.1000064, 0.1000084  (x, y)
## extent      : 990784.6, 991028.1, 1035822, 1036044  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       :  b1,  b2,  b3,  b4,  b5 
## min values  :   1,   1,   1,   9,   1 
## max values  : 124, 128, 127, 129, 137
##Extensión del raster con las clases de entrenamiento
rtrain
## class       : RasterLayer 
## dimensions  : 2214, 2435, 5391090  (nrow, ncol, ncell)
## resolution  : 0.08614146, 0.07168819  (x, y)
## extent      : 990786.9, 990996.6, 1035833, 1035992  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : trainClass 
## values      : 1, 8  (min, max)
## attributes  :
##        ID
##  from:  1
##  to  :  8

Ahora se crea un dataframe que contiene los datos de la imagen multiespectral y de la capa de entrenamiento.

extent(rtrain) <- extent(rst)

rst_df <- na.omit(values(stack(rtrain, rst)))


rst_df[,"trainClass"] <- as.factor(as.character(rst_df[,"trainClass"]))

A continuación se inicializan las variables necesarias para el almacenamiento de la información

# Number of holdout evaluation rounds
nEvalRounds <- 20

# Proportion of the data used for training the classifier
pTrain <- 0.5

n <- nrow(rst_df)
nClasses <- length(unique(rst_df[,"trainClass"]))
#nClasses


# Initialize objects

confMats <- array(NA, dim = c(nClasses,nClasses,nEvalRounds))

evalMatrix<-matrix(NA, nrow=nEvalRounds, ncol=3, 
                   dimnames=list(paste("R_",1:nEvalRounds,sep=""), 
                                 c("Accuracy","Kappa","PSS")))

pb <- txtProgressBar(1, nEvalRounds, style = 3)
library(SDMTools)
## 
## Attaching package: 'SDMTools'
## The following objects are masked from 'package:caret':
## 
##     sensitivity, specificity
## The following object is masked from 'package:raster':
## 
##     distance
library(randomForest )
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
  # Create the random index for row selection at each round
  sampIdx <- sample(1:n, size = round(n*pTrain))
  #rst_df[sampIdx, -1]
  #rst_df[sampIdx, "trainClass"] <- as.factor(rst_df[sampIdx, "trainClass"])
  #rst_df[-sampIdx,]

  # Calibrate the RF classifier
  rf <- randomForest(y = rst_df[sampIdx, "trainClass"], 
                     x = rst_df[sampIdx, -1], 
                     ntree = 200)

  # Predict the class to the test set
  testSetPred <- predict(rf, newdata = rst_df[-sampIdx,], type = "response")
  #plot(testSetPred)
  testSetPred <-as.factor(as.integer(as.character(testSetPred)))
  #testSetPred
  #as.integer(testSetPred)
  #testSetPred <- factor(testSetPred)
  #testSetPred <- table(testSetPred)
  nrowtestPred <- nrow(testSetPred)
  #nrowtestPred

  # Get the observed class vector
  testSetObs <- as.factor(rst_df[-sampIdx,"trainClass"])
  #testSetObs
  #testSetObs <- factor(testSetObs)
  #testSetObs <- table(testSetObs)
  nrowtestObs <- nrow(testSetObs)
  #nrowtestObs
  
  confMatrix <- confusionMatrix(testSetPred,testSetObs)
## Warning in levels(reference) != levels(data): longitud de objeto mayor no
## es múltiplo de la longitud de uno menor
## Warning in confusionMatrix.default(testSetPred, testSetObs): Levels are not
## in the same order for reference and data. Refactoring data to match.
  #confMatrix
  
  # Classify the whole image with raster::predict function
  rstPredClass <- predict(rst, model = rf, 
                             factors = levels(rst_df[,"trainClass"]))
  #rstPredClass

Resultados

Clasificación del mosaico multiespectral.

plot(as.integer(rstPredClass))

confMatrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     1     2     3     4     5     6     7
##          1  5575   123   133     1     1    36     3
##          2  1583  1596   812    24    13   318    33
##          3   749   581  1561    59    50   696   126
##          4   354   171   607    51   165  1640   276
##          5   102    38   226    70    51 12088   752
##          6     2     0     8    51     0  3033   581
##          7     0     0     0     0     0     0     0
## 
## Overall Statistics
##                                           
##                Accuracy : 0.3456          
##                  95% CI : (0.3406, 0.3506)
##     No Information Rate : 0.5187          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.2553          
##  Mcnemar's Test P-Value : <2e-16          
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3 Class: 4 Class: 5 Class: 6
## Sensitivity            0.6665  0.63611  0.46639 0.199219 0.182143  0.17029
## Specificity            0.9886  0.91257  0.92705 0.905730 0.610206  0.96116
## Pos Pred Value         0.9494  0.36447  0.40842 0.015625 0.003827  0.82531
## Neg Pred Value         0.9020  0.96953  0.94148 0.993403 0.989101  0.51807
## Prevalence             0.2436  0.07307  0.09747 0.007455 0.008154  0.51868
## Detection Rate         0.1624  0.04648  0.04546 0.001485 0.001485  0.08833
## Detection Prevalence   0.1710  0.12752  0.11130 0.095052 0.388101  0.10702
## Balanced Accuracy      0.8275  0.77434  0.69672 0.552474 0.396174  0.56572
##                      Class: 7
## Sensitivity           0.00000
## Specificity           1.00000
## Pos Pred Value            NaN
## Neg Pred Value        0.94843
## Prevalence            0.05157
## Detection Rate        0.00000
## Detection Prevalence  0.00000
## Balanced Accuracy     0.50000

Discusión

En la imagen clasificada se puede observar que la clase 7, correspondiente a papa con tizón tardío, no logró ser identificada, ya que es confundida con otras coberturas, como el suelo de la vía (clase 2) o el suelo del cultivo (clase 4), de hecho, la pérdida de área foliar en el momento de la captura de las imágenes era muy alta, por lo cual estos resultados en la clasificación eran esperados. La clase 1, correspondiente al agua fue bien identificada aunque tambien tuvo falsos positivos en zonas de pasto.

En general los resultados no son muy buenos, se tiene una precisión de 0.33 y un kappa de 0.24 es decir que el resultado de usar el clasificador empleado acá es apenas 24% mejor que utilizar un clasificador aleatorio que asigne los píxeles al azar.

Estos resultados tienen una causa principal, que es la calibración radiométrica de las imágenes aéreas obtenidas, en este caso el mosaico empleado proviene de imágenes aéreas sin calibración radiométrica, se esperaría que una buena calibración radiométrica pudiera generar resultados mucho mejores.