INTRODUCCIÓN

Este cuaderno ilustra cómo calcular estadísticas de imágenes unibanda y multibanda utilizando la libreria raster.

La imagen satelital usada se denomina “Landsat 8 "LC08_L1TP_022031_20200831_20200906_01_T1.TIF”, pertenece a la zona sur del lago Michigan, cubre parcialmente los estados de Illionis, Indiana y Michigan.

Se cargan las librerias que se necesitan para elaborar este proyecto:

library(raster)
## Loading required package: sp

Se carga la ruta donde se encuentra la imagen y se visualiza

setwd("G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1")
getwd()
## [1] "G:/PERCEPCION_REMOTA/TALLER/LABORATORIO_2/LC08_L1TP_022031_20200831_20200906_01_T1"

Se descarga cada banda y se guarda en una variable SpatRaster.

# Ultra azul: Estudios costeros y de aerosoles
b1 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF')
# Azul:Cartografía batimétrica, que distingue el suelo de la vegetación y la vegetación caducifolia de la vegetación de coníferas.
b2 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF')
# Verde: Destaca los picos de máxima vegetación, que son útiles para evaluar el vigor de las plantas.
b3 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF')
# Rojo: Distingue las laderas de vegetación
b4 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B4.TIF')
# Infrarrojo Cercano (NIR): Destaca el contenido de biomasa y las costas
b5 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B5.TIF')
# Infrarrojo de Onda Corta 1 (SWIR 1):Distingue la humedad del suelo y de la vegetación; penetra a través de nubes finas
b6 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B6.TIF')
# Infrarrojo de Onda Corta 2 (SWIR 2):Mejora de la lectura de la humedad del suelo y la vegetación y la penetración a través de nubes finas
b7 <- raster('G:\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B7.TIF')

Se crea un _Rasterstac_que permite crear un es un objeto ráster de varias capas, en este caso las capas son las bandas de la imagen satelital. Se revisa las caracteristicas del vector “landsat20” donde almacenamos las 7 bandas

landsat20 <-stack(b1, b2, b3,b4,b5,b6,b7)
landsat20
## class      : RasterStack 
## dimensions : 7871, 7751, 61008121, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 403485, 636015, 4504485, 4740615  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## names      : LC08_L1TP//6_01_T1_B1, LC08_L1TP//6_01_T1_B2, LC08_L1TP//6_01_T1_B3, LC08_L1TP//6_01_T1_B4, LC08_L1TP//6_01_T1_B5, LC08_L1TP//6_01_T1_B6, LC08_L1TP//6_01_T1_B7 
## min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values :                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535

A cada banda se le asigna un nombre

names(landsat20) <- c('Ultra Azul', 'Azul', 'Verde', 'Rojo', 'NIR', 'SWIR1', 'SWIR2')
landsat20
## class      : RasterStack 
## dimensions : 7871, 7751, 61008121, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 403485, 636015, 4504485, 4740615  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## names      : Ultra.Azul,  Azul, Verde,  Rojo,   NIR, SWIR1, SWIR2 
## min values :          0,     0,     0,     0,     0,     0,     0 
## max values :      65535, 65535, 65535, 65535, 65535, 65535, 65535

Calculo de la reflectancia TOA

Los valores digitales están en el rango de 0 a 65535. Los valores de reflectancia aparente se calculan a partir de los metadatos que vienen junto a las imágenes (archivo txt). El cálculo de reflectancia a tope de atmósfera (TOA) utiliza el factor multiplicativo específico por banda y un factor de re escalamiento aditivo específico de la banda. En este caso las 7 bandas tienen la misma resolución (30 metros) y la misma proyección (Zona 16) es así que todas las 7 bandas se multiplican por los mismos valores.

toa20 <-  -0.1 + landsat20*2.0E-5
toa20
## class      : RasterBrick 
## dimensions : 7871, 7751, 61008121, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 403485, 636015, 4504485, 4740615  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/anton/AppData/Local/Temp/RtmpuC8c0X/raster/r_tmp_2021-03-24_160341_6556_20648.grd 
## names      : Ultra.Azul,    Azul,   Verde,    Rojo,     NIR,   SWIR1,   SWIR2 
## min values :       -0.1,    -0.1,    -0.1,    -0.1,    -0.1,    -0.1,    -0.1 
## max values :    0.86540, 0.92802, 1.03274, 1.07986, 1.21070, 1.21070, 1.21070

Calculo de estadísticas unibanda

Media

Es el promedio de los valores de reflectancia por banda

toa20_mean <- cellStats(toa20, 'mean')
toa20_mean
##  Ultra.Azul        Azul       Verde        Rojo         NIR       SWIR1 
## 0.037096503 0.024998087 0.013370702 0.002873349 0.098064068 0.039153471 
##       SWIR2 
## 0.006372616

Desviación estándar

Indica que tan concentrados esan los valores de reflectancia.

toa20_sd <- cellStats(toa20, 'sd')
toa20_sd
## Ultra.Azul       Azul      Verde       Rojo        NIR      SWIR1      SWIR2 
## 0.09687234 0.08946602 0.08272626 0.07788505 0.17263015 0.11450179 0.08404920

Visualización de imágenes

Para la visualización de la reflectancia, el histograma permite observar los valores de reflectancia con su frecuencia de aparición, es decir, la cantidad de veces que aparece un valor en una banda. La función strech (estiramiento lineal de valores en un objeto ráster) para su representación y su visualización.

par(mfrow=c(1, 2))

tmp <- stretch(toa20[[1]], minv=0, maxv=255, minq=0, maxq=1)
plot(tmp, col = gray(0:100 / 100), main='Landsat 8  B1 - linear stretch', legend=FALSE, axes=FALSE)
tmp <- stretch(toa20[[2]], minv=0, maxv=255, minq=0, maxq=1)
plot(tmp, col = gray(0:100 / 100), main='Landsat 8 B2 - linear stretch', legend=FALSE, axes=FALSE)

El histogramas de la banda 1 es graficado con el 3% del total de los Niveles Digitales de reflectancia, se encuentran concentrados asi: -0.2 y -0.1 con una frecuencia de 3200, entre 0.09 y 0.1 con una frecuencia mayor a 4000, entre 0.1 y 0.2 con una frecuencia de 2000 y entre 0.15 a 0.3 frecuencias menores a 100.

El histogramas de la banda 2 es graficado con el 3% del total de los Niveles Digitales de reflectancia, se encuentran concentrados asi: -0.2 y -0.1 con una frecuencia de 3100, entre 0.09 y 0.1 con una frecuencia a 6000, entre 0.1 y 0.2 con una frecuencia de 2000 y entre 0.15 a 0.3 frecuencias menores a 1000.

hist(toa20, c(1,2), maxpixels=2000000, plot=TRUE, main="Histograma de la B1 y B2", xlab= "TOA reflectancia", ylab="Frecuencia", col="violetred1")
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 3% of the raster cells were used. 2000000 values used.

## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 3% of the raster cells were used. 2000000 values used.

par(mfrow=c(1, 2))
tmp <- stretch(toa20[[3]], minv=0, maxv=255, minq=0, maxq=1)
plot(tmp, col = gray(0:1000 / 1000), main='Landsat 8 B3 - linear stretch', legend=FALSE, axes=FALSE)

tmp <- stretch(toa20[[4]], minv=0, maxv=255, minq=0, maxq=1)
plot(tmp, col = gray(0:1000 / 1000), main='Landsat 8 B4 - linear stretch', legend=FALSE, axes=FALSE)

Histogramas de las bandas 3, 4, 5, 6 y 7

En los histogramas 3,4,5 el comportamiento de los datos es similar, su distribución los niveles digitales aparecen en un puntos en concreto (bimodales).

El histograma de la banda NIR y los infrarrojos de onda corta tienen una mayor acomulación de pixeles por zonas y con un rango mas amplio a comparación de las anteriores bandas.

par(mfrow=c(3, 3))
hist(toa20, c(3,4), maxpixels=10000, plot=TRUE, main="Histogramas de las B3 y B4", xlab= "TOA reflectancia", ylab="Frecuencia", col="violetred1")
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.

## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.
hist(toa20, c(4,5), maxpixels=10000, plot=TRUE, main="Histogramas de las B5 y B5", xlab= "TOA reflectancia", ylab="Frecuencia", col="violetred1")
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.

## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.
hist(toa20, c(6,7), maxpixels=10000, plot=TRUE, main="Histogramas de las B6 y B7", xlab= "TOA reflectancia", ylab="Frecuencia", col="violetred1")
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.

## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 10000 values used.

Curtosis

Representa el achatamiento de la distribucion de los datos de la variable TOA 20. Forma de la distribucón de la probabilidad.

library(moments)
toa20kurtosis<-kurtosis(toa20)
toa20kurtosis
## [1] 3.667091

ASIMETRIA

Devuelve la asimetria de los datos de la variable TOA20, es decir, la distribución de la probabilidad de una variable aleatoria de valor real sobre su media.

toa20skewness<-skewness(toa20)
toa20skewness
## [1] 0.6289674

Cálculo de estadísticas multibanda

Se grafican dos imágenes RGB, la primera con las bandas 432 y la segunda con las bandas 564.

Las estadísticas multibanda incluyen: la matriz de covarianza la matriz de correlación diagramas de dispersión por pares, suele ir acompañado del coeficiente de correlación.

par(mfrow=c(1,2))

plotRGB(toa20, r = 4, g = 3, b = 2,main='L8 2020 RGB432', stretch='lin',  colNA='white')
plotRGB(toa20, r = 5, g = 6, b = 4,main='L8 2020 RGB564', stretch='lin', colNA='white')

Correlación de Pearson

La correlación de Pearson mide la relación estadística entre dos bandas. Puede tomar valores entre el rango -1 y 1.

pairs(toa20[[1:3]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=10000)

En este grafico se observa la correlación de asociación positiva entre las bandas ultra azul, azul y verde.

pairs(toa20[[4:6]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=10000)

En la correlación de pearson entre las bandas rojo, NIR, SWIR1. se observa que los valores de correlación son positivos en su mayoría. Entre la banda NIR y la banda Rojo los puntos no se comportan de manera lineal y algunos de estos ND están en valores negativos.

Covarianza

La covarianza la medida de varianza entre dos varibles.

covmat<-cov(as.matrix(toa20))
round(covmat*1000,2)
##            Ultra.Azul  Azul Verde  Rojo   NIR SWIR1 SWIR2
## Ultra.Azul       9.38  8.66  7.92  7.34 13.01  9.58  7.47
## Azul             8.66  8.00  7.35  6.84 12.01  8.93  6.97
## Verde            7.92  7.35  6.84  6.41 11.75  8.68  6.67
## Rojo             7.34  6.84  6.41  6.07 10.80  8.21  6.34
## NIR             13.01 12.01 11.75 10.80 29.80 18.47 12.61
## SWIR1            9.58  8.93  8.68  8.21 18.47 13.11  9.44
## SWIR2            7.47  6.97  6.67  6.34 12.61  9.44  7.06
#boxplot(toa20[[4:6]], maxpixels=100000)
cmat <- cor(as.matrix(toa20))
round(cmat, 2)
##            Ultra.Azul Azul Verde Rojo  NIR SWIR1 SWIR2
## Ultra.Azul       1.00 1.00  0.99 0.97 0.78  0.86  0.92
## Azul             1.00 1.00  0.99 0.98 0.78  0.87  0.93
## Verde            0.99 0.99  1.00 0.99 0.82  0.92  0.96
## Rojo             0.97 0.98  0.99 1.00 0.80  0.92  0.97
## NIR              0.78 0.78  0.82 0.80 1.00  0.93  0.87
## SWIR1            0.86 0.87  0.92 0.92 0.93  1.00  0.98
## SWIR2            0.92 0.93  0.96 0.97 0.87  0.98  1.00
cor.test(x=as.matrix(toa20[[4]]), y=as.matrix(toa20[[5]]), method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  as.matrix(toa20[[4]]) and as.matrix(toa20[[5]])
## t = 10540, df = 61008119, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8033358 0.8035137
## sample estimates:
##       cor 
## 0.8034248

Gráficos

Este gráfico representa la relación de índices de reflectancia entre bandas. La banda NIR muestra una correlación más baja al ser comparado entre bandas (tono más claro).

library(corrplot)
## corrplot 0.84 loaded
ncmat <- round(cmat, 2)
corrplot(ncmat, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

tmp1 <- as.matrix(sample(toa20,500000))
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following objects are masked from 'package:moments':
## 
##     kurtosis, skewness
## The following object is masked from 'package:graphics':
## 
##     legend

La matriz de correlación puede ser representada por un diagrama de dispersión. El valor de la correlación mas el nivel de significacia se ve representado por estrellas. Se observa que en general todas las bandas tienen una correlación positiva.

chart.Correlation(tmp1, histogram=TRUE, pch=19)

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] PerformanceAnalytics_2.0.4 xts_0.12.1                
## [3] zoo_1.8-9                  corrplot_0.84             
## [5] moments_0.14               raster_3.4-5              
## [7] sp_1.4-5                  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        quadprog_1.5-8    codetools_0.2-18  lattice_0.20-41  
##  [5] digest_0.6.27     grid_4.0.4        magrittr_2.0.1    evaluate_0.14    
##  [9] highr_0.8         rlang_0.4.10      stringi_1.5.3     rmarkdown_2.7    
## [13] rgdal_1.5-23      tools_4.0.4       stringr_1.4.0     xfun_0.22        
## [17] yaml_2.2.1        compiler_4.0.4    htmltools_0.5.1.1 knitr_1.31