Clasificación Basada en Objetos

Introducción

En este taller se abordará la clasificación basada en objetos, usando segmentación y Random Forest (RF). Para la elaboración de este taller se tendrá en cuenta la guía de Ned Horning _ (American Museum of Natural History) _ titulada Training Guide to Classify Satellite Images using Segmentation and Random Forests.

Esta técnica usa segementos de la imagen en lugar de pixeles individuales, el resultado de este tipo de clasificación tiene menos speckle (efecto de sal y pimienta), así mismo tiene la ventaja de ser mucho más rápida que la clasifición por pixeles, ya que el numero de los segmanetos es mucho menor a los pixeles en una imagen.

Metodología

Este taller se llevará a cabo en tres tiempos:

  • Segmentación de la imagen
  • Cálculo de las características de la segmentación
  • Clasificación de los segmentos en la Imagen

Para esto se usa el software RStudio for R, QGIS 3.4, Orfeo Toolbox (OTB). En RStudio se deben instalar las siguientes librerias para lo cual se puede usar el siguente código:

#install.packages("raster")
#install.packages("data.table")
#install.packages("maptools")
#install.packages("ramdomForest")
#install.packages("rgdal")
#install.packages("sp")

** Segmentación

Como indiqué antes, la segmentación se lleva acabo con OTB, para esto se descargó la versión 6.6 para Windows X64, y se utilizó la GUI de Monteverdi, con el algoritmo Mean-Shift, con los parametros que se muestran en la siguiente imagen:

· Calcular las características de los segmentos

Estas caracteristicas nos ayudan a escoger las variables que podrían ser usadas para crear un modelo RF. Las caracteríscas podrían ser la media de pixeles por bandas o las formas acarcaterísticas de un área; estas característcas además podrísan contener información de la relación de los segmentos con sus vecinos.

El Script que voy a aplicar en este ejercicio fue escrito por Ned Horning y tiene la capacidad de calcular las siguientes características estadísticas:

  • Media para cada banda de la imagen
  • Desviación estándar para cada banda
  • Coeficiente de variación para cada banda
  • Media, desviación estándar, y coeficiente de varaiación para cada banda luego de remover el 40% de los pixeles
  • Mediana de cada banda
  1. Llamar las librerias necesarias:
library(raster)
## Loading required package: sp
library(data.table)
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:raster':
## 
##     shift
library(maptools)
## Checking rgeos availability: TRUE
library(sp)
library(rgeos)
## rgeos version: 0.3-28, (SVN revision 572)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE
  1. Insertar las variables
setwd('D:/R-Scrips/')
satImage <- "./data-raw/mocoa.tif" #Imagen que será usada para calcular las características
segVector <- "" #Vector de segmentación
idAttributeLabel <- "DN"
segImage <- "./Mocoa/MeanShift_20_15_100.tif" #Raster Segmentado
outFeatures <- "./Mocoa/segmentFeatures_SubRas3.csv" #Ruta para la salida del archivo CSV
ndSat <- 0 #Rellenar datos vacios o nulos en el raster original
ndSeg <- -1 #Rellenar datos vacios o nulos en el raster segmentado
mrst <- brick(satImage)
srst <- raster(segImage)
plotRGB(mrst, r=1, g=2, b=3, scale(10000), stretch="lin", axes=TRUE, main="satImage")

plot(srst)

3. Inicio del proceso

startTime <- Sys.time()
cat("Start time", format(startTime),"\n")
## Start time 2018-11-21 16:47:39
## Define las funciones para calcular las características para los segmentos
# Función para calcular la media luego de quitar el 20% de los valores superiores e inferiores
clippedMeanFun <- function(nums,...) {
  a=quantile(nums, probs=seq(0,1,0.2), na.rm = TRUE)
  mean(nums[which(nums >= a[2] & nums <= a[5])])
}

# Función para calcular la desviación estándar luego de quitar el 20% de los valores superiores e inferiores
clippedSDFun <- function(nums,...) {
  a=quantile(nums, probs=seq(0,1,0.2), na.rm = TRUE)
  sd(nums[which(nums >= a[2] & nums <= a[5])])
}

# Función para calcular el coeficiente of variación luego de quitar el 20% de los valores superiores e inferiores
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])]))
}

# Función para calcular la mediana debido a que la librería data.table tiene problemas con datos enteros
medianFun <- function(nums,...) {
  median(as.double(nums), na.rm = TRUE)
}

# Función para crear los nombres de columnas para que coincida con el número 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 imagen que será usada para crear las características de los segmentos
cat("Reading satellite image\n")
## Reading satellite image
satelliteImage <- stack(satImage)
for (b in 1:nlayers(satelliteImage)) { NAvalue(satelliteImage@layers[[b]]) <- ndSat }

## Determinar si será usado per-segment o la imagen completa para el proceso
# Este bloque se procesará si se definió la ruta y nombre segVector (El vector de segmentación)
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],]
    
    # Calcular las estadísticas zonales
    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)))

    cat("Processing segment", x, "of",numRows, "\r")
  }
  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)
  cat("\nFinished processing polygons\n")
  
  # Renombrar los nombres de las cabeceras de la tabla
  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

# Este bloque se procesará si se definió la ruta y nombre segImage (El raster de segmentación)
} else if (segVector == '' & segImage != '') {
  cat("Reading raster segments\n")
  segmentImage <- raster(segImage)
  NAvalue(segmentImage) <- ndSeg
  
  # Prepare segment values for zonal statistics calculations
  cat("Calculating zonal statistics\n")
  vals <- getValues(satelliteImage)
  zones <- round(getValues(segmentImage), digits = 0)
  rasterDT <- data.table(vals, z=zones)
  setkey(rasterDT, z)
  
  # Calcula las estadíticas zonales
  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]
  
  # Renombrar los nombres de las cabeceras de la tabla
  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)
 
# Si los parámetros segImage y/o segVector no están configurados muestra un mensaje y detenga el procesamiento
} else {
  cat("\n\n*****************************************************************************")
  cat("\nUna y sólo una de las variables segImage o segVector tienen una ruta\n")
  cat("Y la otra variable debe estar en comillas (simple o doble) sin espacio entre ellas\n\n")
  stop("Cambie los valores para segImage o segVector y reinicie el script\n\n", call.=FALSE)
}
## Reading raster segments
## Calculating zonal statistics
# Unir los resultados en una tabla para las variables de entrenamiento
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 CSV
write.csv(trainVals, file=outFeatures, row.names=FALSE)
# Calcular el tiempo del proceso
timeDiff <- Sys.time() - startTime
cat("Processing time", format(timeDiff), "\n")
## Processing time 1.169052 mins

· Clasificación de los segmentos de la Imagen

Para clasificar los segmentos es necesario seleccionar los datos de entrenamiento para informar al modelo de RF como debería estar clasificado. Lo datos de entrenamiento deben estar en aurchivo vector.

  1. Para esto se crea un shape de puntos en el que se nombran los tipos de cobertura, y se correrá el script creado por Ned Horning llamado SegmentRF_Classification.

  2. Llamar las librerías necesrias:

library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
library(raster)
library(rgdal)
## rgdal: version: 1.3-3, (SVN revision 759)
##  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: D:/aamejial/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: D:/aamejial/Documentos/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(sp)
  1. Fijar las variables:
setwd('D:/R-Scrips/')
#Nombre y Ubicación del CSV con las características de los segmentos
segCsv <- "./Mocoa/segmentFeatures_SubRas3.csv"
segImage <- "./Mocoa/MeanShift_20_15_100.tif" #Raster Segmentado
nd <- 1 # valor que van a tomar los datos nulos
# Nombre y ruta de la salida de la imagen clasificada
outImage <- "./Mocoa/Forest_NonforestMapTest_v9.tif"

# Nombre y ruta donde se guardará Shapefile de puntos que será creado. NOTA: Si este archivo existe fallará la escritura y arrojará el mensaje "Error en creación de archivo de salida"
# Si este archivo no es necesario deje las comillas dobles o simples sin nada en el medio _(“” o '')_
outMarginFile <- './Mocoa/testAccuracyPoints8.shp'

#Nombre y ruta del shape que contiene los datos de entrenamiento
trainingDsn <- "./Mocoa/Training_BKP.shp"

# Nombre de la capa con el archivo vector este a menudo es un archovo sin nunga extensión. 
trainingLayer <- 'Training_BKP'

# Ingrese el número de columna o el nombre de la columna que contene las clases 
classNum <- "type_id"

#  Salida del archivo csv con class infromación de las clases. Si no se necesita este archivo se pueden ingresar un par de comillas sin espacio entre ellas (“” o '')
outClassMapCSV <- "./Mocoa/classMapping_v2.csv"

classImage <- TRUE
probImage <- TRUE
threshImage <- TRUE
probThreshold <- 75
  1. Empezar el Proceso
## Hora de inicio del Proceso
startTime <- Sys.time()
cat("Hora de Inicio", format(startTime),"\n")
## Hora de Inicio 2018-11-21 16:48:49
# Leer el archivo Vector 
cat("Leyendo el archivo vector\n")
## Leyendo el archivo vector
vec <- readOGR(trainingDsn, trainingLayer)
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R-Scrips\Mocoa\Training_BKP.shp", layer: "Training_BKP"
## with 267 features
## It has 3 fields
## Integer64 fields read as strings:  type_id id
pts <- slot(vec, "data")
# Cargar el Raster de la imagen segmentada 
segImg <- raster(segImage)

# Extraer los segmentos del Shapefile
cat("Extracting segment IDs under the vector features\n")
## Extracting segment IDs under the vector features
exSegNums <- extract(segImg, vec, cellnumbers=TRUE)
if (is.matrix(exSegNums)) {
  exSegNums <- as.list(data.frame(t(exSegNums)))
}

# Remover los valores nulos
exSegNums <- exSegNums[!sapply(exSegNums, is.null)]
# Asociar el identificador a cada capa
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 el nombre de las columnas
rownames(trainSegs) <- NULL
colnames(trainSegs) <- c("response", "cellNums", "segNums")

# Eliminar los valores nulos de la lista de segmentos
trainSegs_no_na <- as.data.frame(na.omit(trainSegs))

# Leer el archivo csv
segAtr <- read.csv(segCsv, header=TRUE)

# Eliminar los valores nulos de la tabla
segAtr <- na.omit(segAtr)

# Crear un set de entrenamiento que coincida con con la ID del cada segmneto del  set de entreamiento
train <- segAtr[match(trainSegs_no_na$segNums, segAtr$segnumber),]

# Remover los valores nulos del set de entrenamiento
train_no_na <- as.data.frame(na.omit(train))

# Crea una respuesta a la variable del data frame
response_no_na <- trainSegs_no_na[match(train_no_na$segnumber, trainSegs_no_na$segNums), c(1,3)]
cat("Iniciando el cálculo \n")
## Iniciando el cálculo
#Inicio del algoritmo _RF_
randfor <- randomForest(as.factor(response_no_na$response) ~. , data=train_no_na[,-1], proximity=TRUE)

#Se guarda el cálculo de _RF_
bs <- blockSize(segImg)

# se calcula cuantas bandas debe tener la imagen de salida
numBands <- classImage + probImage + threshImage

# Crea el raster de salida y lo guarda
img.out <- brick(segImg, values=FALSE, nl=numBands)
img.out <- writeStart(img.out, outImage, overwrite=TRUE, datatype='INT1U')

# Se Predice la salida
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)
}

# Loop sobre los bloques de los raster de salida desde eCognition y escritura de los nuevos valores clasificados
# Este método permite procesar raster grandes sin problemas de memoria
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)
  # Verificar que los datos nulos se puedan predecir
  is.na(img) <- img == nd
  if (classImage) {
    # Convert the segment ID to the predicted (numeric) class so that a nodata value can be set.
    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 
# Finalmente salva y cierra la conección con la imagen
img.out <- writeStop(img.out)
# Se Guarda el archivo CSV
if (outClassMapCSV != "") {
  write.csv(predValuesDF, file=outClassMapCSV, row.names=FALSE)
}

# Plotea el comportamiento de las variables
varImpPlot(randfor, main ="Comportamiento de las variables")

# Imprime el radio de errror de la matriz de confusión para esta clasificación
confMatrix <- randfor$confusion
cat("\n#################################################################################\n")
## 
## #################################################################################
cat("OOB error  estimado\n", 1 - (sum(diag(confMatrix)) / sum(confMatrix[,1:ncol(confMatrix)-1])), "%\n", sep="")
## OOB error  estimado
## 0.1910112%
cat("#################################################################################\n")
## #################################################################################
cat("Matriz de confusión\n")
## Matriz de confusión
print(randfor$confusion)
##    1  2 3  4  5  6  7  8 class.error
## 1 34  1 1  1  1  5  0  1   0.2272727
## 2  4 17 0  0  0  5  0  0   0.3461538
## 3  1  2 5  0  0  0  0  0   0.3750000
## 4  0  0 0 45  0  0  5  0   0.1000000
## 5  0  0 0  0 14  2  0  0   0.1250000
## 6  1  3 0  2  0 62  3  1   0.1388889
## 7  0  0 0  5  0  2 25  0   0.2187500
## 8  1  0 0  0  0  4  0 14   0.2631579
cat("\nClases de la Matriz de Confusión:\n")
## 
## Clases de la Matriz de Confusión:
cat("\ 1- Viviendas\n 2- Vía Destapada\n 3- Carretera (Asfaltada)\n 4- Árboles\n 5- Clastos (Roca grandes)\n 6- Zona inundada\n 7- Pasto\n 8- Agua\n")
##  1- Viviendas
##  2- Vía Destapada
##  3- Carretera (Asfaltada)
##  4- Árboles
##  5- Clastos (Roca grandes)
##  6- Zona inundada
##  7- Pasto
##  8- Agua
if (outMarginFile != "") {
  # Calcular el excedente
  marginData <- margin(randfor)
  trainingAccuracy <- cbind(marginData[order(marginData)], trainSegs_no_na[order(marginData),])
  
  # Agregar nombres de columnas a la tabal de atributos
  colnames(trainingAccuracy) <- c("margin", "classNum", "cellNum", "segID")
  
  # Calcular cordenadas X, Y para los puntos de los datos de entrenamiento
  xyCoords <- matrix(ncol=2, nrow=0)
  for (z in 1:nrow(trainingAccuracy)) {
    xyCoords <- rbind(xyCoords, xyFromCell(segImg, trainingAccuracy[z,3]))
  }
  
  # Crea un shape file de puntos con la información para que ayude a 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)
}

# Calculara el tiempo de procesamiento
timeDiff <- Sys.time() - startTime
cat("Tiempo de procesamiento:", format(timeDiff), "\n")
## Tiempo de procesamiento: 10.17788 secs
rstout <- raster (outImage)
rstout
## class       : RasterLayer 
## band        : 1  (of  3  bands)
## dimensions  : 3023, 3927, 11871321  (nrow, ncol, ncell)
## resolution  : 0.5, 0.5  (x, y)
## extent      : 1046183, 1048146, 619981, 621492.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.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 : D:\R-Scrips\Mocoa\Forest_NonforestMapTest_v9.tif 
## names       : Forest_NonforestMapTest_v9
colors<-c("#f2f3f4", "#fef9e7", "#aab7b8", "#145a32", "black", "#641e16", "#f4d03f", "#2980B9")
plot(rstout, col=colors, main="Mapa de clasificación")

plotRGB(mrst, r=1, g=2, b=3, scale(10000), stretch="lin", axes=TRUE, main="Imagen Original")

· Edición de los resultados de clasificación

La salida del mapa de clasificación puede tener algunos error, una de las formas de mejorar la clasificación podría ser agregar más capas a la clasificación, otra opción cambiarlas manualmente (corregir la clasificación), Ned Horning preparó un script que nos ayuda a en e esta tarea.

Para esto vamos usar el mapa clasificado, el archivo csv de la clasificación y un nuevo shape en el cual vamos a corregir las clases, para así mejorar la predicción.

  1. Es necesario cargar la liberías necesarias:
library(raster)
library(rgdal)
  1. Se deben configurar las variables para el modelo
setwd("D:/R-Scrips/")
segImage <- "./Mocoa/MeanShift_20_15_100.tif" # Nombre y Ruta del mapa Segmentado
nd <- 1 # valor que van a tomar los datos nulos
# Name and location of the edited classified map
outImagen <- "./Mocoa/Forest_NonforestMapTest_v10.tif" # Nombre y Ruta del mapa de clasicación
inClassMapCSV <- './Mocoa/classMapping_v2.csv' # Nombre y ruta del archivo CSV de clasificación

# Ubicación del shape de puntos que contiene los segmentos que serán midicados 
editPointsDsn <- './Mocoa/reclasificacion.shp'
editPointsLayer <- 'reclasificacion' #Nombre de la capa
newClassNum <- "newClass2" # Nombre o número de las clases que van a ser revisadas
  1. Procesamiento de los datos
## Inicio del Proceso
startTime <- Sys.time()
cat("Hora de Inicio", format(startTime),"\n")
## Hora de Inicio 2018-11-21 16:49:04
cat("Leyendo el Shape\n")
## Leyendo el Shape
vec <- readOGR(editPointsDsn, editPointsLayer) # Lee el Vector
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R-Scrips\Mocoa\reclasificacion.shp", layer: "reclasificacion"
## with 267 features
## It has 5 fields
## Integer64 fields read as strings:  newClass
newClass <- slot(vec, "data")

cat("Leyendo raster de segmentación\n")
## Leyendo raster de segmentación
segImg <- raster(segImage) # Asigna el raster segmentado a una variable

cat("Leyendo archivo CSV\n")
## Leyendo archivo CSV
classMap <- read.csv(inClassMapCSV) # Abre el archivo CSV

# Extrae las IDs de los egmentos que se van a modificar
cat("Extrayendo los puntos que se van a modificar\n")
## Extrayendo los puntos que se van a modificar
changeSegIDs <- cbind(newClass[,newClassNum], extract(segImg, vec))

# Cambia las clases en el mapa de acuerdo al nuevo punto
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]
}
 
# Guarda el mapa de salida
bs <- blockSize(segImg)

# Crea el raster de salida y empieza a guardar
img.out <- raster(segImg)
img.out <- writeStart(img.out, outImagen, overwrite=TRUE, datatype='INT1U')

# Loop sobre los bloques de salida y reclasifica
# Este método permite procesar archivos raster grandes sin problemas de memoria
for (i in 1:bs$n) {
  cat("Procesando bloques", i, "of", bs$n, "\r")
  # library(raster)
    img <- getValues(segImg, row=bs$row[i], nrows=bs$nrows[i])
    # Asigna los valores a los datos nulos
    is.na(img) <- img == nd
    # Convierte los datos sin valor alos datos predichos
  img.match <- as.numeric(classMap$pred[match(img, classMap$segAtr.segnumber)])
  # Asigna los valores a los datos nulos en la imagen de salida
  img.match[is.na(img.match) == TRUE] <- NAvalue(img.out)
  writeValues(img.out, img.match, bs$row[i])
}
## Procesando bloques 1 of 4 
Procesando bloques 2 of 4 
Procesando bloques 3 of 4 
Procesando bloques 4 of 4 
# Guarda y termina la conección a la imagen
img.out <- writeStop(img.out)
  
# Calcula la duración del proceso
timeDiff <- Sys.time() - startTime
cat("Tiempo del Proceso:", format(timeDiff), "\n")
## Tiempo del Proceso: 4.115944 secs
rstout <- raster (outImagen)
colors<-c("#f2f3f4", "#fef9e7", "#aab7b8", "#145a32", "black", "#641e16", "#f4d03f", "#2980B9")
plot(rstout, col=colors, main="Mapa de clasificación")

· Discusión

Para hacer el taller pude probar diferentes configuraciones y diferentes algoritmos de segmentación, decidí hacer el taller con Meanshift por que fue el que mejor se acopló a mi imagen, con los otros o bien creaba muchas particiones o no encontraba algo acorde a la imagen, a parte de eso las posibiliades de configuración a de Meanshift eran muchisimas más que las que ofrecía otro algoritmo como watreshed, esto en un principio no necesariamente failita el proceso de segmentación al haber mas cosas por configurar, pero el proceso es muy intuitivo y permite ir configurando la salida de acuerdo a mi necesidad, un ejemplo es que Ned Horning recomendaba unos valores que al final no servían para mi segmentación y tocó ajustarlos a una agrupación que coincideiera más con mi imagen.

En cuanto a los cálculos estadísticos el algoritmo ofrece una evalucaión del error total y una matriz de error que muestra el error de acuerdo a los puntos de entrenamiento, esto puede ayudarnos a ver donde se equivoca el modelo y tratar de corregir el error por otros medios, en el caso de mi ejercicio por ejemplo muestra en el número de valores que coincide con la predicción y a la vez en a que variable fue predicha.

el mapa de clases coincide de manera bastante acertada con la imagen que se proceso, salvo muy pequeños errores, el error podría deberese a la cantodad de puntos de contro entregados para el modelo, como parte del ejercicio se agregaron diferentes cantidades de puntos de control para ver si hay un patrón entre la cantidad y el error calculado, no se encontró relación en el número de datos tomados con el error calculado, en cambio si hubo una relación entre el tipo de segmento y el error por ejemplo la clase 2 (Caminos destapados) y la clase 3 (Carrertera) obtuvieron el error más alto (35 y 38%), el error podría deberse a la segmentación y al tamaño del área, pero habría que hacer más pruebas para llegar a esta conclusión.

· Conclusiones

Como conclusión es importante destacar la utilidad de este tipo de herramientas a la hora de hacer clasificación, se puede observar que con la clasificación basada en objetose obtiene un mapa más preciso de las zonas.

A pesar de que el error estuvo cercano al 20% el mapa de clases resultado es bastante similar a la imagen.

Se tuvo problemas con la corrigiendo los datos que estaban errado, y el mapa resultado en cambio fue un mapa muy diferente a la imagen, para futuros ejerecicios, es recomendable buscar otros métodos para poder mejorar el error, incluido los recomendados por Ned Horning en su Gúia para clasificar imágenes de satélites usando segmentación y Random Forest, que fue la guía que se usó en este ejercico.

· Bibliografía

Horning N. (2013).Training Guide to Classify Satellite Images using Segmentation and Random Forests
Grolemund G. & Wicham H., 2017 R for Data Science. O’Reilly Media.
R Language Definition, 2018 (Recuperado en línea) disponible en: https://cran.r-project.org/doc/manuals/R-lang.html#Indexing