library(raster)
## Loading required package: sp
raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- stack(raslist)
landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]
# INDICE DE VEGETACION.
#Para calcular el indice de vegetacion se utiliza el objeto raster img, Es un objeto raster de multiples capas.
#Los indices de las capas son i y k
vi <- function(img, k, i)
{bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)}
#
# Para Landsat Infrarojo = 5, Rojo = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "LANDSAT - INDICE DE VEGETACION")
#Otra manera de conseguirlo es:
vi2 <- function(x, y)
{(x - y) / (x + y)}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="LANDSAT - INDICE DE VEGETACION")
# Para calcular la distribucion de datos contenido en un raster, utilizamos la funcion hist.
# El histograma permite encontar datos atipicos e incorrectos
hist(ndvi,
main = "Distribucion de valores NDVI",
xlab = "Indice de Vegetacion - NDVI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
#UMBRALIZACION
# Los indices de vegetacion estan estandarizados y pueden variar entre -1 y 1.Los valores mas altos indican mayor catidad de covertura o mas verde.
#Los valores por encima de 0.4, indican que definitivamente existe covertura verde.
veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetacion')
#Los datos por debajo de 0,4. que se encontaron ene el histograma (2,25 - 0,3)
land <- reclassify(ndvi, c(-Inf, 0.25, NA, 0.25, 0.3, 1, 0.3, Inf, NA))
plot(land, main = 'Desconocido')
# Estas areas desconocidas, pueden ser areas abiertas, para obtener un mayor conociminto sobre lo que es posible encontrar en ellas, se pude trazar land sobre landsatRGB.
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
plot(land, add=TRUE, legend=FALSE)
#También se pude hacer una clasificación de acuerdo al tipo de vegetación.
vegc <- reclassify(ndvi, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI based thresholding')
#ANALISIS DE COMPONENTES MULTIESPECTRALES
#Los datos multiespectrales a veces se transforman para ayudar a reducir la dimensionalidad y el ruido en los datos. La transformación de componentes principales es un método genérico de reducción de datos que se puede utilizar para crear algunas bandas no correlacionadas a partir de un conjunto más grande de bandas correlacionadas.
set.seed(1)
sr = sampleRandom(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")
#Diagrama de vegetacion
pca <- prcomp(sr, scale = TRUE)
pca
## Standard deviations (1, .., p=11):
## [1] 2.52687022 1.40533428 1.09901821 0.92507423 0.53731672 0.42657919
## [7] 0.28002273 0.12466139 0.09031384 0.04761419 0.03609857
##
## Rotation (n x k) = (11 x 11):
## PC1 PC2 PC3 PC4
## LC08_044034_20170614_B1 0.2973469 -0.3516438 -0.29454767 -0.06983456
## LC08_044034_20170614_B2 0.3387629 -0.3301194 -0.16407835 -0.03803678
## LC08_044034_20170614_B3 0.3621173 -0.2641152 0.07373240 0.04090884
## LC08_044034_20170614_B4 0.3684797 -0.1659582 0.10260552 0.03680852
## LC08_044034_20170614_B5 0.1546322 0.1816252 0.71354112 0.32017620
## LC08_044034_20170614_B6 0.3480230 0.2108982 0.23064060 0.16598851
## LC08_044034_20170614_B7 0.3496281 0.2384417 -0.11662258 0.07600209
## LC08_044034_20170614_B8 0.3490464 -0.2007305 0.08765521 0.02303421
## LC08_044034_20170614_B9 0.1314827 0.1047365 0.33741447 -0.92325315
## LC08_044034_20170614_B10 0.2497611 0.4912132 -0.29286315 -0.02950655
## LC08_044034_20170614_B11 0.2472765 0.4931489 -0.29515754 -0.03176549
## PC5 PC6 PC7 PC8
## LC08_044034_20170614_B1 0.49474685 0.175510232 -0.23948553 0.215745065
## LC08_044034_20170614_B2 0.22121122 0.094184121 0.06447037 0.216537517
## LC08_044034_20170614_B3 0.08482031 0.009040232 0.30511210 -0.518233675
## LC08_044034_20170614_B4 -0.33835490 -0.066844529 0.60174786 0.012437959
## LC08_044034_20170614_B5 0.51960822 -0.059561476 -0.07348455 -0.083217504
## LC08_044034_20170614_B6 -0.29437062 0.317984598 -0.02106132 0.632178645
## LC08_044034_20170614_B7 -0.25404931 0.525411720 -0.40543545 -0.478543437
## LC08_044034_20170614_B8 -0.31407992 -0.673584139 -0.52642131 0.003527306
## LC08_044034_20170614_B9 0.03040161 0.059642905 -0.03152221 -0.002775800
## LC08_044034_20170614_B10 0.16317572 -0.243735973 0.14341520 0.041736319
## LC08_044034_20170614_B11 0.19294569 -0.241611777 0.11997475 -0.022446494
## PC9 PC10 PC11
## LC08_044034_20170614_B1 0.122812108 0.535959306 0.1203473847
## LC08_044034_20170614_B2 0.091063964 -0.773112627 -0.1817872036
## LC08_044034_20170614_B3 -0.644305383 0.070860458 0.0540175730
## LC08_044034_20170614_B4 0.543822097 0.218324141 0.0135097158
## LC08_044034_20170614_B5 0.209682702 -0.040186292 -0.0004965182
## LC08_044034_20170614_B6 -0.397210135 0.089423617 -0.0045305608
## LC08_044034_20170614_B7 0.248871690 -0.074907393 0.0003958921
## LC08_044034_20170614_B8 -0.012642132 0.002411975 0.0007749022
## LC08_044034_20170614_B9 0.003176077 -0.000832654 0.0019986985
## LC08_044034_20170614_B10 -0.003574004 -0.158727672 0.6900281980
## LC08_044034_20170614_B11 -0.043475408 0.148102091 -0.6878990264
screeplot(pca)
pci <- predict(landsat, pca, index = 1:2)
plot(pci[[1]])
#Para entender lo que destaca el segundo componente principal. Vamos a intentar umbral de nuevo:
pc2 <- reclassify(pci[[2]], c(-Inf,0,1,0,Inf,NA))
par(mfrow = c(1,2))
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plot(pc2, legend = FALSE, add = TRUE)
