Ejercicio 3 - Estadística de Imagenes
Ejercicio desarrollado en Marzo 18 de 2020, en el curso de Sensores remotos, y orientado por la practica de la página https://rspatial.org/raster/rs/3-basicmath.html
Se cargan los datos para el analisis, (se utilizan los mismo del ejercicio 2)
library(raster)## Loading required package: splibrary(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)]]“————————————” INDICES DE VEGETACIÓN “————————————”
Se construye una función para la determinación del “Vegetation Index” - Ruta 1
vi <- function(img, k, i) { #"k" e "i" corresponden a las bandas a utilizar (usualmente el NIR y R)
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}Analizaremos el indice de diferencia normalizado de la vegetación para la imagen en cuestion (por lo cual se utilizara la banda NIR, Inflarojo cercano (5) y Roja (4))
# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")
Se construye una función para la determinación del “Vegetation Index” - Ruta 2
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-NDVI")
“————————————” HISTOGRAMAS “————————————”
Construimos el histograma de los valores obtenidos en el NDVI
# view histogram of data
hist(ndvi,
main = "Distribution of NDVI values",
xlab = "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))
“————————————” THRESHOLDING (UMBLAL) “————————————”
Podemos obtener una estimación de la extensión espacial de diferentes elementos en la superficie, teniendo en cuenta los valores NDVI previamente obtenidos se conoce que para valores altos (en un rango de -1 a 1) nos indican coberturas verdes (Las celdas con valores de NDVI superiores a 0.4 son definitivamente vegetación").
La siguiente operación enmascara (Mask) todas las celdas que no cumplan con ese criterio.
veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetation')
tambien se puede analizar el histograma, se graficaran los valores NVDI para el rango 0.25 y 0.3.
land <- reclassify(ndvi, c(-Inf, 0.25, NA, 0.25, 0.3, 1, 0.3, Inf, NA))
plot(land, main = 'What is it?')
Podemos superponer estos rangos sobre el raster original
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
plot(land, add=TRUE, legend=FALSE)
Se pueden crear rangos para la clasificación de vegetación basados en el indice NVDI
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')
“————————————” ANÁLISIS DE COMPONENTES PRINCIPALES “————————————”
La transformación de componentes principales es un método genérico de reducción de datos, que ayudara a crear bandas no correlacionadas de ese conjunto más grande de bandas correlacionadas.
set.seed(1)
sr <- sampleRandom(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")
Esta grafica es conocida como diagrama de vegetación y línea de suelo
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.6878990264screeplot(pca)
pci <- predict(landsat, pca, index = 1:2)
plot(pci[[1]])
Se considera que este componente principal fue poco claro, por lo cual se analiza el segundo componente principal
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)