Introducción

La información de cobertura de la tierra es uno de los componentes más importantes en varios estudios ambientales como: cambios en los sistemas ecológicos (Vitousek, 1994), cambio climático (Skole, 1994), calentamiento global, entre otros. Es una variable dinámica ya que refleja la interacción entre las actividades socioeconómicas y los cambios ambientales regionales, siendo necesario esta información actualizarla con frecuencia(Rujoiu-Mare & Mihai, 2016). En los últimos años, entre las principales aplicaciones de la observación de la tierra mediante imágenes satélites se encuentra el mapeo y monitoreo de las coberturas terrestres la cual permite conocer cambios del suelo durante diferentes periodos de tiempo (Rodriguez-Galiano, Ghimire, Rogan, Chica-Olmo, & Rigol-Sanchez, 2012), el avance de la teledetección permite cada vez más adquirir mediciones de la superficie terrestre a diversas escalas espaciales y temporales(Huang, C., L. S. Davis, 2002). El uso de imágenes satelitales de resolución media ha sido utilizado en diferentes estudios de clasificación de coberturas, sin embargo, en escalas espaciales finas, el tamaño de pixel puede llegar a ser más grande que el elemento mapeado en la cobertura, causando una mezcla de clases (Blaschke et al., 2004)

En consecuencia, se han desarrollado diversos métodos que permiten mejorar la precisión de la clasificación de la cobertura del suelo mediante el uso de imágenes satélites, uno de las más novedosas en el uso de imágenes o fotografías satelitales o aéreas es el basado en objetos geográficos (GEOBIA por sus siglas en inglés), que es ampliamente utilizado en diferentes tipos de investigación (Cánovas, 2012; García-Pedrero et al, 2015; Perea et al, 2009; Vásquez, 2011).

La clasificación basada en objetos se fundamenta principalmente en la capacidad de identificar los diferentes elementos en una imagen, estos se representan por agrupaciones de pixeles con las mismas características que forman finalmente al objeto (Chen et al, 2018). La interpretación de un conjunto de pixeles que conforman un objeto dentro de la imagen (GEOBIA) representa de mejor forma a los componentes en la escena, estos a su vez son dependientes de la escala (Blaschke et al, 2014). Este sistema se ha convertido en una impresionante herramienta de análisis y ha superado las complicaciones que se generaban en torno a los análisis de pixeles y ha resultado en una interpretación más precisa (Blaschke et al, 2014; Chen et al 2018). GEOBIA se ha convertido en una de las herramientas más utilizadas en diversos campos de investigación y ha ayudado a resolver problemas y necesidades de diversas índoles.

En el análisis espacial, las técnicas y herramientas utilizadas en combinación nos permiten mejorar nuestra percepción de las áreas estudiadas. Una herramienta fundamental para esto es la clasificación mediante la técnica de “Random Forest” (Breiman, 2001), usada para clasificar y agrupar un conjunto de datos de forma aleatoria convirtiéndose en un excelente sistema para la clasificación, regresión y estimación (Liaw & Wiener, 2002; Louppe 2014). El método de Breiman selecciona en cada nodo un grupo de variables, utilizadas para dividir, en base a la mejor división a partir de estos nodos se construyen árboles y así consecutivamente (Biau, 2012).

Random forest es uno de los clasificadores supervisados más utilizado en solución de problemas de reconocimiento de patrones (Loyola-Gonazales, 2012), consiste en una combinación de clasificadores de árbol donde cada clasificador se genera usando un vector aleatorio muestreado independientemente del vector de entrada, y cada árbol arroja un voto unitario para la clase más popular para clasificar un vector de entrada (Breiman 1999).

Entre las ventajas del uso de random forest (Breiman, 2001)

El objetivo de este estudio es evaluar la exactitud de la clasificación basada en objetos utilizando el metodo Random forest, utilizando la guia Ned Horning https://bitbucket.org/rsbiodiv/

Metodología

Area de estudio

La zona de estudio se encuentra ubicada en la región de los Andes Sur del Ecuador. Comprende parte del canton Cuenca, tiene una superficie de aproximandamente 80 km2 con una altitud que va desde los 2500 a 3600 msnm. A escala muy general presenta 5 coberturas principalmente: pajonal del páramo, bosque alto andino, pastizal, agua (lagunas) e infraestructura.(Fig 1)

Obtención de datos

La imagen fue descargada de la plataforma glovis de la NASA https://glovis.usgs.gov/app. La imagen es una Landsat 8, con código LC80100622016325LGN01 capturada el 21 de noviembre del 2016. Se utilizaron las bandas entre la 1 y la 7 y se recortaron con la zona de estudio. Los niveles digitales de las bandas fueron convertidos a valores de reflectancia utilizando la libreria Landsat8 de RStudio.

Se utilizó la guia de Ned Horning y se siguió los pasos del tutorial de Orfeo toolbox y los scripts de RStudio para la obtención de los datos. https://bitbucket.org/rsbiodiv/

Para la segmentación de la imagen se utilizó la herramiento orfe toolbox. Esta herramienta permite generar grupos continuos de píxeles con caracteristicas multiespectrales similares de la imagen, se utilizó el algoritmo mean-shift y por capacidad del computador en procesos posteriores se utilizaron los siguientes parámetros:

  • Spatial Radius: 6

  • Range radius: 10

  • Mode covergence threshold: 0.1

  • Maximum number of iterations: 100

  • Minimum region size: 60

Se obtuvo una capa vectorial con los segmentos, la cual fue transformada a formato raster con un tamaño de pixel de 30 m. Se calculó las siguientes estadisticas de las imagen por segmento:

  • Promedio para cada banda
  • Desviación estandar para cada banda
  • coeficiente de variación para cada banda.

Este procedimiento implica gran poder computacional, en el caso de este tutorial el tiempo de procesamiento para obtener las estadisticas fue alrededor de 8 horas, por lo cual solo se mostrara el codigo y no se correrá como parte del script. Las estadisticas se presentan en un formato CSV

Codigo

Imagen mutiespectral
satImage <- "D:/maestria/percepcion remota avanzada/cla_landsat8/imagen.tif"

Se puede utilizar en formato Raster o vector la capa de segmentos. 
segVector <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/segmentada.shp'
segImage <- ""

Nombre del campo de atributo que tiene la ID de los segmentos. 
idAttributeLabel <- "DN"

Nombre del archivo CSV de salida que contendra las estadisticas por segmento
outFeatures <- "D:/maestria/percepcion remota avanzada/cla_landsat8/segmentFeatures.csv"

ndSat <- 0
""
ndSeg <- -1

Se define las funciones para calcular características de los segmentos
Función para calcular la media después de recortar el 20% superior e inferior de los valores
clippedSDFun <- function(nums,...) {
a=quantile(nums, probs=seq(0,1,0.2), na.rm = TRUE)
sd(nums[which(nums >= a[2] & nums <= a[5])])
}

Funcion para calcular el coef de variacion
clippedCVFun <- function(nums,...) {
a=quantile(nums, probs=seq(0,1,0.2), na.rm = TRUE)
as.double(cv(nums[which(nums >= a[2] & nums <= a[5])]))
}

Funcion para calcular el promedio
medianFun <- function(nums,...) {
median(as.double(nums), na.rm = TRUE)
}

Funcion para crear etiquetas de columnas con el mismo numero de bandas de la imagen
createColLabels <- function(numBands, baseName) {
labels <- "segnumber"
for (i in 1:numBands) {
labels[i+1] <- paste(baseName, "Band", i, sep='')
}
labels
}

Cargar la imagen multiespectral
satelliteImage <- stack(satImage)
for (b in 1:nlayers(satelliteImage)) { NAvalue(satelliteImage@layers[[b]]) <- ndSat }

Las siguientes funciones calcula las estadisticas por cada segmento. Se puede utilizar la capa de segmentos en formato .shp o .tif, el codigo esta generado para cualquiera de los dos, parte del codigo se bloquea segun el formato que se utiliza.

CUANDO SE UTILIZA LOS SEGMENTOS EN FORMATO VECTOR

if (segVector != '' & segImage == '') { 
cat("Reading vector segments\n")
vec <- readShapePoly(segVector)
allAtt <- slot(vec, "data")
numRows <- nrow(allAtt)

zonalMean <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)  
zonalSD <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)
zonalClippedMean <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)
zonalClippedSD <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)
zonalMedian <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)
zonalCV <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)
zonalClippedCV  <- matrix(nrow=0, ncol=nlayers(satelliteImage)+1)

for (x in 1:numRows) {
for (x in 1:5) {
segPoly<- vec[vec[[idAttributeLabel]]==allAtt[x,idAttributeLabel],]

zonalMean <- rbind(zonalMean, c(x, extract(satelliteImage, segPoly, fun=mean, na.rm = TRUE)))
zonalSD <- rbind(zonalSD, c(x,extract(satelliteImage, segPoly, fun=sd, na.rm = TRUE)))
zonalClippedMean <- rbind(zonalClippedMean, c(x,extract(satelliteImage, segPoly, fun=clippedMeanFun, na.rm = TRUE)))
zonalClippedSD <- rbind(zonalClippedSD, c(x,extract(satelliteImage, segPoly, fun=clippedSDFun, na.rm = TRUE)))
zonalMedian <- rbind(zonalMedian, c(x,extract(satelliteImage, segPoly, fun=medianFun, na.rm = TRUE)))
zonalCV <- rbind(zonalSD, c(x,extract(satelliteImage, segPoly, fun=cv, na.rm = TRUE)))
zonalClippedCV <- rbind(zonalClippedMean, c(x,extract(satelliteImage, segPoly, fun=clippedCVFun, na.rm = TRUE)))
}

zonalMean <- as.data.frame(zonalMean)
zonalSD <- as.data.frame(zonalSD)
zonalClippedMean <- as.data.frame(zonalClippedMean)
zonalClippedSD <- as.data.frame(zonalClippedSD)
zonalMedian <- as.data.frame(zonalMedian)
zonalCV <- as.data.frame(zonalCV)
zonalClippedCV <- as.data.frame(zonalClippedCV)

Volver a etiquetar los encabezados de las columnas de la tabla de datos
meanLab <- createColLabels(nlayers(satelliteImage), "mean")
colnames(zonalMean) <- meanLab
sdLab <- createColLabels(nlayers(satelliteImage), "sd")
colnames(zonalSD) <- sdLab
clipMeanLab <- createColLabels(nlayers(satelliteImage), "clpMean")
colnames(zonalClippedMean) <- clipMeanLab
clipSDLab <- createColLabels(nlayers(satelliteImage), "clpSD")
colnames(zonalClippedSD) <- clipSDLab
medianLab <- createColLabels(nlayers(satelliteImage), "median")
colnames(zonalMedian) <- medianLab
cvLab <- createColLabels(nlayers(satelliteImage), "cv")
colnames(zonalCV) <- cvLab
clipCVLab <- createColLabels(nlayers(satelliteImage), "clpCV")
colnames(zonalClippedCV) <- clipCVLab

CUANDO SE UTILIZA LOS SEGMENTOS EN FORMATO RASTER
} else if (segVector == '' & segImage != '') {
cat("Reading raster segments\n")
segmentImage <- raster(segImage)
NAvalue(segmentImage) <- ndSeg

Preparar valores de segmento para cálculos de estadísticas
vals <- getValues(satelliteImage)
zones <- round(getValues(segmentImage), digits = 0)
rasterDT <- data.table(vals, z=zones)
setkey(rasterDT, z)

Calcular las estadisticas
zonalMean <- rasterDT[, lapply(.SD, "mean", na.rm = TRUE), by=z]
zonalSD <- rasterDT[, lapply(.SD, "sd", na.rm = TRUE), by=z]
zonalClippedMean <- rasterDT[,lapply(.SD, clippedMeanFun), by =z]
zonalClippedSD <- rasterDT[,lapply(.SD, clippedSDFun), by =z]
zonalMedian <- rasterDT[, lapply(.SD, medianFun), by=z]
zonalCV <- rasterDT[, lapply(.SD, "cv", na.rm = TRUE), by=z]
zonalClippedCV <- rasterDT[,lapply(.SD, clippedCVFun), by =z]

Volver a etiquetar los encabezados de las columnas de la tabla de datos
meanLab <- createColLabels(nlayers(satelliteImage), "mean")
setnames(zonalMean, 1:(nlayers(satelliteImage)+1), meanLab)
sdLab <- createColLabels(nlayers(satelliteImage), "sd")
setnames(zonalSD, 1:(nlayers(satelliteImage)+1), sdLab)
clipMeanLab <- createColLabels(nlayers(satelliteImage), "clpMean")
setnames(zonalClippedMean, 1:(nlayers(satelliteImage)+1), clipMeanLab)
clipSDLab <- createColLabels(nlayers(satelliteImage), "clpSD")
setnames(zonalClippedSD, 1:(nlayers(satelliteImage)+1), clipSDLab)
medianLab <- createColLabels(nlayers(satelliteImage), "median")
setnames(zonalMedian, 1:(nlayers(satelliteImage)+1), medianLab)
cvLab <- createColLabels(nlayers(satelliteImage), "cv")
setnames(zonalCV, 1:(nlayers(satelliteImage)+1), cvLab)
clipCVLab <- createColLabels(nlayers(satelliteImage), "clpCV")
setnames(zonalClippedCV, 1:(nlayers(satelliteImage)+1), clipCVLab)

Union los resultados en una sola tabla de variables 
meanSDMerge <- merge(zonalMean, zonalSD)
clippedmeanSDMerge <- merge(zonalClippedMean, zonalClippedSD)
trainValsTemp1 <- merge(meanSDMerge, clippedmeanSDMerge)
trainValsTemp2 <- merge(trainValsTemp1, zonalMedian)
cvMerge <- merge(zonalCV, zonalClippedCV)
trainVals <- merge(trainValsTemp2, cvMerge)

Guardar las estadisticas en un archivo CSV
write.csv(trainVals, file=outFeatures, row.names=FALSE)

En este estudio tomaron 179 muestas (aleatorios estratificados), y a partir de estas se crearon poligonos de entrenamiento. La tabla 1 muestra el numero de muestras por cobertura y el codigo asignado.

Cobertura # poligonos Codigo
Pajonal de paramo 44 1
infraestructura 40 2
Agua 15 3
Bosque alto andino 37 4
Pastizales 43 5

Los segmentos que no fueron clasificados correctamente se reasignaron editandolos. Se generaron puntos en los segmentos mal clasificacos con etiquetas de la clase correcta.

Analisis

Se evaluo la exactitud de la clasificacion basada en objetos, es decir el grado de concordancia entre las clases asignadas por el clasificador y los datos del terreno asignados de referencia, mediante la matriz de confusion

La matriz de confusion es es una matriz cuadrada de n x n, donde n es el número de clases, esta muestra la relación entre dos series de medidas correspondientes al área en estudio; las columnas corresponden a los datos de referencia, y las filas corresponden a las asignaciones del clasificador. La figura muestra un ejemplo de matriz de confusion e indices utilizados para la evalucion

Resultados

#Cargar las librerias
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.4.4
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: D:/Lore/Documents/R/win-library/3.4/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: D:/Lore/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
library(sp)

#Variables
#CSV de las estadisticas de las bandas por segmento
segCsv <- "D:/maestria/percepcion remota avanzada/cla_landsat8/segmentFeatures.csv"

# Segmentos en formato raster 
segImage <- "D:/maestria/percepcion remota avanzada/cla_landsat8/seg_raster.tif"

# Segmento raster con valores nodata 
nd <- 1

# Imagen clasificada con random fores
outImage <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/RF_2.tif'

# Nombre del puntos de evaluacion de exactitud 
outMarginFile <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/puntos_exactitud_2.shp'

#Poligonos de entrenamiento
trainingDsn <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/train.shp'

#nombre del archivo de los poligonos de entrenamiento 
trainingLayer <- 'train'

#nombre de la columna con el codigo de la categoria de clasificacion  
classNum <- "Id"

# Salida con informacion de la clasificacion
outClassMapCSV <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/classmapping_2.csv'

#parametros
classImage <- TRUE
probImage <- TRUE
threshImage <- TRUE
probThreshold <- 75

# Leer los poligonos de entrenamiento 
vec <- readOGR(trainingDsn, trainingLayer)
## OGR data source with driver: ESRI Shapefile 
## Source: "D:/maestria/percepcion remota avanzada/cla_landsat8/train.shp", layer: "train"
## with 179 features
## It has 1 fields
pts <- slot(vec, "data")

# Cargar los segmentos
segImg <- raster(segImage)

#Extraer los IDs de los segmentos segun las caracteristicas de los poligonos
#Extraiga los ID de segmento bajo las características de punto, línea o polígono
exSegNums <- extract(segImg, vec, cellnumbers=TRUE)
if (is.matrix(exSegNums)) {
  exSegNums <- as.list(data.frame(t(exSegNums)))
}

# Eliminar los valores NULL 
exSegNums <- exSegNums[!sapply(exSegNums, is.null)]

#Seleccione las ID de segmento únicas debajo de las características de cada vector y asocie la respuesta
trainSegs <- matrix(ncol=3, nrow=0)
for (i in 1:length(exSegNums)) {
  lineResponse <- pts[i,classNum]
  if (is.matrix(exSegNums[[i]])) { 
    segNums <- exSegNums[[i]][which(duplicated(exSegNums[[i]][,2]) == FALSE),]
  } else {
    segNums <- exSegNums[[i]]
  }
  
  if (is.matrix(segNums)) {
    trainSegs <- rbind(trainSegs, cbind(lineResponse, segNums))
  }
  else {
    trainSegs <- rbind(trainSegs, cbind(lineResponse, segNums[1], segNums[2]))
  }  
}

# Eliminar los nombres de las filas y agregar nombres a las columnas 
rownames(trainSegs) <- NULL
colnames(trainSegs) <- c("response", "cellNums", "segNums")

# Elimine los valores de NA de la lista de ID de segmento únicos
trainSegs_no_na <- as.data.frame(na.omit(trainSegs))

# Leer CSV de los segmentos
segAtr <- read.csv(segCsv, header=TRUE)

# Eliminar los NAs de las tabla 
segAtr <- na.omit(segAtr)

#Crear un conjunto de datos de entrenamiento haciendo coincidir las ID de segmento de entrenamiento con la información de características de los segmento 
train <- segAtr[match(trainSegs_no_na$segNums, segAtr$segnumber),]

# Eliminar NAs
train_no_na <- as.data.frame(na.omit(train))

# crear variable de respuesta
response_no_na <- trainSegs_no_na[match(train_no_na$segnumber, trainSegs_no_na$segNums), c(1,3)]

# Correr el algoritmo de random forest
randfor <- randomForest(as.factor(response_no_na$response) ~. , data=train_no_na[,-1], proximity=TRUE)

# Guardar los resultados.
bs <- blockSize(segImg)

# Calcular cuántas bandas debe tener la imagen de salida.
numBands <- classImage + probImage + threshImage

# Raster de salida
img.out <- brick(segImg, values=FALSE, nl=numBands)
img.out <- writeStart(img.out, outImage, overwrite=TRUE, datatype='INT1U')

# Predicciones
predValues <- predict(randfor, segAtr, type='response')
predValuesDF <- data.frame(segAtr$segnumber, predValues)

if (probImage || threshImage) { 
  predProbs <- predict(randfor, segAtr, type='prob')
  maxProb <- round(apply(predProbs, 1, max) * 100)
  maxProbDF <- data.frame(segAtr$segnumber, maxProb)
}

# Nuevo valor clasificado.
for (i in 1:bs$n) {
  cat("processing block", i, "of", bs$n, "\r")
  img <- getValues(segImg, row=bs$row[i], nrows=bs$nrows[i])
  outMatrix <- matrix(nrow=length(img), ncol=0)
  
  # Cambiar los datos vacios por NAs
  is.na(img) <- img == nd
  if (classImage) {
    # Convierta el ID de segmento a la clase predicha para poder establecer un valor de nodata.
    outMatrix <- cbind(outMatrix, predValuesDF$predValues[match(img, predValuesDF$segAtr.segnumber)])
  }
  if (probImage) { 
    outMatrix <- cbind(outMatrix, maxProbDF$maxProb[match(img, maxProbDF$segAtr.segnumber)])
  }
  if (threshImage) {
    threshValues <- as.numeric(predValuesDF$predValues[match(img, predValuesDF$segAtr.segnumber)])
    threshValues[which(maxProbDF$maxProb[match(img, maxProbDF$segAtr.segnumber)] <= probThreshold)] <- 0
    outMatrix <- cbind(outMatrix,threshValues)
  }
  writeValues(img.out, outMatrix, bs$row[i])
}
## processing block 1 of 4 
processing block 2 of 4 
processing block 3 of 4 
processing block 4 of 4 
# Guardar el mapa de clasificacion 
img.out <- writeStop(img.out)


# Guardar en CSV datos de la clasificacion
if (outClassMapCSV != "") {
  write.csv(predValuesDF, file=outClassMapCSV, row.names=FALSE)
}

# Diagrama de puntos de importancia variable medido de Random Forest
varImpPlot(randfor)

# tasa de error y matriz de confusion
confMatrix <- randfor$confusion
cat("\n#################################################################################\n")
## 
## #################################################################################
cat("OOB error \n", 1 - (sum(diag(confMatrix)) / sum(confMatrix[,1:ncol(confMatrix)-1])), "%\n\n", sep="")
## OOB error 
## 0.01117318%
cat("matriz de confusion\n")
## matriz de confusion
print(randfor$confusion)
##    1  2  3  4  5 class.error
## 1 44  0  0  0  0   0.0000000
## 2  0 40  0  0  0   0.0000000
## 3  2  0 13  0  0   0.1333333
## 4  0  0  0 37  0   0.0000000
## 5  0  0  0  0 43   0.0000000
cat("\n")
if (outMarginFile != "") {
  # Calcular el margen (proporción de votos para la clase correcta menos la proporción máxima de votos para otras clases)
  marginData <- margin(randfor)
  trainingAccuracy <- cbind(marginData[order(marginData)], trainSegs_no_na[order(marginData),])
  
   # Agregar nombres a la columna de la tabla de atributos
  colnames(trainingAccuracy) <- c("margin", "classNum", "cellNum", "segID")
  
  #  Calcular coordenadas xy de los puntos de entrenamiento
  xyCoords <- matrix(ncol=2, nrow=0)
  for (z in 1:nrow(trainingAccuracy)) {
    xyCoords <- rbind(xyCoords, xyFromCell(segImg, trainingAccuracy[z,3]))
  }
  
  #   Crear y guardar los puntos de para mejorar los datos de entrenamiento
  pointVector <- SpatialPointsDataFrame(xyCoords, trainingAccuracy[, c(1,2,4)], coords.nrs = numeric(0), proj4string = segImg@crs)
  writeOGR(pointVector, outMarginFile, "layer", driver="ESRI Shapefile", check_exists=TRUE)
}

  cols <- c("aliceblue",  "gray47","navy","chartreuse1", "darkgreen")

  
  plot(img.out$RF_2.1,col=cols,main = "Clasificacion de cobertura Random Forest", legend=TRUE)

plot(img.out$RF_2.1, col=cols,main = "Puntos de verificacion de exactitud")
plot(pointVector, add=TRUE)

# Edicion de los segmentos mal clasificados

#LIbrerias
require(raster)
require(rgdal)


# Segmentos en formato raster
segImage <- "D:/maestria/percepcion remota avanzada/cla_landsat8/seg_raster.tif"
nd <- 1

# Nombre de la capa con los segmentos con los segmentos mal clasificados corregidos 
outImage <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/nueva_clasificacion.tif'

# classmapping CSV 
inClassMapCSV <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/classmapping.csv'

#Puntos reclasificados
editPointsDsn <- 'D:/maestria/percepcion remota avanzada/cla_landsat8/nuevos_puntos.shp'

# Nombre de la capa 
editPointsLayer <- 'nuevos_puntos'

# nombre de la columna con el codigo de la categoria de reclasificacion 
newClassNum <- "Id"

# Leer los puntos 
vec <- readOGR(editPointsDsn, editPointsLayer)
## OGR data source with driver: ESRI Shapefile 
## Source: "D:/maestria/percepcion remota avanzada/cla_landsat8/nuevos_puntos.shp", layer: "nuevos_puntos"
## with 69 features
## It has 1 fields
newClass <- slot(vec, "data")

# Cargar los segmentos
segImg <- raster(segImage)

#abrir el csv
classMap <- read.csv(inClassMapCSV)

# Extraer ID de segmento bajo  los puntos para editar
cat("Extracting segment IDs under the points to modify\n")
## Extracting segment IDs under the points to modify
changeSegIDs <- cbind(newClass[,newClassNum], extract(segImg, vec))
for (i in 1:nrow(changeSegIDs)) {
  #classMap[changeSegIDs[i, 2], 2] <- changeSegIDs[i,1]
  classMap[which(classMap[,1]==changeSegIDs[i, 2]), 2] <- changeSegIDs[i,1]
}
 
# Guardar el nuevo raster.
# cargar el raster.
bs <- blockSize(segImg)

# Crea el ráster de salida 
img.out <- raster(segImg)
img.out <- writeStart(img.out, outImage, overwrite=TRUE, datatype='INT1U')

for (i in 1:bs$n) {
    img <- getValues(segImg, row=bs$row[i], nrows=bs$nrows[i])
    is.na(img) <- img == nd
  img.match <- as.numeric(classMap$pred[match(img, classMap$segAtr.segnumber)])
  img.match[is.na(img.match) == TRUE] <- NAvalue(img.out)
  writeValues(img.out, img.match, bs$row[i])
}

img.out <- writeStop(img.out)

#PLOTEAR EL MAPA DE CLASIFICACION FINAL
cols <- c("aliceblue",  "gray47","navy","chartreuse1", "darkgreen")
plot(img.out,col=cols,main = "Rectificacion de clasificacion de cobertura Random Forest", legend=TRUE)

Discusion

Este trabajo puede ser visto como una prueba inicial del uso del clasificador de random forest en la zona de estudio, el método en general tuvo buenos resultados en la clasificación, obteniendo un valor de error 0BB igual 0.0055%. Todas las coberturas a excepcion de agua no presentaron error de clasificacion.

Es importante recalcar que el uso de imágenes satélites y software libre, permite generar información de gran calidad a bajo costo, muy importante en zonas donde los presupuestas de investigación y generación de información base es muy reducido