💡 MATERIALES PARA LA ACTIVIDAD:

Los materiales para el desarrollo de esta actividad son los siguientes.

Imágenes satélite.

script_laboratorio_5_operaciones_locales

PASOS INICIALES

Para llevar a término este ejercicio, se requiere el uso de los siguientes paquetes.

library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21
# install.packages("kableExtra")
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.2

El paquete kableExtra se utiliza en el diseño de tabla con datos para su incorporación a informes o trabajos de investigación.

También deberá establecerse el directorio de trabajo.

setwd("D:/G174_2025/LABORATORIO_5_Algebra_espacial/")

Desde ese directorio, se importará el mismo fichero multibanda en formato GeoTIFF de la sesión anterior, solicitando a R que muestre sus atributos.

imagen <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_globales/L8_operaciones_globales.tif")
imagen
## class       : SpatRaster 
## dimensions  : 1395, 1848, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source      : L8_operaciones_globales.tif 
## names       : L8_oper~bales_1, L8_oper~bales_2, L8_oper~bales_3, L8_oper~bales_4 
## min values  :            9287,            8398,            7472,            7127 
## max values  :           54511,           43380,           54694,           65535

Al igual que en la sesión anterior, se debe cambiar el nombre de las bandas

names(imagen) <- c("blue", "green", "red", "NIR")
names(imagen)
## [1] "blue"  "green" "red"   "NIR"

En esa escena Landsat aparece zonas con valores ausentes, codificados como 0. Si no se avisa al programa de esta circunstancia, considerará ese 0 como un valor numérico real. Por ello es necesario especificar esta circunstancia con la función terra::NAflag(). Además, dado que se realizarán diferentes cálculos estadísticos, conviene evitar la conversión de los valores originales a notación científica (función options()).

NAflag(imagen) <- 0                                          
options(scipen = 999)

Como es frecuente, conviene representar gráficamente el objeto importado mediante una imagen en falso color:

plotRGB(imagen, 
        r=4, g=3, b=2, 
        stretch = "lin", 
        axes = TRUE)

ANÁLISIS ESTADÍSTICO

La aplicación de diversos procedimientos estadísticos para el procesamiento de imágenes satelitales permite conocer, cuantificar y mejorar las características de la imagen, bien sea para ajustar la visualización o para orientar cualquier proceso posterior.

Los estadísticos que pueden extraerse son tanto de tendencia central como de dispersión. Los de tendencia central muestran alrededor de qué valores se sitúan la mayoría de los valores de los píxeles de los objetos ráster. El fundamental es la Media o promedio, valor que expresa el promedio de los niveles digitales de todos los pixeles de la imagen (Ecuación (2.1)).

\[\bar{x} = \frac{1}{n} \sum_{i=1}^{n}x_{i}\]

donde, \(x_i\) son los niveles digitales de los píxeles, y \(n\) el número total de pixeles en la imagen.

También es posible caracterizar las imágenes mediante estadísticos de dispersión. Uno de los básicos es la *Varianza**, un estadístico que muestra la dispersión de los niveles digitales de los pixeles alrededor de la media. Si los valores se concentran alrededor de la media la varianza es baja, por lo que su histograma será estrecho; en caso contrario, la varianza será. La fórmula de la varianza es la siguiente:

\[\sigma^{2} = \frac{\sum\limits_{i=1}^{n} \left(x_{i} - \bar{x}\right)^{2}} {n-1}\] donde, \(x_i\) es el nivel digital del pixel, \(n\) el número total de pixeles en la imagen y \(\bar{x}\) la media aritmética de la imagen.

Otro estadístico de dispersión ampliamente utilizado es la Desviación típica o estándar, equivalente a la raíz cuadrada de la varianza

\[\sigma = \sqrt{\frac{\sum\limits_{i=1}^{n} \left(x_{i} - \bar{x}\right)^{2}} {n-1}}\]

Cada banda por separado

El análisis puede comenzar con la función raster::summary.

summary(imagen)
## Warning: [summary] used a sample
##       blue           green            red             NIR       
##  Min.   : 9389   Min.   : 8606   Min.   : 7741   Min.   : 7483  
##  1st Qu.:11095   1st Qu.:10994   1st Qu.:11484   1st Qu.:15839  
##  Median :11701   Median :11837   Median :12830   Median :17211  
##  Mean   :11817   Mean   :11967   Mean   :12939   Mean   :17475  
##  3rd Qu.:12346   3rd Qu.:12784   3rd Qu.:14258   3rd Qu.:18946  
##  Max.   :30426   Max.   :42671   Max.   :49305   Max.   :53159  
##  NA's   :1365    NA's   :1365    NA's   :1365    NA's   :1365

El cálculo de estos valores estadísticos se puede también obtener mediante la función terra::global. Hay que recordar que existen valores NA en el ráster, y que, si no se avisa a R de esa circunstancia devolverá siempre NaN

global(imagen, "mean", na.rm=TRUE)
##           mean
## blue  11817.60
## green 11968.05
## red   12940.29
## NIR   17474.32
global(imagen, "std" , na.rm=TRUE)
##            std
## blue  1129.481
## green 1440.882
## red   1992.220
## NIR   2259.486
global(imagen, "rms", na.rm=TRUE)                                                                     # Media cuadrática
##            rms
## blue  11871.45
## green 12054.47
## red   13092.75
## NIR   17619.80

También es posible el cálculo de cuantiles mediante la función `quantile

quantile(imagen)
## class       : SpatRaster 
## dimensions  : 1395, 1848, 5  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## names       :    q0,    q0.25,    q0.5,    q0.75,    q1 
## min values  :  7127,  7385.75,  8002.0,  8673.75,  9287 
## max values  : 43300, 50527.00, 53723.5, 57326.25, 65535

A su vez, estos valores estadísticos pueden transformarse en dataframes. Primero se obtendrá la media y la desviación típica.

imagen_promedio <- global(imagen, 
                          "mean", 
                          na.rm=TRUE)
imagen_desviacion_tipica <- global(imagen, 
                                   "sd", 
                                   na.rm=TRUE)

estadisticos <- cbind(imagen_promedio, imagen_desviacion_tipica)
colnames(estadisticos) <- c("Media", "Desviación típica")

A continuación se creará la tabla

kable(round(estadisticos,4), 
      caption="Resumen estadísticas de las bandas")
Resumen estadísticas de las bandas
Media Desviación típica
blue 11817.60 1129.481
green 11968.05 1440.883
red 12940.29 1992.221
NIR 17474.32 2259.487

La función acepta concatenar diferentes argumentos.

global(imagen,                                      # objeto SpatRaster
       c("sum", "mean", "sd"),                    # Estadísticos a calcular
       na.rm=TRUE)                                # No tiene en cuenta los valores ausentes NA
##               sum     mean       sd
## blue  30047048710 11817.60 1129.481
## green 30429573652 11968.05 1440.883
## red   32901572015 12940.29 1992.221
## NIR   44429658063 17474.32 2259.487

Terra acepta funciones aplicadas a todos los píxeles de un objeto SpatRaster. La sintaxis siempre debe enunciarse de la siguiente manera: *funcion(variable) {formula matemática}*. Para ello se recurre a la funciónterra::app()[https://rspatial.github.io/terra/reference/app.html], que aplica esa función a todas los píxeles de un ráster. Es recomendable usar la funciónterra::app()`[https://rspatial.github.io/terra/reference/app.html] para cálculos complejos, porque una llamada a esta función no implica la generación de archivos temporales.

La función terra::app() Una función personalizada podría escribirse así

borrame2 <- app(imagen, 
                 fun = function(x) {
                   2*x + 0.5*x^2 - 10})
plot(borrame2)

A continuación se define una función multiplicando todos los píxeles por 100, que es aplicada posteriormente al objeto ráster.

func <- function(x) {x * 100}

imagen_multiplicada <- app(imagen, func)
imagen_multiplicada
## class       : SpatRaster 
## dimensions  : 1395, 1848, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## names       :    blue,   green,     red,     NIR 
## min values  :  928700,  839800,  747200,  712700 
## max values  : 5451100, 4338000, 5469400, 6553500
plot(imagen_multiplicada)

Incluso se pueden sumar por separado alguna de las bandas que constituyen el ráster. Por ejemplo, sumar dos bandas requiere que cada píxel en el raster A corresponda a un píxel en el ráster B; ambos serán sumados y almacenados en un nuevo ráster de las mismas dimensiones:

borrame1 <- imagen[[1]] + imagen[[2]]
borrame1
## class       : SpatRaster 
## dimensions  : 1395, 1848, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## varname     : L8_operaciones_globales 
## name        :  blue 
## min value   : 17766 
## max value   : 97811
plot(borrame1)

Todas las bandas simultáneamente

También pueden calcularse esos mismos estadísticos a nivel de píxel de todas las capas de una imagen multibanda, creando una nueva capa raster que almacena esos nuevos datos. Existen varios procedimientos.

Por un lado, se pueden obtener de manera individual estos estadísticos aplicando una única función a todas las bandas de un objeto multibanda

min(imagen)
## class       : SpatRaster 
## dimensions  : 1395, 1848, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## name        :   min 
## min value   :  7127 
## max value   : 43300
max(imagen)
## class       : SpatRaster 
## dimensions  : 1395, 1848, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## name        :   max 
## min value   :  9287 
## max value   : 65535
mean(imagen)
## class       : SpatRaster 
## dimensions  : 1395, 1848, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## name        :    mean 
## min value   :  8146.5 
## max value   : 54070.5
plot(mean(imagen))

⚠️ ATENCIÓN:

Otros argumentos son: log(), sqrt(), cumsum() o sign(). Algunas funciones, como cumsum(), devuelven más de una capa.

Además, se puede aplicar la función terra::app. El estadístico deseado se incorpora con el argumento fun. Por ejemplo:

promedio_imagen <- app(imagen, fun = mean, na.rm = TRUE)
promedio_imagen
## class       : SpatRaster 
## dimensions  : 1395, 1848, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source(s)   : memory
## name        :    mean 
## min value   :  8146.5 
## max value   : 54070.5
plot(promedio_imagen, main = "promedio de todas las bandas, axes = TRUE")

ANÁLISIS GRÁFICO

La obtención de los valores estadísticos de forma numérica puede ser completada con un análisis gráfico, usando para ello las diferentes herramientas disponibles en R.

Por ejemplo, los histogramas muestran la distribución de frecuencias de los valores de los píxeles de una imagen raster. Con su ayuda se puede cuantificar la frecuencia de niveles digitales, ya que el histograma indica el número de pixeles presentes en cada rango. Además, se puede extraer una primera aproximación de las posibles clases temáticas y definir la calidad de la imagen por medio del contraste: una imagen que genere un histograma muy estrecho corresponde a una imagen con un pobre contraste visual.

hist(promedio_imagen, 
    col = "wheat", 
    main = "Histograma de valores medios de todas las bandas")
## Warning: [hist] a sample of 39% of the cells was used (of which 1% was NA)

El rango de los valores que puede tener un pixel se muestra en el eje X, mientras que su frecuencia aparece en el eje Y.

Lo más usual es representar simultáneamente todas las bandas disponibles en el objeto ráster.

par(mfrow=c(2,2))
hist(imagen[[1]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="darkblue", xlab=" Banda blue")
## Warning: [hist] a sample of 39% of the cells was used (of which 1% was NA)
## Warning in plot.window(xlim, ylim, "", ...): "maxpixels" is not a graphical
## parameter
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "maxpixels" is not a graphical parameter
## Warning in axis(1, ...): "maxpixels" is not a graphical parameter
## Warning in axis(2, at = yt, ...): "maxpixels" is not a graphical parameter
hist(imagen[[2]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="green",xlab=" Banda green")
## Warning: [hist] a sample of 39% of the cells was used (of which 1% was NA)
## Warning in plot.window(xlim, ylim, "", ...): "maxpixels" is not a graphical
## parameter
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "maxpixels" is not a graphical parameter
## Warning in axis(1, ...): "maxpixels" is not a graphical parameter
## Warning in axis(2, at = yt, ...): "maxpixels" is not a graphical parameter
hist(imagen[[3]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="red", xlab="Banda red")
## Warning: [hist] a sample of 39% of the cells was used (of which 1% was NA)
## Warning in plot.window(xlim, ylim, "", ...): "maxpixels" is not a graphical
## parameter
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "maxpixels" is not a graphical parameter
## Warning in axis(1, ...): "maxpixels" is not a graphical parameter
## Warning in axis(2, at = yt, ...): "maxpixels" is not a graphical parameter
hist(imagen[[4]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="darkred", xlab="Banda NIR")
## Warning: [hist] a sample of 39% of the cells was used (of which 1% was NA)
## Warning in plot.window(xlim, ylim, "", ...): "maxpixels" is not a graphical
## parameter
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "maxpixels" is not a graphical parameter
## Warning in axis(1, ...): "maxpixels" is not a graphical parameter
## Warning in axis(2, at = yt, ...): "maxpixels" is not a graphical parameter

par(mfrow=c(1,1))

Una interesante alternativa para la visualización de objetos ráster es el paquete rasterVis.

library(rasterVis)
## Warning: package 'rasterVis' was built under R version 4.3.2
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.2

La función rasterVis::histogram() dibuja los histogramas de todas al bandas simultáneamente.

histogram(imagen, 
          main = "Histograma de todas las bandas", 
          xlab = "Valores de los píxeles", 
          ylab = "Frecuencia", 
          col = "wheat")

Un gráfico de densidad muestra la distribución de los valores de una imagen raster con intervalos continuos (a diferencia del histograma), y es muy útil para determinar la forma de la distribución de esa imagen. Rl pico del gráfico de densidad muestra el intervalo de los datos con una mayor frecuencia (la moda). Para ello se usa la función rastervis::densityplotdel paquete rasterVis

densityplot(imagen,
            xlab = "Valores",
            main = "Gráfico de densidad de cada banda")

El gráfico de violín es una combinación del gráfico de caja y bigotes y del gráfico de densidad rotado. El primero muestra la distribución de los datos a partir de los cuartiles.

● La caja incluye todos los datos entre los cuartiles superior e inferior.

● Los bigotes representan la variabilidad de los datos más alllá de los cuartiles superior e inferior.

Además, el gráfico de caja y bigotes incluye otros estadísticos como la media, el rango intercuartílico, etc..

El gráfico de densidad muestra la distribución de los datos en términos de densidad de píxeles con respecto a sus valores. Por lo tanto, un gráfico de violín muestra a la vez un resumen de los datos y la distribución completa de las variables. Para realizarlo sobre datos raster acudiremos a la funcion bwplotdel paquete rasterVis.

bwplot(imagen, 
       main = "Gráfico de violín")

CUANTIFICACIÓN DE LA RELACIÓN ENTRE BANDAS

La aplicación de otros estadísticos permite conocer la relación entre las bandas de una imagen. Los estadísticos multibanda más habituales son la covarianza y la correlación, normalmente acompañados de un diagrama de dispersión por pares.

Covarianza

Mide la variación conjunta de los niveles digitales de los pixeles de dos bandas alrededor de una media común. Es esencialmente una medida de la varianza entre dos variables. Toma valores positivos cuando ambas variables tienden a aumentar o disminuir a la vez y valores negativos cuando una variable tiene a incrementar mientras la otra disminuye o viceversa. Sin embargo, no evalúa la dependencia entre variables y se mide en las mismas unidades que las que aparecen medidas dichas variables: la varianza puede tomar cualquier valor positivo o negativo.

  • Covarianza positiva: indica que dos variables tienden a moverse en la misma dirección.

  • Covarianza negativa: revela que dos variables tienden a moverse en direcciones inversas.

La fórmula de cálculo es la siguiente

\[cov_{x,y} = \frac{\sum\limits_{i=1}^{n}{(x_i-\overline{x}) \cdot (y_i-\overline{y})} }{n-1}\]

donde

\(X_i\) es el nivel digital del pixel en la banda \(X\),

\(Y_i\): Nivel digital del pixel en la banda \(Y\),

\(\bar{x}\) la media aritmética de la banda \(X\),

\(\bar{y}\)y la media aritmética de la banda \(Y\) y

\(m\) el número total de píxeles en la imagen

Para el cálculo de la matriz de covarianza entre las bandas de una misma imagen es necesario transformar cada una de ellas en una matriz.

cov <- cov(as.matrix(imagen), 
           use = "complete.obs")                                                # Utiliza los pares completos de ambas variables

round(cov*1000,2)
##             blue      green        red        NIR
## blue  1275728032 1555661556 1941097249 1398264024
## green 1555661556 2076142997 2753069621 2236621479
## red   1941097249 2753069621 3968943740 3294359200
## NIR   1398264024 2236621479 3294359200 5105280315

Coeficiente de correlación

También mide la relación entre los niveles digitales de dos bandas, pero toma valores entre entre -1 y 1, siendo \(r_{xy}=1\) una relación perfecta entre las dos bandas, \(r_{xy}=- 1\) una relación inversa perfecta entre las dos bandas y \(r_{xy}=0\) no existe relación entre las dos bandas.

\[\rho = \frac{\text{cov}(X,Y)}{\sigma_x \sigma_y}\]

siendo

\({cov}(X,Y)\) la covarianza de las bandas X y Y,

\(sigma_x\) la desviación estándar de la banda \(X\) y

\(sigma_y\) la desviación estándar de la banda \(Y\).

Para el cálculo del coeficiente de correlación también debe transformarse la imagen ráster en una matriz de datos

cor <- cor(as.matrix(imagen), method="spearman",
           use = "complete.obs")                                                # Utiliza los pares completos de ambas variables
cor
##            blue     green       red       NIR
## blue  1.0000000 0.9551156 0.8714753 0.5585544
## green 0.9551156 1.0000000 0.9639095 0.6937604
## red   0.8714753 0.9639095 1.0000000 0.7254449
## NIR   0.5585544 0.6937604 0.7254449 1.0000000

El coeficiente de correlación puede calcularse también transformando el fichero ráster se transforme en un dataframe.

x <- as.data.frame(na.omit(imagen))
cor(x, method="spearman", use='complete.obs')                                   # Todas las bandas simultáneamente
##            blue     green       red       NIR
## blue  1.0000000 0.9551156 0.8714753 0.5585544
## green 0.9551156 1.0000000 0.9639095 0.6937604
## red   0.8714753 0.9639095 1.0000000 0.7254449
## NIR   0.5585544 0.6937604 0.7254449 1.0000000

Esto tiene ciertas ventajas. En primer lugar, permite reducir el número de bandas a calcular si fuera necesario.

cor(x$blue, x$NIR, method="spearman", use='complete.obs')                       # Sólo dos bandas
## [1] 0.5585544

Existen otras opciones para la representación gráfica de la relación entre capas de una imagen. Una de ellas es el gráfico de correlaciones o correlograma.

library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.3
## corrplot 0.95 loaded
matriz_correlacion <- round(cor(x, method="spearman", use='complete.obs'), 2)
corrplot(matriz_correlacion, 
         type = "upper", 
         order = "hclust", 
         tl.col = "black", 
         tl.srt = 45)

tmp1 <- as.matrix(spatSample(imagen,500000))
library(PerformanceAnalytics)
## Warning: package 'PerformanceAnalytics' was built under R version 4.3.3
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.3.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
chart.Correlation(tmp1, 
                  histogram=TRUE, 
                  pch=19)

Un diagrama de dispersión por pares constituye una buena alternativa a los valores numéricos. La comparacion entre los valores espectrales de dos bandas de una misma imagen raster puede realizarse mediante la función terra::pairs. Esta función devuelve un gráfico de dispersión entre las dos bandas especificadas (que muestra el valor de un determinado pixel en ambas bandas y el grado de relacion entre ambas), sus histogramas y el valor del coeficiente de correlación entre ambas.

pairs(imagen[[1:4]], 
      hist=TRUE,                                # Muestra el histograma
      cor=TRUE,                                 # Muestra el coeficiente de correlación
      main = "Bandas Red versus NIR",           # Título
      use="pairwise.complete.obs",              # Correlaciones calculadas usando todas las parejas de variables (número de parejas diferente)
      maxpixels=100000) 

La funcion rasterVis::splom[https://rdrr.io/cran/rasterVis/man/splom-methods.html] del paquete rasterVis también ofrece el gráfico de dispersión y la curva de densidad de todas las bandas al mismo tiempo.

splom(raster::brick(imagen),                                                    # Atención, rasterVis sólo funciona con el paquete `raster`.
      main = "Gráfico de dispersion con todas las bandas")

📝 ACTIVIDAD DE EVALUACIÓN OPERACIONES GLOBALES:

Los materiales para esta actividad están incluidos en el mismo fichero que incluye las actividades locales, zonales y globales.

ACTIVIDAD 1. En el presente apartado trabajaremos con varios ficheros ráster, que son los siguientes:

  • lst.tif

  • ndvi.tif

  • savi.tif

  • evi.tif

  • Importa cada uno de estos ficheros por separado, creando luego un objeto ráster multibanda que agrupe a todos los ficheros.

  • Multiplica los valores originales del objeto ráster ndvi por un factor, en este caso, por 1000, y represéntalo gráficamente

  • Suma los índices NDVI y EVI, multiplica el resultado por 1000 y divide todo el conjunto entre 2 para obtener un nuevo índice.

Si desea realizar alguna operación matemática que involucre todos los píxeles de una capa ráster, es conveniente definir una función. En el siguiente ejemplo, se aplica álgebra de mapas para cambiar la escala de los valores del índice EVI 0 a 100.

reescalado <- function(x) {
  if (max(x, na.rm = TRUE) == min(x, na.rm = TRUE)) {
    return(rep(0, length(x)))  # or another appropriate constant
  }
  y <- (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)) * 100
  return(y)
}
  • A continuación, aplica esta función al índice EVI y dibuja en un lienzo tanto el mapa correspondiente al índice EVI original como el reescalado.

  • Calcula el coeficiente de correlación de Pearson entre la banda lst y la banda ndvi

  • Crea un correlograma calculado a partir del coeficiente de correlación de Pearson aplicado a todas las bandas del objeto indices.

ACTIVIDAD 2. En esta actividad se trabajará con un único índice medido en diferentes momentos del tiempo. Los ficheros son los siguientes:

  • NBR_2022-01-05.tif

  • NBR_2022-06-14.tif

  • NBR_2022-07-16.tif

  • NBR_2022-08-18.tif

  • NBR_2022-09-18.tif

  • NBR_2022-11-05.tif

  • Crea una objeto ráster multibanda que se llamará nbr, conteniendo todos esos ficheros mediante una lista. El patrón de nombres que puedes utilizar es el siguiente: pattern = “^NBR.*\.tif$”

  • Cambia los nombres de las variables sustituyéndolos por el nombre del mes en el que se tomó la imagen

  • Calcula el valor medio del índice NBR entre todas las imágenes y represéntalo en forma de mapa.

  • Divide los valores del índice NBR entre 10000

  • ¿En qué mes se registraron los valores más altos de NBR?

  • ¿En qué mes se registrató una mayor dispersión de los valores de NBR?

  • Utiliza la función hist() para representar la distribución de valores del mes con mayor índice NBR

  • Utiliza la función rasterVis::bwplot() para mostrar la distribución del índice NBR en todos los meses .

  • Calcula una matriz con los valores del coeficiente de correlación de Spearman entre las diversas fechas en las que se tomaron las imágenes. ¿Cúal es la pareja de meses con el valor más alto del coeficiente de correlación? ¿y la pareja con el coeficiente más bajo?

  • Elabora un diagrama de dispersión con la función pairs() entre el índice NBR de los meses enero, junio, septiembre y noviembre.

  • Elabora un correlograma que muestre la relación enter todos los meses.