Clasificación Supervisada de la cobertura del suelo en Zona Oriental (Colombia) sobre cultivos de Palma de Aceite usando el software R

Introducción

Contexto

En Colombia el cultivo de la palma de aceite (E. guineensis) representa una gran importancia en la economía nacional, teniendo en cuenta su adaptación y fácil manejo dadas las condiciones climáticas del país, motivo por el cual, Colombia es el principal cultivador y productor de palma de aceite en América y ocupa el cuarto lugar a nivel mundial (Fedepalma, 2018). La palma de aceite se considera de alto rendimiento por la cantidad de aceite que produce su fruto por hectárea y por la variedad de productos obtenidos no solo de su fruto sino de otras partes de la planta, de ella se extrae aceite comestible, aunque cuenta con diversos usos como productos alimenticios, medicinales, fabricación de fibras, entre otros. En Colombia, los departamentos productores se dividen en diferentes zonas, las cuales son: Zona Oriental, Norte, Centro y Suroccidental; en donde se tienen grandes extensiones de tierra, por lo que innovaciones tecnológicas que se logren implementar puede hacer que este cultivo se desarrolle de una forma óptima.

Ademas enfermedades como la ML son uno de los retos fitosanitarios de la palmicultura en la Zona Oriental de Colombia, debido a que es una de las enfermedades más limitantes para este cultivo en esta región y cuyo impacto económico impacto económico se estimó en más de US 24.203.507, durante el período comprendido entre 1994 y 2010 en el Bajo Upia. La ML provocó la destrucción de casi mil hectáreas en el Bajo Upía, en la región Oriental entre 1994 y 2013. Por eso es de gran importancia hacer seguimiento del cultivo, ya que se presentan diferentes factores abioticos y bioticos que pueden afectar su crecimiento, desarollo y asi mismo su produccion de aceite.

Objetivos

General

Realizar una clasificacion supervisada sobre un cultivo de palma aceite en Zona Oriental mediante imagenes multiespectrales adquiridas de sensores transportados por Aeronaves Remotamente Tripuladas (ART).

Especificos

• Analizar las estadisticas descriptivas de las imagenes adquiridas del sensor multiespetral sobre la zona de interes.

• Diferenciar la cobertura de palma de aceite del suelo y de la vegetacion secundaria para visualizar alteraciones en el cultivo de palma de aceite debido a factores bioticos y/o abioticos.

• Evaluar el metodo de clasificacion supervisada de Decision Trees (DT) sobre la zona de estudio.

Datos y Métodos

Las imagenes que se utilizaron en el presente proyecto fueron adquiridas por medio del sensor multiespetral RedEdge (MicaSense) adaptado a una Aeronave Remotamente Tripulada (ART), especificamente un multirotor, para sobrevolar un lote de Palma de Aceite en el departamento del Meta (Zona Oriental - Colombia, La principal actividad económica de la zona de estudio se basa en la ganadería y producción agrícola, encontrando grandes extensiones de pastos, palma de aceite, bosques, entre otros). Esto con el fin de obtener imagenes de alta resolucion y asi poder identificar o diagnosticar alteraciones dentro del cultivoLos por factores bioticos o abioticos. Los diferentes vuelos para la toma de imagenes se realizaron a una altura de 100 metros sobre la superficie para obtener un Ground Sampled Distance (GSD) de 6.8 cm/pixel aproximadamente (MicaSense).

Para realizar los diferentes procedimiento del flujo de trabajo se usaron las librerias dismo, glcm, sp, raster, rasterVis, rgdal, rgeos y rpart.

Indices Multiespectrales de Vegetacion

Para el calculo de diferentes indices multiespetrales de vegetacion, se tuvieron en cuenta indices reportados en la literatura (Index DataBase) e indices ya evaluadas en el cultivo con el fin de observar si estos indices nos indican o muestran diferencias entre las diferentes coberturas que se pueden encontrar en un lote de palma.

Textura

En cuanto a la determinacion de textura en las imagenes, se determinaron las diferentes metricas de textura y se realizaron pruebas con diferentes kernel teniendo en cuenta la resolucion de la imagenes y el tamaño del objeto de estudio (99, 1717, 3535, 7171, y 93*93), con el fin de observar cual de estas ventanas mostraban mejores resultados en la discriminacion de coberturas.

Clasificacion Supervisada

Para realizar la clasificacion supervisada se hizo uso del metodo de Decision Trees (DT) ya que nos permite ver el diseño del arbol de decision que se utilizo para realizar la clasificaion. Para esta se determinaron 6 clases.

Resultados y Discusión

Objetivo I

Las imagenes multiespectrales adquiridas por el sensor RedEdge tienen 4188 filas y 5040 columnas dando un numero total de 211070520 pixeles en la escena adquirida en el mes de Mayo del 2019, presente una resolucion de 6.4 cm/pixel aproximadamente y un pero aproximado de 0.84 gigabytes (GB).

## [1] 4188 5040    5

## [1] 21107520

## [1] 0.0636323 0.0636323

## 844314968 bytes

Para evitar errores de los bordes de los ortomosaicos en las diferentes bandas del sensor se realizo un recorte teniendo en cuenta abarcar la mayor cantidad de palmas de aceite. Ademas se evitaron coberturas como arbustos, matorrales o bosques, ya que el interes del trabajo es diferenciar las palmas de vegetacion secundaria como gramineas y leguminosas como frijol.

# Extent of interest
E <- shapefile("E:/Master_Project_Data/Further_Info/Zona")
E

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 1
## names       : Id 
## value       :  0

# Raster plot
par(mfrow = c(1,1))
plotRGB(PLC_MayI, r=3,g=2,b=1, axes=TRUE, stretch="lin",
        main="MicaSense True Color Composite")
plot(E, add=TRUE)

# Crop MicaSense by the extent
E_PLC_MayI <- crop(PLC_MayI, E); E_PLC_MayI

## class      : RasterBrick 
## dimensions : 2492, 3534, 8806728, 5  (nrow, ncol, ncell, nlayers)
## resolution : 0.0636323, 0.0636323  (x, y)
## extent     : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      :        Blue,       Green,         Red,     RedEdge,         NIR 
## min values : 0.003112840, 0.010345617, 0.003387503, 0.022583352, 0.115815976 
## max values :   0.1207294,   0.2369116,   0.2265660,   0.5106584,   0.7437858

#writeRaster(E_PLC_MayI, "Extent_PLC_MayI.tif", format="GTiff", overwrite=TRUE)
plot(E_PLC_MayI, col=gray(0:100 / 100))

Segun la estadistica descriptiva la banda del infrarrojo cercano (NIR) es la banda que presenta mayor desviacion estandar y varianza por lo que se puede inferir que es una de las bandas que nos proporciona mayor informacion, logrando diferenciar o discriminar las coberturas que estan presentes en el area seleccionada.

# Correlation between bands
pairs(E_PLC_MayI, main="ALL BANDS vs ALL BANDS")

# Calculate layer statistics
summary(E_PLC_MayI)

##               Blue      Green         Red    RedEdge       NIR
## Min.    0.00311284 0.01034562 0.003387503 0.02258335 0.1158160
## 1st Qu. 0.02166781 0.06399634 0.026489662 0.15356680 0.3366445
## Median  0.02752728 0.07858396 0.033691920 0.18335241 0.3971313
## 3rd Qu. 0.03430228 0.09674220 0.043305104 0.22043183 0.4657664
## Max.    0.12072938 0.23691157 0.226565957 0.51065843 0.7437858
## NA's    0.00000000 0.00000000 0.000000000 0.00000000 0.0000000

Ref_sd <- cellStats(E_PLC_MayI, stat='sd', na.rm=TRUE); Ref_sd

##       Blue      Green        Red    RedEdge        NIR 
## 0.01098986 0.02553556 0.02294592 0.05354388 0.08849730

Ref_var <- cellStats(E_PLC_MayI, stat='var', na.rm=TRUE); Ref_var

##         Blue        Green          Red      RedEdge          NIR 
## 0.0001207771 0.0006520647 0.0005265152 0.0028669471 0.0078317728

Ref_skew <- cellStats(E_PLC_MayI, stat='skew', na.rm=TRUE); Ref_skew

##      Blue     Green       Red   RedEdge       NIR 
## 1.6578077 0.8769409 2.6616931 0.8211554 0.1279812

# Correlation matrix
Cor <- layerStats(E_PLC_MayI, 'pearson', na.rm = TRUE); Cor

## $`pearson correlation coefficient`
##                Blue     Green        Red   RedEdge         NIR
## Blue     1.00000000 0.6480373  0.8666857 0.4448153 -0.02471353
## Green    0.64803729 1.0000000  0.4440459 0.9050439  0.42952828
## Red      0.86668570 0.4440459  1.0000000 0.1737107 -0.36772783
## RedEdge  0.44481525 0.9050439  0.1737107 1.0000000  0.67286580
## NIR     -0.02471353 0.4295283 -0.3677278 0.6728658  1.00000000
## 
## $mean
##       Blue      Green        Red    RedEdge        NIR 
## 0.02916007 0.08237979 0.03968080 0.19136677 0.40137466

# Covariance matrix
Cov <- layerStats(E_PLC_MayI, 'cov', na.rm = TRUE); Cov

## $covariance
##                  Blue        Green           Red      RedEdge
## Blue     1.207771e-04 0.0001818602  0.0002185543 0.0002617471
## Green    1.818602e-04 0.0006520647  0.0002601829 0.0012374419
## Red      2.185543e-04 0.0002601829  0.0005265152 0.0002134233
## RedEdge  2.617471e-04 0.0012374419  0.0002134233 0.0028669471
## NIR     -2.403572e-05 0.0009706600 -0.0007467273 0.0031883672
##                   NIR
## Blue    -2.403572e-05
## Green    9.706600e-04
## Red     -7.467273e-04
## RedEdge  3.188367e-03
## NIR      7.831773e-03
## 
## $mean
##       Blue      Green        Red    RedEdge        NIR 
## 0.02916007 0.08237979 0.03968080 0.19136677 0.40137466

# OIF
Matriz_cor = Cor[[1]]
nrow(Matriz_cor) # numero de bandas = 5

## [1] 5

choose(5, 3)     # numero de combinaciones de 5 bandas, de 3 en 3 = 10

## [1] 10

id_comb = combn(1:5, 3) # matriz de las 10 posibles combinaciones de indices (i,j,k de la formula)
OIF = 0 #inicia variable en 0
for(n in 1:10){
        i = id_comb[1,n]
        j = id_comb[2,n]
        k = id_comb[3,n]
        OIF[n]=(Ref_sd[[i]]+Ref_sd[[j]]+Ref_sd[[k]])/(Matriz_cor[i,j]+Matriz_cor[i,k]+Matriz_cor[j,k])
}
ID_OIF <- rbind(id_comb, OIF)
ID_OIF <- t(ID_OIF); ID_OIF

##                    OIF
##  [1,] 1 2 3 0.03036159
##  [2,] 1 2 4 0.04508207
##  [3,] 1 2 5 0.11874672
##  [4,] 1 3 4 0.05890047
##  [5,] 1 3 5 0.25816458
##  [6,] 1 4 5 0.14001427
##  [7,] 2 3 4 0.06699850
##  [8,] 2 3 5 0.27079126
##  [9,] 2 4 5 0.08347792
## [10,] 3 4 5 0.34454959

data_ID_OIF <- data.frame(ID_OIF); data_ID_OIF

##    V1 V2 V3        OIF
## 1   1  2  3 0.03036159
## 2   1  2  4 0.04508207
## 3   1  2  5 0.11874672
## 4   1  3  4 0.05890047
## 5   1  3  5 0.25816458
## 6   1  4  5 0.14001427
## 7   2  3  4 0.06699850
## 8   2  3  5 0.27079126
## 9   2  4  5 0.08347792
## 10  3  4  5 0.34454959

Order_OIF <- sort(data_ID_OIF$OIF, decreasing = T); Order_OIF

##  [1] 0.34454959 0.27079126 0.25816458 0.14001427 0.11874672 0.08347792
##  [7] 0.06699850 0.05890047 0.04508207 0.03036159

En cuanto al Factor de Indice Optimo (OIF, por sus siglas en inglés) se puede observar que los mayores valores se presentaron en las combinaciones 543, 542, 531, 541 y 521; en las cuales se observa la presencia de la banda NIR y RedEdge; lo cual indica una baja correlación y valores de desviación estándar altos. Un valor alto en el OIF también nos da información sobre las bandas que presentan mayor contraste entre si, e igualmente nos indica que la información que se encuentra en cada una de las bandas no se repita; por lo que las bandas 5(NIR) y 4(RedEdge) generan o proporcionan datos de gran importancia sobre la superficie en la cual fue capturadas las imagenes.

par(mfrow=c(1,2))
plotRGB(PLC_MayI, r=3,g=2,b=1, axes=TRUE, stretch="lin",
        main="MicaSense True Color Composite")
plotRGB(PLC_MayI, r=5,g=4,b=3, axes=TRUE, stretch="lin",
        main="MicaSense True Color Composite")

Para el calculo de los perfiles espetrales a partir de los ortomosaicos se determinaron 6 clases (1.Palma Aparentemente Sana 2.Palma Enferma 3.Palma Erradicada 4.Palma Inyetada 5.Suelo y 6.Vegetacion Secundaria), para observar como era el comportamiento de estos y determinar que banda del sensor podria dar valores de reflectancia diferentes para las coberturas determinadas.

par(mfrow=c(1,1))
Tra_Pol<-shapefile('E:/Master_Project_Data/Further_Info/Training_Polygons.shp')
Tra_Pol<-crop(Tra_Pol, E); Tra_Pol

## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 1
## names       : Land_Cover 
## min values  :  Palm_ApSa 
## max values  :    Veg_Sec

plot(Tra_Pol, axes=TRUE, stretch="lin",
     main="Training Polygons")
plot(E, add=TRUE)

# Generate 1000 point samples from the polygons
Pt_Samp <- spsample(Tra_Pol, 1000, type='regular')
# Add the land cover class to the points
Pt_Samp$Land_Cover <- over(Pt_Samp, Tra_Pol)$Land_Cover
# extract values with points
DF <- extract(E_PLC_MayI, Pt_Samp)
# To see some of the reflectance values
head(DF)

##            Blue      Green        Red   RedEdge       NIR
## [1,] 0.02490272 0.07449454 0.04303044 0.1864347 0.2723430
## [2,] 0.02172885 0.06570535 0.03039597 0.1245747 0.2208591
## [3,] 0.02920577 0.09411765 0.03833066 0.2128023 0.3227893
## [4,] 0.03906310 0.07272450 0.08010986 0.1649805 0.2359045
## [5,] 0.03048753 0.06860456 0.05197223 0.1558251 0.2298314
## [6,] 0.02679484 0.08084230 0.03930724 0.1871366 0.3300832

En el siguiente grafico se puede observar que las bandas del RedEdge y del Infrarrojo Cercano (NIR) son las que presentan mayor variabilidad en la refletancia espetral como se pudo inferir anteriormente segun la estadisticas descriptivas de zona de interes en cada una de las bandas.

Mean_Spec <- aggregate(DF, list(Pt_Samp$Land_Cover), mean)
Median_Spec <- aggregate(DF, list(Pt_Samp$Land_Cover), median)
# instead of the first column, we use row names
rownames(Mean_Spec) <- Mean_Spec[,1]
Mean_Spec <- Mean_Spec[,-1]; Mean_Spec # Mean

##                   Blue      Green        Red   RedEdge       NIR
## Palm_ApSa   0.02789555 0.06970967 0.03225004 0.1702003 0.4199402
## Palm_Enf    0.02749777 0.08174843 0.03363266 0.2104331 0.4039671
## Palm_Erra   0.03762657 0.06626339 0.07186781 0.1174422 0.2070431
## Palm_Inject 0.03273470 0.07108620 0.08485728 0.1663598 0.3004903
## Suelo       0.03856090 0.08977948 0.06659217 0.1887751 0.3130314
## Veg_Sec     0.02622310 0.09214685 0.03243056 0.2145819 0.4466799

rownames(Median_Spec) <- Median_Spec[,1]
Median_Spec <- Median_Spec[,-1]; Median_Spec # Median

##                   Blue      Green        Red   RedEdge       NIR
## Palm_ApSa   0.02720684 0.06692607 0.03187610 0.1656367 0.4177920
## Palm_Enf    0.02661173 0.07840085 0.03292897 0.2039979 0.3970397
## Palm_Erra   0.03860533 0.06243992 0.07104601 0.1078508 0.1901732
## Palm_Inject 0.03137255 0.06466773 0.08078126 0.1559472 0.2870527
## Suelo       0.03521782 0.08749523 0.05587854 0.1837797 0.3102464
## Veg_Sec     0.02453651 0.07928588 0.02969406 0.1895171 0.4377508

# Create a vector of color for the land cover classes for use in plotting
MyColor <- c('brown', 'green', 'gray', 'yellow', 'black', 'blue')

par(mfrow=c(1,2))
# Transform Mean or Median from a data.frame to a matrix
Mean_Spec <- as.matrix(Mean_Spec)
# First create an empty plot
plot(0, ylim=c(0,0.5), xlim=c(1,5), type='n',
     xlab="Bands", ylab="Reflectance")
# Add the different classes
for (i in 1:nrow(Mean_Spec)){
        lines(Mean_Spec[i,], type="l", lwd=2, lty=1, col=MyColor[i])
}
# Title
title(main="Spectral Profile (Mean)", font.main=2)
# Legend
legend("topleft", rownames(Mean_Spec), inset=0.01,
       cex=0.9, col=MyColor, lty=1, lwd=3, bty="n")
Median_Spec <- as.matrix(Median_Spec)
# First create an empty plot
plot(0, ylim=c(0,0.5), xlim=c(1,5), type='n',
     xlab="Bands", ylab="Reflectance")
# Add the different classes
for (i in 1:nrow(Median_Spec)){
        lines(Median_Spec[i,], type="l", lwd=2, lty=1, col=MyColor[i])
}
# Title
title(main="Spectral Profile (Median)", font.main=2)
# Legend
legend("topleft", rownames(Median_Spec), inset=0.01,
       cex=0.9, col=MyColor, lty=1, lwd=3, bty="n")

Objetivo II

Para diferenciar la cobertura de palma de aceite del suelo y de la vegetacion secundaria se realizo el calculo de diferentes indices multiespetrales de vegetacion, pero en muchos de ellos no se observaron grandes diferencias como el mas usado NDVI, por lo cual no se tuvieron en cuenta ya no brindaban mucha informacion al momento de discriminar coberturas.

ARI_II <- overlay(E_PLC_MayI$Green, E_PLC_MayI$RedEdge,
                  fun=function(x,y){(1/x)-(1/y)})
GM1 <- overlay(E_PLC_MayI$RedEdge, E_PLC_MayI$Green,
               fun=function(x,y){x/y})
NDRE <- overlay(E_PLC_MayI$NIR, E_PLC_MayI$RedEdge,
                fun=function(x,y){(x-y)/(x+y)})

NGRDI <- function(IMG, x,y) {
        G<-IMG[[x]]
        R<-IMG[[y]]
        NGRDI<- (G-R)/(G+R)
        return(NGRDI)
}
TTVI <- function(IMG, x,y) {
        NIR<-IMG[[x]]
        R<-IMG[[y]]
        TTVI<- sqrt(abs((NIR-R)/(NIR+R)))
        return(TTVI)
}
PLC_MayI <- stack('E:/Master_Project_Data/Further_Info/Extent_PLC_MayI.tif')
NGRDI_MayI <- NGRDI(PLC_MayI, 2, 3)
TTVI_MayI <- TTVI(PLC_MayI, 5, 3)

par(mfrow = c(2,3))
plot(ARI_II, col=gray(0:100/100), main = "MicaSense - ARI_II")
plot(GM1, col=gray(0:100/100), main = "MicaSense - GM1")
plot(NDRE, col=gray(0:100/100), main = "MicaSense - NDRE")
plot(NGRDI_MayI, col=gray(0:100/100), main="MicaSense - NGRDI")
plot(TTVI_MayI, col=gray(0:100/100), main="MicaSense - TTVI")
par(mfrow = c(2,3))

hist(ARI_II, main = "Distribution of ARI_II values")
hist(GM1, main = "Distribution of GM1 values")
hist(NDRE, main = "Distribution of NDRE values")
hist(NGRDI_MayI, main = "Distribution of NGRDI values")
hist(TTVI_MayI, main = "Distribution of TTVI values")

# View histogram of data
par(mfrow = c(1,1))
hist(NDRE, main = "Distribution of NDRE values",
     xlab = "NDRE", ylab= "Frequency",
     col = "wheat", xlim = c(0, 1), breaks = 35, xaxt = 'n')
axis(side=1, at=seq(0,1, 0.05), labels=seq(0,1, 0.05))

Adicionalmente se hizo una correlacion entre todos los indices y se seleccionaron aquellos que presentaran menor correlacion entre ellos, ya que estos nos proporcionarian mas informacion al momento de observarlos y analizarlos individualmente.

index<-stack(ARI_II,GM1,NDRE,NGRDI_MayI,TTVI_MayI)
writeRaster(index, filename="Extent_MayI_INDEX.tif", format="GTiff", overwrite=TRUE)
pairs(index, main="INDEX vs INDEX")

Tambien para intentar dicriminar la vegetacion secundario y el suelo de las palmas de aceite se hizo un thresholding de la cada del Modelo de Superficie Digital (DSM) generado en el procesamiento de imagenes en el software MetaShape (Agisoft). Y en el cual se observan diferencias entre las palmas y las otras coberturas por los cambios en la altura.

# Load single raster
PLC_MayI_DSM <- raster('E:/Master_Project_Data/Muestreo_4 (May-II)/PLC_DSM_20190528_1045.tif')
E <- shapefile("E:/Master_Project_Data/Further_Info/Zona")
E_PLC_MayI_DSM <- crop(PLC_MayI_DSM, E); E_PLC_MayI_DSM

## class      : RasterLayer 
## dimensions : 2496, 3540, 8835840  (nrow, ncol, ncell)
## resolution : 0.06353734, 0.06353734  (x, y)
## extent     : 1078365, 1078590, 912173.4, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +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 
## source     : memory
## names      : PLC_DSM_20190528_1045 
## values     : 261.7552, 271.4909  (min, max)

writeRaster(E_PLC_MayI_DSM, filename="Extent_MayI_DSM.tif", format="GTiff", overwrite=TRUE)

par(mfrow=c(1,2))
plot(E_PLC_MayI_DSM, col=gray(0:100/100), main='Extent_DSM')
E_DSM_Mask <- reclassify(E_PLC_MayI_DSM, cbind(-Inf, 264, NA))
plot(E_DSM_Mask, main='Extent_DSM_Mask')

Adicionalmente se hizo el calculo de la metricas de textura para ver algun cambio entre las coberturas o definicion de las diferentes formas que pueden tomar las coberturas.

plotRGB(E_PLC_MayI, 3,2,1, stretch="lin", axes=TRUE)

plot(E_PLC_MayI, col=gray.colors(10, start = 0.0, end = 1.0, gamma = 2.3, alpha = NULL),
     axes=TRUE)

RE_glcm <- glcm(E_PLC_MayI$RedEdge, window=c(17,17), shift=c(1,1),
                 statistics=c("homogeneity","dissimilarity","entropy","second_moment"))

par(mfrow=c(2,2))
plot(RE_glcm$glcm_homogeneity, col=gray(0:100/100))
plot(RE_glcm$glcm_dissimilarity, col=gray(0:100/100))
plot(RE_glcm$glcm_entropy, col=gray(0:100/100))
plot(RE_glcm$glcm_second_moment, col=gray(0:100/100))

NIR_glcm <- glcm(E_PLC_MayI$NIR, window=c(17,17), shift=c(1,1),
                 statistics=c("homogeneity","dissimilarity","entropy","second_moment"))
par(mfrow=c(2,4))
plot(RE_glcm$glcm_homogeneity, col=gray(0:100/100))
plot(RE_glcm$glcm_dissimilarity, col=gray(0:100/100))
plot(RE_glcm$glcm_entropy, col=gray(0:100/100))
plot(RE_glcm$glcm_second_moment, col=gray(0:100/100))

plot(NIR_glcm$glcm_homogeneity, col=gray(0:100/100))
plot(NIR_glcm$glcm_dissimilarity, col=gray(0:100/100))
plot(NIR_glcm$glcm_entropy, col=gray(0:100/100))
plot(NIR_glcm$glcm_second_moment, col=gray(0:100/100))

RE_glcm_RIT <- glcm(E_PLC_MayI$RedEdge, window=c(17,17), shift=list(c(0,1),c(1,1),c(1,0),c(1,-1)),
                    statistics=c("homogeneity","contrast","dissimilarity","mean",
                              "variance","entropy","second_moment"))
plot(RE_glcm_RIT)

## Principal Component Analysis (PCA)
set.seed(1)
sr <- sampleRandom(RE_glcm_RIT, 10000)
pairs(RE_glcm_RIT, main="Texture vs Texture")

PCA <- prcomp(sr, scale=TRUE); PCA

## Standard deviations (1, .., p=7):
## [1] 2.40451819 0.93475008 0.51925966 0.21720419 0.13846386 0.09102256
## [7] 0.01640081
## 
## Rotation (n x k) = (7 x 7):
##                           PC1         PC2        PC3         PC4
## glcm_homogeneity    0.4068867 -0.09272348  0.1850938 -0.72606976
## glcm_contrast      -0.3944445  0.02673206 -0.5694936 -0.44251028
## glcm_dissimilarity -0.4069474  0.05797045 -0.3765647  0.15596832
## glcm_mean          -0.3483257 -0.55682464  0.2933321  0.09089883
## glcm_variance      -0.3525752 -0.54594373  0.2333129 -0.21419242
## glcm_entropy       -0.3820953  0.38370478  0.1823458 -0.43872957
## glcm_second_moment  0.3487315 -0.48165296 -0.5708658 -0.07782424
##                            PC5         PC6         PC7
## glcm_homogeneity    0.21661811 -0.16578754  0.43588952
## glcm_contrast       0.38463876 -0.17699195 -0.37985204
## glcm_dissimilarity  0.02015536  0.03314619  0.81449443
## glcm_mean          -0.08081637 -0.68380704  0.01363398
## glcm_variance       0.07344529  0.68672996 -0.01818201
## glcm_entropy       -0.69262708 -0.01669564 -0.03208100
## glcm_second_moment -0.55952267  0.02491634 -0.02767518

screeplot(PCA)

PCi <- predict(RE_glcm_RIT, PCA, index = 1:2)
writeRaster(PCi[[1]], filename="Extent_MayI_PC1.tif", format="GTiff", overwrite=TRUE)

par(mfrow=c(1,2))
plot(PCi[[1]], col=gray(0:100/100), main="Principal Component I")
plot(PCi[[2]], col=gray(0:100/100), main="Principal Component II")

Objetivo III

Finalmente se realizo una clasificacion supervisada teniendo en cuenta los ortomosaicos originales en valores de reflectancia, los diferentes indices multiespetrales de vegetacion y la CP1 del analisis de componentes principales de las variables de la Rotacion Inversa de Textura calculadas sin tener en cuenta la Correlacion ya que esta no proporcionaba informacion para discriminar las coberturas.

Supervised Classification

Tra_Pol<-shapefile('E:/Master_Project_Data/Further_Info/Training_Polygons.shp')
Tra_Pol<-crop(Tra_Pol, E); Tra_Pol

## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 1
## names       : Land_Cover 
## min values  :  Palm_ApSa 
## max values  :    Veg_Sec

Tra_Pol$Land_Cover

## [1] "Palm_ApSa"   "Palm_Enf"    "Palm_Erra"   "Palm_Inject" "Suelo"      
## [6] "Veg_Sec"

table(Tra_Pol$Land_Cover)

## 
##   Palm_ApSa    Palm_Enf   Palm_Erra Palm_Inject       Suelo     Veg_Sec 
##           1           1           1           1           1           1

plot(Tra_Pol, axes=T, main="Distribution of Training Sites - shapefile")

# Load raster to convert Training shapefile
Tra_Sam <- rasterize(Tra_Pol, E_PLC_MayI); Tra_Sam

## class      : RasterLayer 
## dimensions : 2492, 3534, 8806728  (nrow, ncol, ncell)
## resolution : 0.0636323, 0.0636323  (x, y)
## extent     : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 6  (min, max)
## attributes :
##        ID Land_Cover
##  from:  1  Palm_ApSa
##   to :  6    Veg_Sec

plot(Tra_Sam, axes=T, main="Distribution of Training Sites - raster")

names(Tra_Sam) <- c("Land_Cover")
Tra_Sam$Land_Cover

## class      : RasterLayer 
## dimensions : 2492, 3534, 8806728  (nrow, ncol, ncell)
## resolution : 0.0636323, 0.0636323  (x, y)
## extent     : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : Land_Cover 
## values     : 1, 6  (min, max)
## attributes :
##        ID Land_Cover
##  from:  1  Palm_ApSa
##   to :  6    Veg_Sec

# The class names and colors for plotting
table(Tra_Pol$Land_Cover)

## 
##   Palm_ApSa    Palm_Enf   Palm_Erra Palm_Inject       Suelo     Veg_Sec 
##           1           1           1           1           1           1

Class <- c("OP_ApSa", "OP_Enf", "OP_Erra", "OP_Inj", "Suelo", "Veg_Sec")
ClassDF <- data.frame(ClassValue=c(1,2,3,4,5,6), ClassNames=Class)
# Hex codes of colors
ClassColor <- c("darkgreen", "red", "darkgray", "blue", "black", "green")
# Now we ratify (RAT="Raster Attribute Table") the Class (define RasterLayer as a categorical variable)
Land_Cover <- Tra_Sam[[1]]
Land_Cover <- ratify(Land_Cover)
RAT <- levels(Land_Cover)[[1]]; RAT

##   ID
## 1  1
## 2  2
## 3  3
## 4  4
## 5  5
## 6  6

#
RAT$ID <- Class
levels(Land_Cover) <- RAT; RAT

## Warning in .checkLevels(levels(x)[[1]], value): the values in the "ID"
## column in the raster attributes (factors) data.frame have changed

## Warning in .checkLevels(levels(x)[[1]], value): NAs introducidos por
## coerción

##        ID
## 1 OP_ApSa
## 2  OP_Enf
## 3 OP_Erra
## 4  OP_Inj
## 5   Suelo
## 6 Veg_Sec

plot(Tra_Sam, axes=TRUE, col=ClassColor,
     main='Distribution of Training Samples')

# Set the random number generator to reproduce the results
set.seed(100)
# Sampling
Samp_MayI <- sampleStratified(Land_Cover, size=1000, na.rm=TRUE, sp=TRUE)
Samp_MayI

## class       : SpatialPointsDataFrame 
## features    : 6000 
## extent      : 1078365, 1078590, 912173.4, 912331.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 2
## names       :    cell, Land_Cover 
## min values  :     646,          1 
## max values  : 8806296,          6

table(Samp_MayI$Land_Cover)

## 
##    1    2    3    4    5    6 
## 1000 1000 1000 1000 1000 1000

# Plot SpatialPointsDataFrame
#install.packages("rasterVis")
plot_SPDF <- levelplot(Tra_Sam$Land_Cover, col.regions=ClassColor,
                 main='Distribution of Training Sites')
print(plot_SPDF + layer(sp.points(Samp_MayI, pch=3, cex=0.5, col=1)))

Extent_PLC_MayI <- stack('E:/Master_Project_Data/Further_Info/Extent_PLC_MayI.tif')
plot(Extent_PLC_MayI)

Extent_MayI_INDEX <- stack('E:/Master_Project_Data/Further_Info/Extent_MayI_INDEX.tif')
plot(Extent_MayI_INDEX)

Extent_MayI_PC1 <- raster('E:/Master_Project_Data/Further_Info/Extent_MayI_PC1.tif')
plot(Extent_MayI_PC1)

E_PLC_MayI <- stack(Extent_PLC_MayI, Extent_MayI_INDEX, Extent_MayI_PC1); E_PLC_MayI

## class      : RasterStack 
## dimensions : 2492, 3534, 8806728, 11  (nrow, ncol, ncell, nlayers)
## resolution : 0.0636323, 0.0636323  (x, y)
## extent     : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## names      : Extent_PLC_MayI.1, Extent_PLC_MayI.2, Extent_PLC_MayI.3, Extent_PLC_MayI.4, Extent_PLC_MayI.5, Extent_MayI_INDEX.1, Extent_MayI_INDEX.2, Extent_MayI_INDEX.3, Extent_MayI_INDEX.4, Extent_MayI_INDEX.5, Extent_MayI_PC1 
## min values :       0.003112840,       0.010345617,       0.003387503,       0.022583352,       0.115815976,        -2.934280145,         0.862184874,        -0.120839711,        -0.401420692,         0.176209277,   -10.746134280 
## max values :         0.1207294,         0.2369116,         0.2265660,         0.5106584,         0.7437858,          78.7731021,           5.4041298,           0.7680310,           0.8459403,           0.9933934,       7.6927341

# Extract the layer values for the locations
# To keep the spatial information you use `sp=TRUE` argument in the `extract` function
SampVals <- extract(E_PLC_MayI, Samp_MayI, df=TRUE)
SampVals <- SampVals[,-1] # drop the ID column
# Combine the class information with extracted values
SampData <- data.frame(ID_Class = Samp_MayI$Land_Cover, SampVals)

#############################################################################
#install.packages("rpart")
# Train the model
cart <- rpart(as.factor(ID_Class)~., data=SampData, method='class', minsplit=5)
# print(model.class)
# Plot the trained classification tree
plot(cart, uniform=TRUE, main="Classification Tree")
text(cart, cex = 0.8)

pr2011 <- predict(E_PLC_MayI, cart, type='class')
pr2011

## class      : RasterLayer 
## dimensions : 2492, 3534, 8806728  (nrow, ncol, ncell)
## resolution : 0.0636323, 0.0636323  (x, y)
## extent     : 1078365, 1078590, 912173.3, 912331.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 6  (min, max)
## attributes :
##        ID value
##  from:  1     1
##   to :  6     6

rat <- levels(pr2011)[[1]]
rat

##   ID value
## 1  1     1
## 2  2     2
## 3  3     3
## 4  4     4
## 5  5     5
## 6  6     6

rat$Legn <- ClassDF$ClassNames; rat

##   ID value    Legn
## 1  1     1 OP_ApSa
## 2  2     2  OP_Enf
## 3  3     3 OP_Erra
## 4  4     4  OP_Inj
## 5  5     5   Suelo
## 6  6     6 Veg_Sec

levels(pr2011) <- rat
levelplot(pr2011, maxpixels = 1e6,
          col.regions = ClassColor,
          main = "Decision Tree classification of MicaSense")

#install.packages("dismo")
set.seed(100)
j <- kfold(SampData, k = 5, by=SampData$ID_Class)
table(j)

## j
##    1    2    3    4    5 
## 1200 1200 1200 1200 1200

## j
##   1   2   3   4   5
## 320 320 320 320 320
x <- list()
for (k in 1:5) {
  train <- SampData[j!= k, ]
  test <- SampData[j == k, ]
  cart <- rpart(as.factor(ID_Class)~., data=train, method = 'class', minsplit = 5)
  pclass <- predict(cart, test, type='class')
  # create a data.frame using the reference and prediction
  x[[k]] <- cbind(test$ID_Class, as.integer(pclass))
}
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
conmat <- table(y)
# change the name of the classes
colnames(conmat) <- ClassDF$ClassNames
rownames(conmat) <- ClassDF$ClassNames
conmat

##          predicted
## observed  OP_ApSa OP_Enf OP_Erra OP_Inj Suelo Veg_Sec
##   OP_ApSa     625    188       3      8    23     153
##   OP_Enf      183    548       2      9    34     224
##   OP_Erra       0     13     876      4    94      13
##   OP_Inj        5     18      60    739   173       5
##   Suelo        38    140      68     71   568     115
##   Veg_Sec     117    122       0      5    18     738

# number of cases
n <- sum(conmat)
n

## [1] 6000

# number of correctly classified cases per class
diag <- diag(conmat)
# Overall Accuracy
OA <- sum(diag) / n
OA

## [1] 0.6823333

# observed (true) cases per class
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa

## [1] 0.6188

# Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc

##         producerAccuracy userAccuracy
## OP_ApSa        0.6456612        0.625
## OP_Enf         0.5325559        0.548
## OP_Erra        0.8681863        0.876
## OP_Inj         0.8839713        0.739
## Suelo          0.6241758        0.568
## Veg_Sec        0.5913462        0.738

Conclusiones

• Los arboles de decision o Deision Trees (DT) nos puede discriminar la palma de aceite y la cobertura seundaria con un indice kappa de 0.62

• Hay que realizar analisis mas detallado del objeto de estudio para poder determinar alguna diferencia significativa entre las palmas aparentemente sanas y palmas enfermas