Lab 3B Basic mathematic operations
Zuly Olaya
20/3/2021
En este cuaderno se prentende explicar la utilización de indices vegetales en imagenes Landsat 8. apoyado en el cuaderno el profesor Ivan Lizarazo: https://rspatial.org/terra/rs/3-basicmath.html
Operaciones matemáticas básicas
Se cargala las libreria Terra la ruta de trabajo.
rm(list=ls())library(terra)## terra version 1.1.4##
## Attaching package: 'terra'## The following object is masked from 'package:knitr':
##
## spinlibrary(wesanderson)
library("rgdal")## Loading required package: sp## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/anton/OneDrive/Documentos/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/anton/OneDrive/Documentos/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/anton/OneDrive/Documentos/R/win-library/4.0/rgdal/proj##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(wesanderson)
knitr::opts_chunk$set(echo = TRUE)
setwd("G:\\LABORATORIO_3B\\LC08_L1TP_022031_20200831_20200906_01_T1_1")
getwd()## [1] "G:/LABORATORIO_3B/LC08_L1TP_022031_20200831_20200906_01_T1_1"list.files()## [1] "B1.TIF"
## [2] "B10.TIF"
## [3] "B11.TIF"
## [4] "B2.TIF"
## [5] "B3.TIF"
## [6] "B4.TIF"
## [7] "B5.TIF"
## [8] "B6.TIF"
## [9] "B7.TIF"
## [10] "BQA.TIF"
## [11] "LABORATORIO_3B.html"
## [12] "LABORATORIO_3B.Rmd"
## [13] "LC08_L1TP_022031_20200831_20200906_01_T1_1"
## [14] "LC08_L1TP_022031_20200831_20200906_01_T1_ANG.txt"
## [15] "MTL.txt"Se cargan las bandas del 1 al 7 la imagen Landsat
Se utiliza la función rast , se guarda en “landsat” luego se confiorma: “landsatRGB” que contiene las bandas Rojo, Verde y Azul y “landsatFCC” que contien las bandas NIR,Rojo y Verde.
rfiles <- paste0('B',1: 7,".TIF")
rfiles## [1] "B1.TIF" "B2.TIF" "B3.TIF" "B4.TIF" "B5.TIF" "B6.TIF" "B7.TIF"list.files()## [1] "B1.TIF"
## [2] "B10.TIF"
## [3] "B11.TIF"
## [4] "B2.TIF"
## [5] "B3.TIF"
## [6] "B4.TIF"
## [7] "B5.TIF"
## [8] "B6.TIF"
## [9] "B7.TIF"
## [10] "BQA.TIF"
## [11] "LABORATORIO_3B.html"
## [12] "LABORATORIO_3B.Rmd"
## [13] "LC08_L1TP_022031_20200831_20200906_01_T1_1"
## [14] "LC08_L1TP_022031_20200831_20200906_01_T1_ANG.txt"
## [15] "MTL.txt"landsat <- rast(rfiles)
landsat## class : SpatRaster
## dimensions : 7871, 7751, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## sources : B1.TIF
## B2.TIF
## B3.TIF
## ... and 4 more source(s)
## names : B1, B2, B3, B4, B5, B6, ...landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]Índices de vegetación
En la función img un objeto SpatRaster de múltiples capas i y k son los índices de las capas (números de capa) que se utilizan para calcular el índice 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)
}Calcular el índice de vegetación de diferencia normalizada (NDVI)
Es un indice que permite estimar el desarrollo de la vegetación, su calidad, determinar su estado en general y cuantificación, ademas de diferenciar la vegetación de otros tipos de cobertura del suelo. Su formula general en dond NIR es luz infrarroja cercana y Red es luz roja visible.
NDVI= (NIR-RED) / (NIR + RED)
# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col=rev(terrain.colors(10)), main = "Landsat-NDVI")
La respuesta espectral muestra la vegetación reasaltada en color verde la absorción de pigmento de clorofila en la banda roja y la alta reflectividad del material de las plantas en la banda infrarroja cercana (NIR).
Una forma alternativa de lograr esto es así:
# Write a general function that can compute 2-layer NDVI type indices
vi2 <- function(x, y) {
(x - y) / (x + y)
}
nir <- landsat[[5]]
red <- landsat[[4]]
ndvi2 <- lapp(c(nir, red), fun = vi2) #superposición
# or in one line
ndvi2 <- lapp(landsat[[5:4]], fun=vi2)#superposición
plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")# 10 tonalidades de colores
Índice Diferencial de Agua Normalizado (NDWI)
El NDWI se utiliza como una medida de la cantidad de agua que posee la vegetación o el nivel de saturación de humedad que posee el suelo. Para obtener este índice la combinación de bandas es la siguiente:
(Green-SWIR1)/(Green+SWIR1)
# For Landsat Green = 3, swir1 = 6.
ndvi <- vi(landsat, 3, 6)
plot(ndvi, col=rev(topo.colors(1000)), main = "Landsat-NDWI")
# Write a general function that can compute 2-layer NDWI type indices
vi3 <- function(x, y) {
(x - y) / (x + y)
}
green <- landsat[[3]]
swir1 <- landsat[[6]]
ndvi3 <- lapp(c(green, swir1), fun = vi3) #superposición
# or in one line
plot(ndvi3, col=rev(topo.colors(1000)), main="Landsat-NDWI")# 1000 tonalidades de colores con la paleta topo.colors ### NDBI o Índice de Diferencia Normalizada Edificada
Permite llevar a cabo la estimación de zonas con superficies edificadas o en desarrollo de construcción frente a las habituales zonas naturalizadas con vegetación o desnuda.Requiere de las bandas de análisis del infrarrojo a través de las bandas SWIR y NIR. NDBI= (SWIR-NIR)/ (SWIR+NIR)
# Write a general function that can compute 2-layer NDBI type indices
vi4 <- function(x, y) {
(x - y) / (x + y)
}
swir1 <- landsat[[6]]
NIR <- landsat[[5]]
ndvi4 <- lapp(c(swir1, NIR), fun = vi4) #superposición
# or in one line
plot(ndvi4, col=rev(topo.colors(8)), main="Landsat-NDBI")
Las edificaciones se ven resaltadas de color verde claro. Se encuentran concentradas al lado izquierdo del lago Michigan.
Histograma
Podemos explorar la distribución de valores contenidos dentro de nuestro ráster usando la función hist() que produce un histograma. Los histogramas a menudo son útiles para identificar valores atípicos y valores de datos incorrectos en nuestros datos ráster.
# view histogram of data
hist(ndvi,
main = "Distribution of NDVI values",
xlab = "NDVI",
ylab= "Frequency",
col = "steelblue",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n',
maxcell=50000000)## Warning: [hist] 82% of the raster cells were used. (of which 32% was NA)axis(side=1, at = seq(-0.5,0.5, 0.05), labels = seq(-0.5,0.5, 0.05))
El histograma refleja los indice de vegetación de diferencia normalizada concentrados entre -0.05 y-0.04 con una frecuencia mayor de 200000.
par(mfrow=c(2, 2))
# view histogram of data
hist(ndvi2,
main = "Distribution of NDWI values",
xlab = "NDWI",
ylab= "Frequency",
col = "steelblue",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning: [hist] 2% of the raster cells were used. (of which 32% was NA)axis(side=1, at = seq(-0.5,0.7, 0.05), labels = seq(-0.5,0.7, 0.05))
# view histogram of data
hist(ndvi3,
main = "Distribution of NDBI values",
xlab = "NDBI",
ylab= "Frequency",
col = "steelblue",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning: [hist] 2% of the raster cells were used. (of which 32% was NA)axis(side=1, at = seq(-0.5,0.7, 0.05), labels = seq(-0.5,0.7, 0.05))
El histograma NDWI refleja los indices concentrados -0.1 y 0 confrecuencia mayor de 200000 y el histograma de NDBI mestra mas picos en el histograma entre los que se resaltan estan entre -0.2 a 0 y entre 0.1 a 0.2 siendo un histograma bimodal.
Umbral
Podemos aplicar reglas básicas para obtener una estimación de la extensión espacial de diferentes características de la superficie terrestre. Tenga en cuenta que los valores de NDVI están estandarizados y oscilan entre -1 y +1. Los valores más altos indican más cobertura verde.
Las celdas con valores de NDVI superiores a 0,4 son definitivamente vegetación. La siguiente operación enmascara todas las celdas que quizás no sean vegetación (NDVI <0.4).
agua <- classify(ndvi, cbind(-Inf, 0.1, NA))
plot(agua, main='Cuerpos de agua')
Tracemos el área que corresponde al pico entre 0.25 y 0.3 en el histograma NDVI.
m <- c(-Inf, 0.25, NA, 0.25, 0.3, 1, 0.3, Inf, NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
land <- classify(ndvi, rcl)
plot(land, main = 'Suelo Desubierto')
Estas pueden ser áreas abiertas. Puede trazar land sobre el landsatFCC ráster original para obtener más información.
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 puede crear clases para diferentes cantidades de presencia de vegetación.
m <- c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
vegc <- classify(ndvi, rcl)
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI based thresholding')
Análisis de componentes principales
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 de reducción de datos genéricos que se puede utilizar para crear algunas bandas no correlacionadas a partir de un conjunto más grande de bandas correlacionadas.
Puede calcular el mismo número de componentes principales que el número de bandas de entrada. El primer componente principal (PC) explica el mayor porcentaje de varianza y otros PC explican adicionalmente la varianza en orden decreciente.
set.seed(1)
sr <- spatSample(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")
Esto se conoce como diagrama de vegetación y línea de suelo (igual que el diagrama de dispersión en la sección anterior).
pca <- prcomp(sr, scale = TRUE)
pca## Standard deviations (1, .., p=7):
## [1] 2.55099844 0.61958475 0.30985799 0.09663712 0.04772286 0.02670859 0.01342208
##
## Rotation (n x k) = (7 x 7):
## PC1 PC2 PC3 PC4 PC5 PC6
## B1 -0.3795621 -0.3564829 0.34567684 -0.4779427 0.008108849 0.13279312
## B2 -0.3812175 -0.3541485 0.24416352 -0.1667111 0.009042713 0.21979631
## B3 -0.3882508 -0.2121425 0.06492825 0.2938115 -0.191212781 -0.82283701
## B4 -0.3866806 -0.2132277 -0.22473073 0.6937218 -0.041270339 0.47938554
## B5 -0.3467903 0.6974748 0.56096893 0.1826765 0.207739470 0.04508059
## B6 -0.3766616 0.3970910 -0.38607205 -0.3055947 -0.669737775 0.09099867
## B7 -0.3849776 0.1134718 -0.54994032 -0.2223606 0.685479782 -0.12988957
## PC7
## B1 0.602686537
## B2 -0.770358001
## B3 0.009636866
## B4 0.204314729
## B5 0.004535443
## B6 -0.034286578
## B7 -0.016851474screeplot(pca)
# function to restrict prediction to the
# first two principal components
pca_predict2 <- function(model, data, ...) {
predict(model, data, ...)[,1:2]
}
pci <- predict(landsat, pca, fun=pca_predict2)
plot(pci)
El primer componente principal destaca los límites entre las clases de uso de la tierra o los detalles espaciales, que es la información más común entre todas las longitudes de onda. Es difícil comprender qué destaca el segundo componente principal. Intentemos establecer el umbral de nuevo:
# quick look at the histogram of second component
hist <- pci[[2]]
m <- c(-Inf,-3,NA, -3,-2,0, -2,-1,1, -1,0,2, 0,1,3, 1,2,4, 2,6,5, 6,Inf,NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
rcl## [,1] [,2] [,3]
## [1,] -Inf -3 NA
## [2,] -3 -2 0
## [3,] -2 -1 1
## [4,] -1 0 2
## [5,] 0 1 3
## [6,] 1 2 4
## [7,] 2 6 5
## [8,] 6 Inf NApcClass <- classify(pci[[2]], rcl)Nuevo gráfico de resultados
par(mfrow=c(1,2))
plotRGB(landsatFCC, r = 1, g = 2, b = 3, stretch = "lin", axes=TRUE, main = "False Color")
plot(pcClass, main="PCA")
Reproductibilidad
Este cuaderno fue basado en el cuaderno realizado por el profesor Ivan Lizarazo. Se puede consultar en el siguiente enlace: https://rpubs.com/ials2un/landsat_image_stats
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Latin America.1252
## [2] LC_CTYPE=Spanish_Latin America.1252
## [3] LC_MONETARY=Spanish_Latin America.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=Spanish_Latin America.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rgdal_1.5-23 sp_1.4-5 wesanderson_0.3.6 terra_1.1-4
## [5] knitr_1.31
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 codetools_0.2-18 lattice_0.20-41 digest_0.6.27
## [5] grid_4.0.4 magrittr_2.0.1 evaluate_0.14 highr_0.8
## [9] rlang_0.4.10 stringi_1.5.3 raster_3.4-5 rmarkdown_2.7
## [13] tools_4.0.4 stringr_1.4.0 xfun_0.22 yaml_2.2.1
## [17] compiler_4.0.4 htmltools_0.5.1.1.