Bases para una Planificación
Pablo Cesar Serrano Agudelo
14 de julio de 2019
###Introduccion
Los ecosistemas de páramo se distinguen por su provisión de servicios hidrológicos, la protección de la biodiversidad y el almacenamiento de carbono (Martínez et al, 2009), destacar estas funciones del ecosistema de páramo en la política internacional de cambio climático, contribuiría a la adopción de incentivos económicos para apoyar un mejor manejo de los recursos naturales en estos ecosistemas (Ward et al, 2015), esto debido a que por la variabilidad climática y el aumento de la actividad antrópica, en especial por el uso de la tierra y el cambio de coberturas; acciones antrópicas tales como la deforestación, forestación, el avance de la frontera agrícola, el pastoreo en áreas inadecuadas o el sobrepastoreo, afectan estos servicios ecosistémicos, razón por la cual se ha incentivado la implantación de medidas de protección y conservación de los páramos resaltando la importancia de estos ecosistemas (Quichimbo et al, 2012).
Una de las bases para poder elaborar planes de manejo para la protección y conservación de estas zonas es la clasificación de coberturas, ya que mediante esta se pueden observar zonas de uso inadecuado de la tierra, o mediante estudios multitemporales desarrollar mapas de cambios de uso de la tierra para identificar zonas de mayor presión por cambios antrópicos.
Las condiciones geográficas de los ecosistemas de páramo no permiten realizar una clasificación de estas coberturas con ideales datos de campo, debido al acceso algunas zonas, también a las grandes extensiones que ocupan estos paramos debido a esto la metodología utilizada y más apropiada para estimar área de cobertura y usos del suelo se basa en equipos de teledetección y percepción remota (IGAC, 1979).
Para poder solucionar este problema se plantean realizar una clasificación no supervisada y supervisada en el páramo de Rabanal con el fin de evaluar si a estas clasificaciones sirven como base para formular estudios de planificacion territorial enfocados a la conservación y protección del ecosistema de páramo. ###Datos y métodos
Se utilizó como dato la imagen Landsat 8 del path y row 8_56, descargada de la página de internet USGSS, con fecha de toma del 2019-01-30, y donde se utilizaron de las bandas 2 a la 7 correspondientes a las bandas ‘blue’ ,‘green’ ,‘red’ ,‘nir’ ,‘swir1’ ,‘swir2’ ,‘cirrus’ ,‘tirs1’ , ‘tirs2’ con una resolución espacial de 30 metros.
####Área de estudio
El área de estudio será el páramo Rabanal delimitado a escala 1:100.000 por el instituto Von Humboldt y el Misterio de Medio Ambiente en la Resolución 1768 del 28 octubre de 2016 y que se encuentre dentro de la jurisdicción del Departamento de Boyacá; esta zona se ubica en la parte central de la cordillera oriental entre las coordenadas geográficas 5°28’23.93“N, 73°38’0.52”W y 5°20’8.20“N, 73°31’8.84”W (X=1049205, Y= 1097011 y X=1061892, Y= 1081794 magna Colombia-Bogotá). La zona se encuentra dentro del municipio de Ráquira, del municipio de Samacá y del municipio de Ventaquemada, y al occidente colinda con los municipios de Guachetá y Lenguazaque del departamento de Cundinamarca para un área total de 7249.6 hectáreas.
Área de estudio
Metodologia
La metodologia utilizada se basa en los siguientes procesos: • Correción radiometrica. • Análisis de bandas corregidas. • Calculo de NDVI • Recorte de la imagen al área de estudio. • Clasificación no supervisada y supervisada. • Matriz de confusion.
Metodologia
###Resultados
Como primera medida se muestra toda la escena de la imagen satelital Landsat la cual como se puede ver gran parte está cubierta de nubes pero la zona de estudio está libre de nubes.
En esta parte se muestra la escena total de la imagen Landsat.
MTL_Landsat <- "LC08_L1TP_008056_20181230_20190130_01_T1_MTL.txt"
MTL_Landsat## [1] "LC08_L1TP_008056_20181230_20190130_01_T1_MTL.txt"metadata <- readMeta(MTL_Landsat)
metadata## $METADATA_FILE
## [1] "LC08_L1TP_008056_20181230_20190130_01_T1_MTL.txt"
##
## $METADATA_FORMAT
## [1] "MTL"
##
## $SATELLITE
## [1] "LANDSAT8"
##
## $SENSOR
## [1] "OLI_TIRS"
##
## $SCENE_ID
## [1] "LC80080562018364LGN00"
##
## $ACQUISITION_DATE
## [1] "2018-12-30 15:12:13 GMT"
##
## $PROCESSING_DATE
## [1] "2019-01-30 GMT"
##
## $PATH_ROW
## path row
## 8 56
##
## $PROJECTION
## CRS arguments:
## +proj=utm +zone=18 +units=m +datum=WGS84 +ellps=WGS84
## +towgs84=0,0,0
##
## $SOLAR_PARAMETERS
## azimuth elevation distance
## 138.8810140 51.0324022 0.9833282
##
## $DATA
## FILES BANDS QUANTITY
## B1_dn LC08_L1TP_008056_20181230_20190130_01_T1_B1.TIF B1_dn dn
## B2_dn LC08_L1TP_008056_20181230_20190130_01_T1_B2.TIF B2_dn dn
## B3_dn LC08_L1TP_008056_20181230_20190130_01_T1_B3.TIF B3_dn dn
## B4_dn LC08_L1TP_008056_20181230_20190130_01_T1_B4.TIF B4_dn dn
## B5_dn LC08_L1TP_008056_20181230_20190130_01_T1_B5.TIF B5_dn dn
## B6_dn LC08_L1TP_008056_20181230_20190130_01_T1_B6.TIF B6_dn dn
## B7_dn LC08_L1TP_008056_20181230_20190130_01_T1_B7.TIF B7_dn dn
## B9_dn LC08_L1TP_008056_20181230_20190130_01_T1_B9.TIF B9_dn dn
## B10_dn LC08_L1TP_008056_20181230_20190130_01_T1_B10.TIF B10_dn dn
## B11_dn LC08_L1TP_008056_20181230_20190130_01_T1_B11.TIF B11_dn dn
## B8_dn LC08_L1TP_008056_20181230_20190130_01_T1_B8.TIF B8_dn dn
## QA_dn LC08_L1TP_008056_20181230_20190130_01_T1_BQA.TIF QA_dn dn
## CATEGORY NA_VALUE SATURATE_VALUE SCALE_FACTOR DATA_TYPE
## B1_dn image NA NA 1 NA
## B2_dn image NA NA 1 NA
## B3_dn image NA NA 1 NA
## B4_dn image NA NA 1 NA
## B5_dn image NA NA 1 NA
## B6_dn image NA NA 1 NA
## B7_dn image NA NA 1 NA
## B9_dn image NA NA 1 NA
## B10_dn image NA NA 1 NA
## B11_dn image NA NA 1 NA
## B8_dn pan NA NA 1 NA
## QA_dn qa NA NA 1 NA
## SPATIAL_RESOLUTION RADIOMETRIC_RESOLUTION
## B1_dn 30 16
## B2_dn 30 16
## B3_dn 30 16
## B4_dn 30 16
## B5_dn 30 16
## B6_dn 30 16
## B7_dn 30 16
## B9_dn 30 16
## B10_dn 30 16
## B11_dn 30 16
## B8_dn 15 16
## QA_dn 30 16
##
## $CALRAD
## offset gain
## B1_dn -64.92558 0.0129850
## B2_dn -66.48459 0.0132970
## B3_dn -61.26498 0.0122530
## B4_dn -51.66208 0.0103320
## B5_dn -31.61462 0.0063229
## B6_dn -7.86227 0.0015725
## B7_dn -2.65001 0.0005300
## B8_dn -58.46726 0.0116930
## B9_dn -12.35571 0.0024711
## B10_dn 0.10000 0.0003342
## B11_dn 0.10000 0.0003342
##
## $CALREF
## offset gain
## B1_dn -0.1 2e-05
## B2_dn -0.1 2e-05
## B3_dn -0.1 2e-05
## B4_dn -0.1 2e-05
## B5_dn -0.1 2e-05
## B6_dn -0.1 2e-05
## B7_dn -0.1 2e-05
## B8_dn -0.1 2e-05
## B9_dn -0.1 2e-05
##
## $CALBT
## K1 K2
## B10_dn 774.8853 480.8883
## B11_dn 1321.0789 1201.1442
##
## attr(,"class")
## [1] "ImageMetaData" "RStoolbox"lsat <- stackMeta(metadata)
lsat## class : RasterStack
## dimensions : 7731, 7571, 58531401, 10 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 480885, 708015, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B1_dn, B2_dn, B3_dn, B4_dn, B5_dn, B6_dn, B7_dn, B9_dn, B10_dn, B11_dn
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535plotRGB(lsat, main = "Landsat 8 composición RGB sin correcion",
r = 4, g = 3, b = 2, scale = 65535, stretch = "lin", axes = TRUE)
El primer paso de la metodologia es hacer la correción
hazeDN <- estimateHaze(lsat, hazeBands = 1:7, darkProp = 0.01,
plot = TRUE)
lsat_sref_sdos <- radCor(lsat, metaData = metadata, method = "sdos",
hazeValues = hazeDN, hazeBands = 1:7)
lsat_sref_sdos## class : RasterStack
## dimensions : 7731, 7571, 58531401, 10 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 480885, 708015, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre, B9_sre, B10_bt, B11_bt
## min values : 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 53.69796, 126.58460
## max values : 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 0.2267014, 114.0384445, 246.5825917names(lsat_sref_sdos) <- c("costalaerosol", "blue", "green",
"red", "NIR", "SWIR1", "SWIR2", "cirrus", "tirs1", "tirs2")
banda1 = lsat_sref_sdos$costalaerosol
banda2 = lsat_sref_sdos$blue
banda3 = lsat_sref_sdos$green
banda4 = lsat_sref_sdos$red
banda5 = lsat_sref_sdos$NIR
banda6 = lsat_sref_sdos$SWIR1
banda7 = lsat_sref_sdos$SWIR2Se realiza la relación entre las bandas para poder analizarlas y se crea el stack de las bandas para poder realizar el cálculo NDVI y recortar la imagen al área de estudio.
pairs(lsat_sref_sdos[[2:7]])
multi_banda_30m = stack(banda1, banda2, banda3, banda4, banda5,
banda6, banda7)
multi_banda_30m## class : RasterStack
## dimensions : 7731, 7571, 58531401, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 480885, 708015, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : costalaerosol, blue, green, red, NIR, SWIR1, SWIR2
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1plotRGB(multi_banda_30m, main = "Landsat 8 composición RGB corregida",
r = 4, g = 3, b = 2, scale = 65535, stretch = "lin", axes = TRUE)
lsat_sref_sdos## class : RasterStack
## dimensions : 7731, 7571, 58531401, 10 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 480885, 708015, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : costalaerosol, blue, green, red, NIR, SWIR1, SWIR2, cirrus, tirs1, tirs2
## min values : 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 53.69796, 126.58460
## max values : 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 0.2267014, 114.0384445, 246.5825917areaestudio_lim <- shapefile("rabanal")
areaes <- spTransform(areaestudio_lim, crs(lsat_sref_sdos))
areaes## class : SpatialPolygonsDataFrame
## features : 1
## extent : 651566.5, 663827.8, 590058.6, 604946.4 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 14
## names : OBJECTID_1, OBJECTID, CODSECTOR, CODIGO, COMPLEJO, DISTRITO, SECTOR, CODCOMPLEJ, CODID, IdNum, AREA, Shape_Leng, Shape_Le_1, Shape_Area
## value : 3, 55, CE, CE-CM-RRB, Rabanal y rÃo Bogotá, Cundinamarca, Cordillera Oriental, RRB, CE-CM-RRB-01, 54, 11096.5483118, 117047.828732, 82456.0101559, 7249.60368708landsatcrop <- crop(multi_banda_30m, areaestudio_lim)
landsat_clip = mask(landsatcrop, areaestudio_lim)
plot(landsat_clip)
Al obtener el grafico que correlaciona las bandas se encuentra alta correlación entre las bandas 2,3,4 que pertenecen al espectro visible azul, verde, rojo, así mismo otra grupo de bandas que se encuentra alta correlación entre la bandas del NIR y SWIR1, SWIR1 y SWIR 2. Es debido a este análisis de correlación de las bandas y las características de la zona de estudio; en donde predomina la vegetación, que se realiza el Índice de Vegetación de Diferencia Normalizada NDVIya que la banda roja y la banda infrarrojo cercano representan al grupo de bandas de interés (por tener la misma resolución espacial).
iv <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
iv <- (bk - bi)/(bk + bi)
return(iv)
}
ndvi_landsat <- iv(lsat_sref_sdos, 5, 4)
ndvi_landsat## class : RasterLayer
## dimensions : 7731, 7571, 58531401 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 480885, 708015, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/USUARIO/AppData/Local/Temp/Rtmp0ovGjN/raster/r_tmp_2019-10-14_225231_7068_96685.grd
## names : layer
## values : -1, 1 (min, max)plot(ndvi_landsat, main = "Landsat-NDVI")
areaes## class : SpatialPolygonsDataFrame
## features : 1
## extent : 651566.5, 663827.8, 590058.6, 604946.4 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 14
## names : OBJECTID_1, OBJECTID, CODSECTOR, CODIGO, COMPLEJO, DISTRITO, SECTOR, CODCOMPLEJ, CODID, IdNum, AREA, Shape_Leng, Shape_Le_1, Shape_Area
## value : 3, 55, CE, CE-CM-RRB, Rabanal y rÃo Bogotá, Cundinamarca, Cordillera Oriental, RRB, CE-CM-RRB-01, 54, 11096.5483118, 117047.828732, 82456.0101559, 7249.60368708ndvicrop <- crop(ndvi_landsat, areaes)
ndvimask = mask(ndvicrop, areaestudio_lim)
plot(ndvimask)
Es a partir de este NDVI que se realizó la clasificación no supervisada (por el método de Loyd), para la zona de estudio, en donde se generaron 10 clases generales de cobertura de tierra tanto para la clasificación no supervisada como para la supervisada, en la siguiente parte se muestra la no supervisada.
nr <- getValues(ndvicrop)
str(nr)## num [1:202864] 0.46 0.574 0.643 0.634 0.727 ...set.seed(99)
kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 200,
nstart = 10, algorithm = "Lloyd")
str(kmncluster)## List of 9
## $ cluster : int [1:202864] 5 1 3 3 4 7 4 1 9 9 ...
## $ centers : num [1:10, 1] 0.564 -0.938 0.637 0.705 0.493 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 11006
## $ withinss : num [1:10] 12.8 39.4 12.9 12.9 13.2 ...
## $ tot.withinss: num 174
## $ betweenss : num 10832
## $ size : int [1:10] 30044 2697 31177 34387 30195 18810 34858 4273 15648 775
## $ iter : int 173
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"knr <- raster(ndvimask)
knr## class : RasterLayer
## dimensions : 496, 409, 202864 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 651555, 663825, 590055, 604935 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0values(knr) <- kmncluster$cluster
knr## class : RasterLayer
## dimensions : 496, 409, 202864 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 651555, 663825, 590055, 604935 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 1, 10 (min, max)# plot(knr)
nosupervisada = mask(knr, areaestudio_lim)
classcolor <- c("lightgreen", "deepskyblue", "darkolivegreen1",
"darkolivegreen", "mediumaquamarine", "dimgray", "firebrick4",
"green1", "lawngreen", "khaki2")
# par(mfrow = c(1,2)) plot(ndvimask, col=
# rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(nosupervisada, col = classcolor)
Clase no supervisada Agricultura 4 Agua 2 Arbustal 7 Bosque 3 Herbazal 1 Humedal 3 Minería 10 Pastos 8 Pino 6 Suelo 5
Es importante resaltar que por ejemplo para la clase suelo se tuvieron en cuento tanto suelos degradados como suelos desnudos y suelos arados con fines productivos agrícolamente hablando, se podría indicar que en el caso de las represas están siendo bien representadas por lo menos visualmente, estas son de gran importancia para la región ya que la represa de Teatinos es la que abastece los acueductos de Tunja entre otros municipios, datos que se comparan cuando se presente la matriz de confusión.
Para realizar la clasificación supervisada se utilizaron datos de entrenamientos recogidos en pantalla, en donde por el método “class” que utiliza arboles de decisión se generaron 10 clases al igual que en la clasificación no supervisada.
#### supervisada
library(rasterVis)
sample <- shapefile("puntosentrena")
sampvals <- extract(landsat_clip, sample, df = TRUE)
sampvals <- sampvals[, -1]
sampdata <- data.frame(classvalue = sample$COBER, sampvals)
library(rpart)
cart <- rpart(as.factor(classvalue) ~ ., data = sampdata, method = "class",
minsplit = 4)
plot(cart, uniform = TRUE, main = "Classification Tree")
text(cart, cex = 0.5)
pr2011 <- predict(landsat_clip, cart, type = "class")
plot(pr2011, col = classcolor)
pr2011 <- ratify(pr2011)
rat <- levels(pr2011)[[1]]
levels(pr2011) <- rat
nlcdclass <- c("Agricultura", "Agua", "Arbustal de paramo", "Bosque de paramo",
"Herbazal de paramo", "Humedal", "Mineria", "Pastos", "Pino",
"Suelo descubierto")
classdf <- data.frame(classvalue1 = c(1, 2, 3, 4, 5, 6, 7, 8,
9, 10), classnames1 = nlcdclass)
classcolor <- c("lightgreen", "deepskyblue", "darkolivegreen1",
"darkolivegreen", "mediumaquamarine", "dimgray", "firebrick4",
"green1", "lawngreen", "khaki2")
plot(pr2011, col = classcolor)
summary(pr2011)## layer
## Min. 1
## 1st Qu. 3
## Median 5
## 3rd Qu. 5
## Max. 10
## NA's 122356freq(pr2011)## value count
## [1,] 1 9319
## [2,] 2 2659
## [3,] 3 20746
## [4,] 4 1561
## [5,] 5 28782
## [6,] 6 4745
## [7,] 8 3261
## [8,] 9 5471
## [9,] 10 3964
## [10,] NA 122356Clase
Agricultura 1 Agua 2 Arbustal 3 Bosque 4 Herbazal 5 Humedal 6 Mineria 7 Pastos 8 Pino 9 Suelo 10
Para poder comparar y evaluar los métodos de clasificación se realiza para cada uno de ellos su respectiva matriz de confusión.
#### **evaluacion**
sampvals <- extract(landsat_clip, sample, df = TRUE)
sampvals <- sampvals[, -1]
classvale = data.frame(classvalue = sample$COBER)
sampdata <- data.frame(classvalue = sample$COBER, sampvals)
library(dismo)
set.seed(99)
j <- kfold(sampdata, k = 10, by = sampdata$classvalue)
table(j)## j
## 1 2 3 4 5 6 7 8 9 10
## 78 76 79 77 79 77 75 79 76 76x <- list()
for (k in 1:10) {
train <- sampdata[j != k, ]
test <- sampdata[j == k, ]
cart <- rpart(as.factor(classvalue) ~ ., data = train, method = "class",
minsplit = 4)
pclass <- predict(cart, test, type = "class")
# create a data.frame using the reference and prediction
x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c("observed", "predicted")
conmatsuper <- table(y)
# colnames(conmatsuper) <- classdf$classvalue
# rownames(conmatsuper) <- classdf$classnames1
conmatsuper## predicted
## observed 1 2 3 4 5 6 8 9 10
## 1 47 0 0 2 3 4 26 0 4
## 2 0 89 0 0 2 0 0 0 0
## 3 0 0 87 2 4 7 0 16 0
## 4 1 0 5 39 0 3 0 0 4
## 5 2 0 18 0 91 9 0 0 0
## 6 1 0 10 0 2 58 0 0 1
## 7 0 2 0 0 3 1 0 0 0
## 8 15 0 0 0 0 0 65 0 0
## 9 0 0 21 0 0 0 0 54 0
## 10 2 0 40 2 5 0 1 0 24library(rasterVis)
sampvals <- extract(knr, sample, df = TRUE)
sampdata <- data.frame(classvalue = sample$COBER, sampvals)
library(dismo)
set.seed(99)
j <- kfold(sampdata, k = 10, by = sampdata$classvalue)
table(j)## j
## 1 2 3 4 5 6 7 8 9 10
## 78 76 79 77 79 77 75 79 76 76x <- list()
for (k in 1:10) {
train <- sampdata[j != k, ]
test <- sampdata[j == k, ]
cart <- rpart(as.factor(classvalue) ~ ., data = train, method = "class",
minsplit = 4)
pclass <- predict(cart, test, type = "class")
# create a data.frame using the reference and prediction
x[[k]] <- cbind(test$classvalue, 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 1 2 3 4 5 6 8 9 10
## 1 85 1 0 0 0 0 0 0 0
## 2 2 88 1 0 0 0 0 0 0
## 3 2 0 113 1 0 0 0 0 0
## 4 2 0 1 48 1 0 0 0 0
## 5 2 0 0 0 118 0 0 0 0
## 6 2 0 0 0 0 69 1 0 0
## 7 6 0 0 0 0 0 0 0 0
## 8 2 0 0 0 0 0 77 1 0
## 9 2 0 0 0 0 0 0 72 1
## 10 2 0 0 0 0 0 0 0 72conmatsuper## predicted
## observed 1 2 3 4 5 6 8 9 10
## 1 47 0 0 2 3 4 26 0 4
## 2 0 89 0 0 2 0 0 0 0
## 3 0 0 87 2 4 7 0 16 0
## 4 1 0 5 39 0 3 0 0 4
## 5 2 0 18 0 91 9 0 0 0
## 6 1 0 10 0 2 58 0 0 1
## 7 0 2 0 0 3 1 0 0 0
## 8 15 0 0 0 0 0 65 0 0
## 9 0 0 21 0 0 0 0 54 0
## 10 2 0 40 2 5 0 1 0 24A partir de la matriz de confusión se genera la exactitud global (EG), exactitud del usuario (EU), y la exactitud del productor (EG), y los resultados se presentan a continuación.
matriz = read.table("libro1.csv", header = T, sep = ";")
matriz## Nosupervisada EU EP EG Supervisada EU.1 EP.1 EG.1
## 1 Agricultura 0,988 0,794 0,988 Agricultura 0,616 0,663 0,616
## 2 Agua 0,967 0,989 0,967 Agua 0,978 0,978 0,978
## 3 Arbustal 0,974 0,983 0,974 Arbustal 0,845 0,53 0,845
## 4 Bosque 0,923 0,98 0,923 Bosque 0,75 0,813 0,75
## 5 Herbazal 0,983 0,992 0,983 Herbazal 0,758 0,82 0,758
## 6 Humedal 0,958 1 0,958 Humedal 0,847 0,813 0,847
## 7 Mineria 0 0 0 Minería 0 0 0
## 8 Pastos 0,963 0,987 0,096 Pastos 0,768 0,733 0,72
## 9 Pino 0,96 0,986 0,96 Pino 0,72 0,844 0,72
## 10 Suelo 0,973 0,986 0,973 Suelo 0,311 0,676 0,311###Discusión
Ya que el objetivo principal de este estudio es comparar los métodos de clasificación, se revisó la información obtenida de la tabla matriz para determinar cuál de los dos métodos muestra mejores resultados, en cuanto a la clasificación no supervisada en los valores de las diferentes exactitudes se tienen que los porcentajes en todos los casos están por encima del 90% determinando en el caso por ejemplo de la exactitud del usuario que más del 90% de los pixeles clasificados en las diferentes clases, realmente lo son en terreno.
Caso contrario el que paso en la clasificación supervisada, aunque los valores de exactitud no son bajos al compararlos si existe una diferencia, esto puede deberse a que como se mencionó para la clasificación no supervisada se utilizó el NDVI por la naturaleza de la zona y del estudio, para la supervisada los datos de entrenamiento se generaron a partir de composiciones de falso color en donde para las clases de vegetación no llegaban al ser de todo representativas para que el método pudiera separarlas y clasificarlas de una mejor manera, también los datos de muestreo no llegan a ser lo suficiente pues para todas la clases se tienen alrededor de 60 a 100 puntos de entrenamiento, siendo idealmente una cifra mucho mayor la que pueda llegar a mejorar este clasificador, siendo tal que en los dos clasificadores la clase minería con solo 6 puntos de entrenamiento no arrojo resultados de estar clasificada, lo que indica que también hay que mejorar los puntos de entrenamientos y re muestrearlos.
Los mejores resultados los obtiene la clasificación no supervisada, esto no significa que sea mejor método clasificador que la supervisada, simplemente que para este estudio por la naturaleza, los datos utilizados tiene mejor desempeño.
###Conclusión
Los métodos basados en percepción remota para la clasificación de coberturas y uso de la tierra son el estudio o base para el estudio de acciones que lleven a una correcta planificacion del territorio ya sea bien por estudios multitemporales de cambio o de presión antrópica por nombrar algunos, en especial en ecosistemas como páramos que por condiciones geográficas, ambientales y sociales son de relevancia en la actualidad.
La clasificación no supervisada por los resultados obtenidos en la evaluación de exactitud muestra que es eficiente para los datos que se utilizaron, en cuanto a la clasificación supervisada se podría concluir que aunque no ser valores bajos en la evaluación de la exactitud, no es el más óptimo en este caso, se puede mejorar con mejores puntos de entrenamiento o con una mejor resolución espacial, también para estudios futuros se podría evaluar qué pasa si los puntos de entrenamiento se sacan a partir de una fusión de imágenes para mejor los resultados del clasificador ya que en este estudio no se tuvo en cuenta esta opción. Ambos clasificadores obtendría mejores resultados si solo se tienen en cuenta 9 clases eliminando la clase de minería ya que aunque se pensaba podría arrojar algunos resultados que hubieran sido interesantes analizar, lo que se reflejó en la pantalla no alcanza a ser significativo para que sea clasificado en ninguno de los dos métodos, pero si sabe de actividad minera en esta zona de estudio pero afortunadamente no es alcanza ser descrita en los resultados probablemente con una mejor resolución espacial esta clase si tenga representatividad en los clasificadores.
###Bibliografía
INSTITUTO GEOGRÁFICO AGUSTÍN CODAZZI – IGAC. (1979). PRORADAM: Proyecto Radargramétrico del Amazonas. República de Colombia, Bogotá. 590p.
Quichimbo, P., Tenorio, G., Borja, P., Cárdenas, I., Crespo, P., Célleri, R., & Célleri, R. (2012). Efectos sobre las propiedades físicas y químicas de los suelos por el cambio de la cobertura vegetal y uso del suelo: páramo de Quimsacocha al sur del Ecuador. Suelos Ecuatoriales, 42(2), 138-153.
Martínez, M. L., Pérez O., Vázquez, G., Castillo G., Garcia J., Mehltreter, K., Equihua, M., & Landgrave, R. (2009). Effects of land use change on biodiversity and ecosystem services in tropical montane cloud forests of Mexico,. In: Forest Ecol. Manag.p 258.
Ward, A., Dargusch, P., Grussu, G. & Romeo, R. (2015) Using carbón finance to support climate policy objectives in high mountain ecosystems. In: Climate Policy. 1-20. Se realiza la relación entre las bandas para poder analizarlas y se crea el stack de las bandas para poder realizar el cálculo NDVI y recortar la imagen al área de estudio. Calculo NDVI.