INTRODUCCION

CONTEXTO

Colombia es un pais con un potencial agricola enorme, cuenta con una gran diversidad de climas, y una gran diversidad de suelos que facilitan la implementacion de gran diversidad de cultivos, como frutales, hortalizas, tuberculos, forestales, adicional a una amplia oferta hidrica. Por su parte el valle del cauca cuenta con psos termicos calidos a paramo, una region montanosa asociada a las cordilleras occidental y central que permite asi una region de valle fisico en medio en donde se ubican principalmente los municipios de Andalucia y Bugalagrande las los cuales basan su economia principalmente en la agricultura y la ganaderia, siendo sus principales cultivos la cana de azucar, limon, maiz, y platano.

PROBLEMA

Por lo cual es fundamental el acompanamiento tecnico continuo con los agricultores, mantener el municipio al tanto de las actividades agricolas predominantes generando un continuo trabajo de compilacion de informacion relacionada con el area agricola que se puede determinar a traves de la cobertura del suelo, el cual debe basarse en la toma de datos in situ, sin embargo con la suficiente informacion se puede plantear el diseno de un modelo que permita obtener datos de cobertura de suelo a traves de fotografias aereas o imagenes satelitales que facilitan la velocidad y la eficiencia de los datos.

OBJETIVOS

Por lo cual se plantea el uso de metodologias de clasificacion supervisada o no supervisada para la determinacion de modelos que permitan obtener rapidamente datos de cobertura del suelo, asi mismo se deben hacer pruebas que permitan reconocer la exactitud de los modelos y lograr identificar que modelo se acerca mas a la realidad y si puede ser utilizado.

REVISION DE CONCEPTOS

cobertura del suelo

La "Cobertura" de la tierra, es la cobertura (bio) fisica que se observa sobre la superficie de la tierra, en un termino amplio no solamente describe la vegetacion y los elementos antropicos existentes sobre la tierra, sino que tambien describen otras superficies terrestres como afloramientos rocosos y cuerpos de agua(IDEAM, 2014).

Uso de la tierra

el uso se relaciona con las actividades humanas o las funciones economicas de una porcion especifica de la Tierra como el uso urbano o industrial, de reserva natural, etc(IDEAM, 2014).

landsat8

es un satelite de observacion terrestre estadounidense lanzado el 11 de febrero de 2013. Es el octavo y mas reciente satelite del proyecto Landsat operado por la NASA y el Servicio Geologico de los Estados Unidos (USGS) desde 1972. El satelite Landsat 8 transporta dos instrumentos OLI y TIRS, que corresponden a las siglas en ingles para Operational Land Imager (OLI) y Thermal Infrared Sensor (TIRS). El sensor OLI provee acceso a nueve bandas espectrales que cubren el espectro desde los 0.433 micras a los 1.390 micras, mientras que TIRS registra de 10.30micras a 12.50micras (USGS. 2013).

##       Banda Nombre         Longitud_de_onda Resolucion
##  [1,] "1"   "Costera"      "0.435 - 0.451"  "30"      
##  [2,] "2"   "Azul"         "0.452 - 0.512"  "30"      
##  [3,] "3"   "Verde"        "0.533 - 0.590"  "30"      
##  [4,] "4"   "Rojo"         "0.636 - 0.673"  "30"      
##  [5,] "5"   "NIR"          "0.851 - 0.879"  "30"      
##  [6,] "6"   "SWIR1"        "1.566 - 1.651"  "30"      
##  [7,] "7"   "SWIR2"        "2.107 - 2.294"  "30"      
##  [8,] "8"   "pancromatica" "0.503 - 0.676"  "15"      
##  [9,] "9"   "cirrus"       "1.363 - 1.384"  "30"      
## [10,] "10"  "TIR1"         "10.60 - 11.19"  "100"     
## [11,] "11"  "TIR2"         "11.50 - 12.51"  "100"

Clasificacion no supervisada

Se basa en el planteamiento de un sistema de clasificacion sin proporcionar ningun dato de respuesta (es decir, no identificamos ningun pixel como perteneciente a una clase en particular). Esto puede parecer extrano, pero puede ser util cuando no tenemos mucho conocimiento previo de un area de estudio. O si tiene un amplio conocimiento de la distribucion de las clases de cobertura terrestre de interes, pero no tiene datos especificos del terreno.

Clasificacion supervisada

tenemos conocimiento previo sobre algunos de los tipos de cobertura del suelo a traves de, por ejemplo, trabajo de campo, datos espaciales de referencia o interpretacion de imagenes de alta resolucion (como las disponibles en los mapas de Google). Se identifican sitios especificos en el area de estudio que representan ejemplos homogeneos de estos tipos de cobertura terrestre conocidos. Estas areas se conocen comunmente como sitios de entrenamiento porque las propiedades espectrales de estos sitios se utilizan para entrenar el algoritmo de clasificacion.

MATERIALES Y METODOS

ZONA DE ESTUDIO

Los datos de referencia fueron tomados en los municipios de Andalucia y Bugalagrande en el departamento del Valle del cauca, Colombia en el ano 2014 a traves del servicio de extension, tomandose datos en formato vector.

Ubicacion

La zona de estudio se encuentra entre los municipios de Andalucia y Bugalagrande, Valle del cauca, Colombia, se delimita a partir de un shapefile entre las coordenadas 4.312815244 a 4.086735824 y -76.288478628 a -76.100159796, donde se identifican fundamentalemente el cultivo de cana de azucar y potreros para ganaderia.

PAIS <- shapefile(input('PAIS','','.shp',T,'shapes'))
DTO <- shapefile(input('DTO','','.shp',T,'shapes'))
MPO <- shapefile(input('MPO','','.shp',T,'shapes'))
ZONA <- shapefile(input('ZONA','','.shp',T,'shapes'))
bt_landsat2014 <- rast(input('LC08_L1TP_009057_20140820_20200911_02_T1_B',1:7,'.tif',T,
                            'landsat8'))
names(bt_landsat2014) <- c('ultrablue','blue','green','red','NIR','SWIR1','SWIR2')
par(mfrow = c(1,2))

plot(ZONA,border='red',axes=T,xlim=c(356513,378235),ylim=c(451125,478191),main='Zona de estudio')

plot(MPO,col='greenyellow',border='White',axes=T,xlim=c(356367,392638),ylim=c(443434,478113))
plot(ZONA,border='red',axes=F,add=T)
text(MPO,MPO$NOMBRE_MPI,cex=0.6)

plot(DTO,col='greenyellow',border=F,axes=T,xlim=c(203614,428384),ylim=c(337273,569428))
plot(MPO,col='green2',border='White',axes=T,add=T)
plot(ZONA,border='red',axes=F,add=T)
text(DTO,DTO$NOMBRE_DPT,cex=0.6)

plot(PAIS,col='darkolivegreen1',border=F,axes=T,xlim=c(-9155,1515773),ylim=c(-488421,1431585))
text(PAIS,PAIS$Pais,cex=0.6)
plot(ZONA,border='red',axes=F,add=T)
plot(DTO,col='greenyellow',border='white',axes=F,add=T)

DATOS UTILIZADOS

Imagenes Landsat

Las imagnes landsat fueron obtenidas a traves de la platafora web https://earthexplorer.usgs.gov/ con el path: 009 y Row: 057 para la fecha del 20 de agosto del 2014 de la coleccion 2 con un nivel de procesamiento 1.

plot(bt_landsat2014$NIR,main='LC08_L1TP_20140820')
plot(ZONA,border='red',add=T)

Poligono de referencia

Los datos describen fundamentalmente 5 coberturas de rio, centro urbano, cultivo de cana de azucar, pasturas y suelo desnudo asi como un poligono que remueve la zona montanosa, dichas coberturas estan agrupadas por lotes que a su vez se agrupan en fincas sin embargo la infomacion no se encuentra completa en relacion al habeas data de la gran mayoria de usuarios.

vallecsv <- read.csv(input('TablaValle','','.csv',T,'shapes'))
vallecsv[10:15,]
shapes <- shapefile(input('valle1','','.shp',T,'Reference'))
shapes <- rbind(shapes[10:97,],shapes[100:102,])
plot(shapes,axes=T,main='Poligonos de referencia')

Herramientas

Se utilizo el programa Rstudio para el analisis de datos, las pruebas de exactitud con las librerias:

  • terra
  • raster
  • sp
  • modeest
  • dplyr
  • rpart

y el complemento Rmarkdown

METODOS

Procesamiento de datos

se genera la funcion "input" para importar datos a partir de 'function'

input <- function(name,num='',type,foldata=F,folder){
  ruta <- getwd()
  if (foldata == T){
    carpeta <- folder
    file <- paste0(ruta,'/',carpeta,'/',name, num, type)
  }else{
    file <- paste0(ruta,'/',name, num, type)
  }
}

posteriormente se realiza un recorte de la imagen satelital con la mascara de la Zona de estudio

MASK <- rast(input('MASK','','.TIF',T,'Reference'))
Extend <- crop(bt_landsat2014, MASK)
Extend<-  -0.1 + Extend*2.0E-5
writeRaster(Extend,filename='landsat2014.tif',overwrite=T)
landsat2014 <- terra::mask(Extend,mask=MASK,filename='Croplandsat2014.tif',overwrite=T)
plot(landsat2014$NIR,main= "recorte zona de estudio")

par algunos analisis es mas conveniente el uso de datos espaciales por capas Stack o brick

br_landsat2014 <- brick(input('landsat2014','','.tif',F))
br_landsat2014
## class      : RasterBrick 
## dimensions : 824, 682, 561968, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 357165, 377625, 452055, 476775  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/LENOVO/Desktop/11/RS/informe1/landsat2014.tif 
## names      : ultrablue,     blue,    green,      red,      NIR,    SWIR1,    SWIR2 
## min values :   0.08276,  0.06332,  0.04064,  0.02462,  0.01154, -0.00640, -0.00310 
## max values :   0.51662,  0.55484,  0.62072,  0.69514,  0.87832,  1.21070,  1.21070

Clasificacion no supervisada

  • Usaremos como capa de referencia el indice NDVI Calculo ndvi
ndvi <- ((br_landsat2014$NIR-br_landsat2014$red)/(br_landsat2014$NIR+br_landsat2014$red))
plot(ndvi,col=rev(terrain.colors(100)),main='landsat2014 NDVI')

los datos del NDVI son obtenidos en formato SpatRaster los cuales son pasados a un formato en una matriz

nr <- values(ndvi)
length(nr)
## [1] 561968

con la matriz contruida se puede dar uso de la clasificacion Kmeans

set.seed(99) #puntos aleatorios
kmncluster <- kmeans(na.omit(nr), #omitir datos nulos
                   centers = 10, #numero de clusters
                   iter.max = 500, #number de iteraciones
                   nstart = 5,
                   algorithm = 'Lloyd' #metodo del algoritmo
                   )
str(kmncluster)
## List of 9
##  $ cluster     : int [1:561968] 6 10 6 6 6 6 6 7 10 10 ...
##  $ centers     : num [1:10, 1] 0.729 0.432 0.494 0.175 -0.192 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:10] "1" "2" "3" "4" ...
##   .. ..$ : NULL
##  $ totss       : num 14244
##  $ withinss    : num [1:10] 30.7 22.7 22.7 46.4 27.2 ...
##  $ tot.withinss: num 266
##  $ betweenss   : num 13978
##  $ size        : int [1:10] 68658 64797 71955 22479 3450 74244 79998 34974 53240 88173
##  $ iter        : int 247
##  $ ifault      : NULL
##  - attr(*, "class")= chr "kmeans"

Retornar nuevamente a datos Spatial

knr <- ndvi
values(knr) <- kmncluster$cluster
knr
## class      : RasterLayer 
## dimensions : 824, 682, 561968  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 357165, 377625, 452055, 476775  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 10  (min, max)

identificacion de clases

mycolor <- c("greenyellow","green","forestgreen","darkblue",
             "dodgerblue4","deepskyblue3","deepskyblue","purple4", "darkorchid1",
             "magenta", "brown1")
plot(knr, main = 'Clasificacion no supervisada', col = mycolor)

Interactuar con datos de referencia

plot(knr, main = 'shapes', col = mycolor)
plot(shapes, lwd=2 ,add=T)

muestreo uniforme de pixeles dentro de los poligonos de referencia

samples <- spsample(shapes, 10000, type='regular')
numero <- length(samples)
df_samples <- as(samples,"SpatialPointsDataFrame")
df_samples@data = data.frame(ID=1:numero, size=1)
plot(knr, main = 'Muestreo uniforme', col = mycolor)
plot(shapes,  add=T)
plot(df_samples, pch = 1, cex = (df_samples$size)/4, add=T)

Extraer valores del raster en los puntos de muestreo

df_samples$estado <- over(df_samples, shapes)$estado
rasdf <- terra::extract(knr, df_samples)

Tabla de Categoria predominante en cada tipo de poligono

modams <- aggregate(rasdf, list(df_samples$estado),mfv)
listms <- aggregate(rasdf, list(df_samples$estado),list)
colnames(modams) <- c("estado","categoria")
colnames(listms) <- c("estado","categoria")
modams

Definir categoria de los poligonos

Lista de datos del clasificador dentro de cada categoria de cobertura

clist <- listms$categoria[[1]]
nlist <- listms$categoria[[2]]
plist <- listms$categoria[[3]]
riolist <- listms$categoria[[4]]
urbanolist <- listms$categoria[[5]]

se plantea una funcion de probabilidad

Identifa <- function(lista,nan=c(),Categoria=Categoria,inicio='c',sens=c(100,10,5/3),diffval=F,hist=T,colhist='black',mapa=F,titulo='mapa de probabilidades', colmap=c('black','bisque4','bisque3','white')) {
  tabla<-aggregate(data.frame(Frecuencia=lista),list(Categoria = lista),length)
  if (length(nan) >= 1){
    for (i in 1:length(nan)) {
      tabla <- rbind(tabla,c(nan[i],0))
    }
  }
  tabla <- rbind(tabla)
  tabla <- tabla[with(tabla, order(tabla$Categoria)),]
  tabla$Categoria <- as.factor(tabla$Categoria)
  dca1<-aov(Frecuencia~Categoria,data=tabla)
  TUK <- TukeyHSD(dca1)
  TUK<-as.list(TUK)
  TUK <- abs(TUK[[1]])
  ListR <- c(inicio)
  if (inicio == 'c'){
    beg <- 1
  }
  if (inicio == 'd'){
    beg <- 1
  }
  if (inicio == 'a'){
    beg <- 0
  }
  if (inicio == 'b'){
    beg <- 0
  }
  if (beg == 1){ 
    for (i in 1:9) {
      a = TUK[i]
      if(a < (max(TUK[1:9])/sens[1])){
        b = 'd'
        ListR<-c(ListR,b)
      }else{
        if(a < (max(TUK[1:9])/sens[2])){
          b = 'c'
          ListR <- c(ListR,b)
        }else {
          if(a < (max(TUK[1:9])/(sens[3]))){
            b = 'b'
            ListR <- c(ListR,b)
          }else {
            b = 'a'
            ListR <- c(ListR,b)
          }
        }
      }
    }
  }else{
    for (i in 1:9) {
      a = TUK[i]
      if(a < (max(TUK[1:9])/sens[1])){
        b = 'a'
        ListR<-c(ListR,b)
      }else{
        if(a < (max(TUK[1:9])/sens[2])){
          b = 'b'
          ListR <- c(ListR,b)
        }else {
          if(a < (max(TUK[1:9])/(sens[3]))){
            b = 'c'
            ListR <- c(ListR,b)
          }else {
            b = 'd'
            ListR <- c(ListR,b)
          }
        }
      }
    }
  }
  
  if(diffval == T){
    print(TUK[1:length(unique(lista))])
  }
  if(hist == T){
    tabla1 <- tabla
    tabla1$Diff <- ListR
    ymax <- max(TUK[1:9]) + 100
    plot.default(tabla,type = "h", ylim = c(0,ymax), col=colhist,lwd=4)
    axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
    text(tabla1,tabla1$Diff,pos=3)
  }
  if(mapa == T){
  atabla <- filter(tabla1,Diff == "a")
  btabla <- filter(tabla1,Diff == "b")
  ctabla <- filter(tabla1,Diff == "c")
  atablavalue <- atabla[,1]
  btablavalue <- btabla[,1]
  ctablavalue <- ctabla[,1]
  colormap<- c(rep(colmap[4],10))
  colormap[atablavalue] <- colmap[1]
  colormap[btablavalue] <- colmap[2]
  colormap[ctablavalue] <- colmap[3]
  plot(knr, main = titulo , col = colormap)
  legend(x ='bottomright',legend=c('Alta','Media','Baja','nula'),fill=colmap)
  }
}

Una vez obtenida la funcion de probabilidad se ejecuta en cada una de las categorias de cobertura para identificar las zonas de mayor probabilidad de cada categoria

  • Identificacion de cuerpos de agua
Identifa(riolist,colhist = 'blue',mapa=T,titulo = 'Probabilidad cuerpos de agua',
         colmap=c('blue','dodgerblue','lightblue1','white'))

  • Identificacion de centros urbanos
Identifa(urbanolist,colhist = 'gray',mapa = T,titulo = 'Probabilidad centros urbanos',
         colmap=c('azure4','azure3','azure2','white'))

  • identificacion de suelo desnudo
Identifa(nlist,nan=c(1,5),inicio='d',sens = c(100,3,1.5),colhist = 'brown', mapa = T, titulo = 'Probabilidad suelos desnudos',colmap = c('brown4','brown2','darksalmon','white'))

  • identificacion de pastos
Identifa(plist,nan=c(4,5,8),inicio='c',sens = c(10,3,1.5),colhist = 'chartreuse4', mapa = T, titulo = 'Probabilidad pasturas',colmap = c('chartreuse4','chartreuse2','greenyellow','white'))

  • identificacion de cultivo de cana de azucar
Identifa(clist,inicio='a',sens = c(5,1.5,1.05),colhist = 'chartreuse4', mapa = T, titulo = 'Probabilidad Cultivo',colmap = c('darkgreen','green3','green','white'))

Evaluacion del modelo no supervisado

samples$estado <- over(samples,shapes)$estado
xy <- as.matrix(geom(samples)[,c('x','y')])
df <- extract(knr, xy)
sampdata <- data.frame(estado = samples$estado, df)
cartmodel <- rpart(as.factor(estado)~., data = sampdata, method = 'class',
minsplit = 5)

classun <- c("cutivo","suelo","pastos","rio","urbano")
mycolorun <- c('darkgreen', 'Brown', 'chartreuse', 'blue', 'gray')
classdfun <- data.frame(classvalue = c(1,2,3,4,5), classnames = classun, color = mycolorun)

set.seed(99)
kun <- 5
jun <- sample(rep(1:kun, each = round(nrow(sampdata))/kun))
xun <- list()
for (kun in 1:5) {
    train <- sampdata[jun!= kun, ]
    test <- sampdata[jun == kun, ]
    cart <- rpart(as.factor(estado)~., data=train, method = 'class',
                  minsplit = 5)
    pclass <- predict(cart, test, na.rm = TRUE)
    pclass <- apply(pclass, 1, which.max)
    xun[[kun]] <- cbind(test$estado, as.integer(pclass))
}
yun <- do.call(rbind, xun)
yun <- data.frame(yun)
colnames(yun) <- c('observed', 'predicted')
# confusion matrix
conmatun <- table(yun)
# change the name of the classes
colnames(conmatun) <- classdfun$classnames
rownames(conmatun) <- classdfun$classnames
nun <- sum(conmatun)
diagun <- diag(conmatun)
OAun <- sum(diagun) / nun
## Kappa
rowsumsun <- apply(conmatun, 1, sum)
pun <- rowsumsun / nun
colsumsun <- apply(conmatun, 2, sum)
qun <- colsumsun / nun
expAccuracyun <- sum(pun*qun)
kappaun <- (OAun - expAccuracyun) / (1 - expAccuracyun)
PAun <- diagun / colsumsun
UAun <- diagun / rowsumsun
outAccun <- data.frame(producerAccuracy = PAun, userAccuracy = UAun)

clasificacion supervisada

Se inicia con un muestreo sobre los poligonos de referencia como base del modelo

set.seed(1)
samples1 <- spsample(shapes, 10000, type='random')
samples1$estado <- over(samples1,shapes)$estado

extraer valores espectrales

xy <- as.matrix(geom(samples1)[,c('x','y')])
df <- extract(Extend, xy)[,-1]
head(df)
sampdata1 <- data.frame(estado = samples1$estado, df)

entrenar al clasificador

cartmodel <- rpart(as.factor(estado)~., data = sampdata1, method = 'class',
minsplit = 5)
print(cartmodel)
## n= 10000 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 10000 6666 c (0.33 0.23 0.15 0.072 0.21)  
##    2) SWIR2< 0.07355 4872 1771 c (0.64 0.045 0.15 0.15 0.024)  
##      4) SWIR1>=0.08671 4188 1097 c (0.74 0.052 0.17 0.011 0.026)  
##        8) green< 0.08579 3187  405 c (0.87 0.055 0.035 0.013 0.024)  
##         16) NIR>=0.24511 2954  216 c (0.93 0.01 0.038 0.0064 0.019) *
##         17) NIR< 0.24511 233   89 n (0.19 0.62 0 0.099 0.094) *
##        9) green>=0.08579 1001  387 p (0.31 0.044 0.61 0.004 0.03)  
##         18) blue>=0.09315 196   77 c (0.61 0.2 0.02 0.02 0.15) *
##         19) blue< 0.09315 805  195 p (0.24 0.0062 0.76 0 0) *
##      5) SWIR1< 0.08671 684   19 rio (0.015 0 0 0.97 0.013) *
##    3) SWIR2>=0.07355 5128 3034 n (0.045 0.41 0.15 0.0014 0.39)  
##      6) NIR< 0.22735 3206 1216 n (0.019 0.62 0.0075 0.0022 0.35)  
##       12) SWIR2< 0.14343 2200  413 n (0.028 0.81 0.011 0.0032 0.15) *
##       13) SWIR2>=0.14343 1006  203 urbano (0 0.2 0 0 0.8)  
##         26) red>=0.11998 283  113 n (0 0.6 0 0 0.4)  
##           52) SWIR2< 0.1777 179   20 n (0 0.89 0 0 0.11) *
##           53) SWIR2>=0.1777 104   11 urbano (0 0.11 0 0 0.89) *
##         27) red< 0.11998 723   33 urbano (0 0.046 0 0 0.95) *
##      7) NIR>=0.22735 1922 1035 urbano (0.089 0.054 0.39 0 0.46)  
##       14) SWIR2< 0.10427 1176  435 p (0.15 0.048 0.63 0 0.18)  
##         28) green>=0.08846 890  172 p (0.12 0.048 0.81 0 0.029) *
##         29) green< 0.08846 286  105 urbano (0.24 0.045 0.08 0 0.63) *
##       15) SWIR2>=0.10427 746   66 urbano (0 0.064 0.024 0 0.91) *

Arbol de clasificacion

a partir de los datos de referencia la herramienta rpart permite generar un arbol de clasificacion que tiene como respuesta final las 5 coberturas iniciales

plot(cartmodel, uniform=TRUE, main="Arbol de clasificacion")
text(cartmodel, cex = 0.6)

#### clasificar Una vez generado el modelo es aplicarlo sobre la imagen landsat de tal manera que identifique la cobertura esperada para cada pixel para este caso el arbol de clasificacion uso las bandas azul, verde, rojo, NIR, SWIR1 y SWIR2

classified <- predict(Extend, cartmodel, na.rm = TRUE)
plot(classified)

class <- c("cutivo","suelo","pastos","rio","urbano")
mycolor <- c('darkgreen', 'Brown', 'chartreuse', 'blue', 'gray')
classdf <- data.frame(classvalue = c(1,2,3,4,5),
                      classnames = class,
                      color = mycolor)
lulc <- app(classified, fun = which.max)
lulcc <- as.factor(lulc)
levels(lulcc) <- as.character(classdf$classnames)

evaluacion del modelo supervisado

para las pruebas de exactitud cuando no se cuenta con la posibilidad de ir a campo o aun es inviable hacerlo se puede ensayar el modelo con muestreos aleatorios en este caso se usaran de 5 muestras

set.seed(99)
k <- 5
j <- sample(rep(1:k, each = round(nrow(sampdata1))/k))
table(j)
## j
##    1    2    3    4    5 
## 2000 2000 2000 2000 2000

prueba del modelo

x <- list()
for (k in 1:5) {
    train <- sampdata1[j!= k, ]
    test <- sampdata1[j == k, ]
    cart <- rpart(as.factor(estado)~., data=train, method = 'class',
                  minsplit = 5)
    pclass <- predict(cart, test, na.rm = TRUE)
    pclass <- apply(pclass, 1, which.max)
    x[[k]] <- cbind(test$estado, as.integer(pclass))
}

Convinar los 5 marco de datos

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(y)
# change the name of the classes
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
print(conmat)
##         predicted
## observed cutivo suelo pastos  rio urbano
##   cutivo   2841   132    310   10     41
##   suelo      85  2005     49    3    170
##   pastos    111    32   1311    0     54
##   rio        22    35      0  661      0
##   urbano    103   406     41    8   1570

exactitud general y "kappa"

n <- sum(conmat)
##exactitud general
diag <- diag(conmat)
OA <- sum(diag) / n
## Kappa
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)

exactitud de usuario

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

RESULTADOS

Finalmente se obtienen dos mapas de cobertura del suelo basados en una clasificacion no supervisada y supervisada correspondientemente, a su vez a cada modelos se le realizaron evaluaciones de exactitud con la exactitud general, el Kappa, exactitud de productor y de usuario.

mapa de clasificacion no supervisada

newcolor <- c('darkgreen','white','greenyellow','brown','blue','greenyellow','greenyellow','brown','brown','chartreuse4')
plot(knr, main = 'Mapa de cobertura del modelo no supervisado',col = newcolor)
plot(shapes[1:87,],axes=T, border='black',add=T)
plot(ZONA,border='red',lwd=2,add=T)
legend(x ='bottomright',legend=c('Agua','Suelo desnudo','pastura','Cultivo joven','Cultivo','Sin definir'), fill=c('blue','brown','greenyellow','chartreuse4','darkgreen','white'))

como se observa en la figura 'Mapa de cobertura del modelo no supervisado' finalmente teniendo en cuenta la maxima probabilidad para cada valor del modelo predictivo, ocupando asi cada categoria de cobertura diferentes valores del modelo, donde se logra identificar claramente el rio como un cuerpo de agua bien diferenciado, asi como algunos lotes de color verde oscuro representando el cultivo de cana, sin embargo el 10 valor del modelo parecia estar entre pastura y el cultivo por lo que se asume que se puede tratar de un lote en barbecho, cana de azuca con poca biomasa u otra cobertura no reconocida en el muestreo. Asi mismo se puede observar la gran cobertura que presenta la categoria pastos que con analisis empirico de la imagen satelital parecen superficies escarpadas en la zona montanosa y otros asociados a cultivos con mayor y menor biomasa por lo que parece ser que el modelo es poco exacto al reconocer pasturas; en ese sentido la cobertura de ciudad no es mostrada en el mapa debido a su clara ausensia de criterio y alta dispariedad en su ubicacion. Finalmente cabe resaltar que el 2do valor del modelo no se ajustaba lo suficiente a ninguna categoria de cobertura por lo cual corresponde a una cobertura sin definir mostrada en el mapa de color blanco.

exactitud del modelo no supervisado

## [1] "n= 10016"
## [1] "exactitud general 0.65055910543131"
## [1] "Kappa = 0.543099746554343"

Se tomaron 2.000 puntos de muestreo en 5 muestreos diferentes para la evaluacion del modelo permitiendo obtener una exactitud general realmente baja en el regimen del 65%, donde por lo regular en los trabajos de campo y de mas se tienen umbrales de inexactitud hasta del 10% lo cual es bastante alejado del modelo, en ese sentido el Kappa con una mayor rigurosidad determina asi un 54% de exactitud claramente. Una vez observada la tabla de exactitud de productor y usuario se puede identificar cuales categorias se ajustan mejor al modelos, por ejemplo los cuerpos de agua mustran la mayor exactitud de productor con un 98% mientras que para la de usuario el que mejor se ajusta es para identificacion de cana de azucar, por el cotrario la categoria mas baja en ambos criterios es urbano lo que coincide con el histograma realizado en la determinacion de centros urbano.

mapa de clasificacion supervisada

plot(lulcc,col=mycolor,main='cobertura clasificacion supervisada')
plot(shapes[1:87,],axes=T, border='black',add=T)
plot(ZONA,border='red',lwd=2,add=T)

A su vez la clasificacion supervisada genera un mapa de cobertura donde nuevamente es facil identificar afluentes, sin embargo parece notarse una gran confusion entre la superficie de suelo desnudo o de baja vegetacion con la categoria urbano, mientras se identifica claramente los lotes con altaproduccion de biomasa de color verde oscuro correspondiente a cana de azucar. Si se asume la categoria urbano como suelo desnudo tambien se puede reconocer claramente que ambos sistemas montanosos son pobres en vegetacion.

exactitud del modelo supervisado

## [1] "n= 10000"
## [1] "exactitud general 0.8388"
## [1] "Kappa = 0.789139998326724"

por su lado las pruebas de exactitud para el caso de clasificacion supervisada tambien usaron 2.000 puntos aleatorios con 5 muestreos donde vemos una exactitud general del 83% con un Kappa del 79%, lo que aun es alejado de un modelo con criterio suficiente para ser utilizado con confianza, sin embargo, observando la exactitud de productor y de usuario ya existe suficiente confienza para usar el modelo si quiera en cuerpos de agua ya que son del 97% y 92% correspondientemente, mientras que para cultivo es del 89% y 85% lo que tiende a acercarce a lo esperado, sin embargo el dato de urbano de 85% y 73% deja una incertidumbre alta con relacion al modelo ya que parecen ser exactitudes muy altas para la interpretacion del mapa de cobertura observado.

DISCUCION

Modelo de clasificacion no supervisado

Una problematica del modelo no supervisado es la dificultad para ser procesado, la poca reproducibilidad y aplicabilidad, ya que debe ser reevaluado en cada mapa ya sea por una investigacion de diferentes zonas o diferentes tiempos ya que los valores arojados no se logra identificar claramente el criterio de clasificacion. Sin embargo en algunos casos permite identificar diferencias dentro de las categorias preestablecidas donde si ejecutaramos el modelo con datos de referencia los datos de dos diferentes astados fenologicas se verian combinados por el modelo supervisado mientras que con este puede lograr diferenciarse y apoyarse con datos de campo para mejorar el nivel de detalle de los datos de cobertura, asi mismo resulta ser un buen metodo para reconocer un territorio con poco conocimiento. lo que corresponde a este modelo se puede identificar que la categorias de rio, urbano, cultivo, pastura y suelo desnudo parecen insuficientes para describir la diversidad de datos mostrados en el mapa de prediciones, y por lo cual incluso la categoria 2 queda completamente vacia.Tambien cabe resaltar que el modelo fue disenado sobre una capa de indice NDVI por lo cual tambien se puede plantear la busqueda de diferentes capas que faciliten la identificacion de coberturas en el territorio.

Modelo de clasificacion supervisado

Al ser un modelo dependiente de datos de campo puede presentar un alto margen de error por la accion humana y que directamente afectara la exactitud general del modelo, por lo cual un tratamiento de la infomacion adecuado permitira tener un modelo adecuado, ya sea con la ayuda de normalizar datos, identificacion de datos atipicos etc. Por otro lado un modelo con la suficiente confianza puede ser utilizado en analisis temporales o de diferentes zonas ya que muestra claramente a traves de un arbol de clasificacion (que para este caso utilizo las 7 bandas suministradas) los criterios que se utilizan facilitando su reproducibilidad. Un problema que presenta este modelo es el reconocimiento de los espacios urbanos como suelo desnudo, esto sucede con facilidad ya que ambas categorias presentan altas reflectancias por el materia en relacion a los minerales del suelo en contraste, podria ajustarce mas la sensibilidad de los cascos urbanos permitiendo descartar zonas verdes dentro del centro urbano por un limite difuso con la seccion rural.

comparacion de modelos

Ambos modelos presentan una gran correspondencia con la categoria de cuerpos de agua aunque basado en las diferentes evaluaciones de exactitud el modelo supervisado es claramente superior, con la ventaja de que puede ser replicado, por lo cual se pueden difernciar en el objetivo a cumplir con cada clasificacion donde la clasificacion no supervisada facilita la genereacion de informacion, puede definir los niveles de detalle de cobertura, puede definir caracteristicas adicionales, mientras que el modelo supervisado se acerca mas al cumplimiento del objetivo de generar una herramienta util para la compilacion de datos de cobertura, a traves del tiempo, sin embargo ambos modelos aun carecen de confianza suficiente, se recomienda tomar mayor cantidad de datos en campo, tambien se puede hacer una busqueda de valores atipicos, pensar en modelos diferentes a kmeans.

CONCLUCIONES

A partir de los datos obtenidos es un acercamiento a la comprension de modelos de clasificacion donde el modelo de clasificacion supervisado se ajusta mejor a los datos de referencia aunque aun no cumple con los criterios de confienza por lo que se debe contruir una matriz de datos mas solida que garantice un mejor input para el modelo, asi como la busqueda de modelos mas riguros basados en inteligencia artificial como analisis neurales o de machine learning.

REFERENCIAS