💡 MATERIALES PARA LA ACTIVIDAD:

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

► Imágenes satélite.

► script_laboratorio_5_operaciones_locales

INTRODUCCIÓN

Son el grupo de operaciones más simples de álgebra de mapas. En esencia, consisten en la aplicación de algún tipo de función matemática a dos o más ráster, utilizando los valores de sus píxeles. Esta función matemática se aplica píxel a píxel, sobre píxeles que coinciden en ubicación y dimensiones.

Ejemplo de unas operación local
Ejemplo de unas operación local

Un caso típico de utilización del álgebra de mapas en el mundo es la estimación de la erosión potencial mediante el modelo USLE. Este modelo utiliza la siguiente ecuación:

\[Ex,y = Rx,y ∗ Kx,y ∗ Lx,y ∗ Sx,y ∗ Cx,y ∗ Px,y\]

donde

Si se dispone de una capa raster para cada uno de estos factores resulta muy simple obtener una estimación de la erosión

ACTIVIDADES INICIALES

Activación de paquetes de R

Para esta actividad sólo es necesario activar el paquete terra.

library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21

Preparación de los datos

Deberá establecerse el directorio de trabajo.

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

Desde ese directorio, se importarán los dos ficheros en formato GeoTIFF, solicitándo a R que muestre sus atributos. Uno de ellos es un fichero

L8 <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_locales/L8_NDVI.tif")
L8
## class       : SpatRaster 
## dimensions  : 1500, 1350, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 229995, 270495, 3964995, 4009995  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 11N (EPSG:32611) 
## source      : L8_NDVI.tif 
## name        :    L8_NDVI 
## min value   : -0.2774381 
## max value   :  0.5495766

A continuación se representará este objeto ráster

plot(L8)

El otro fichero es proviene de la plataforma Sentinel 2.

sentinel2a <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_locales/Sentinel2A.tif")
sentinel2a
## class       : SpatRaster 
## dimensions  : 1242, 1061, 2  (nrow, ncol, nlyr)
## resolution  : 9.997187, 9.997187  (x, y)
## extent      : 570465.7, 581072.7, 6368052, 6380468  (xmin, xmax, ymin, ymax)
## coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154) 
## source      : Sentinel2A.tif 
## names       : T31TCK_20211012T105011_B04_10m, T31TCK_20211012T105011_B08_10m 
## min values  :                              1,                              4 
## max values  :                          16112,                          11273
plot(sentinel2a)

TIPOS DE OPERACIONES LOCALES

Las operaciones locales hacen referencia a cálculos realizados independientemente sobre una célula, a partir de una o numerosas capas. En el caso del análisis de imágenes de satélite, podemos mencionar las siguientes.

Sustitución de valores

Para poder trabajar sobre el fichero imagen duplicaremos este último con el nombre de borrame.

borrame1 <- L8

La imagen muestra el índice NDVI correspondiente al lago Tulare en California, EEUU. Valores inferiores a -0.1 podrían corresponder a agua. A continuación, los valores por encima de ese valor van a ser reemplazados por NA con el fin de identificar los límites del lago.

borrame1[borrame1[[1]] > -0.1]  <- NA
plot(borrame1)

Además, también podría considerarse reemplazar los anteriores NA por 0.

borrame1[is.na(borrame1)] <- 0
plot(borrame1)

Una alternativa es la función terra::clamp() que restringe los valores a un mínimo y a un máximo. Es decir, todos los valores por debajo de un valor umbral inferior y por encima del valor umbral superior se convierten en NA o, si el argumento values=TRUE, se convierten en el valor umbral

borrame2 <- L8
borrame2 <- clamp(borrame2, 
                  -Inf, -0.1, 
                  values= FALSE)                       # Si FALSE los valores fuera se convierten en NA
plot(borrame2)

Operaciones sobre cada pixel

Por ejemplo, puede multiplicarse el valor de cada píxels por 10

borrame2 <-  L8 * 100
plot(borrame2)

También podría sustraerse el valor mediano del índice de cada una de los píxeles

borrame3 <-  borrame2 - global(borrame2, 
                               fun = median, 
                               na.rm = TRUE)[[1]]
plot(borrame3)

Reclasificación

La reclasificación de valores ráster puede utilizarse para discretizar datos cuantitativos o para categorizar datos cualitativos. La función terra::classify() realiza esta clasificación. En su forma más sencilla, requiere la implementación de una matriz con 3 columnas: las dos primeras con los límites inferior y superior de cada clase, y la tercera con la etiqueta (valor) que adquiere en el nuevo ráster.

⚠️ Atención:

Si se desea una clasificación con los intervalos cerrados a la derecha y abiertos a la izquierda deben incluirse los argumentos right=TRUE e include.lowest=TRUE.

▶ 1-5 -> 1 : todos los casos entre 1.01 y 5.

▶ 5-10 -> 2 : todos los casos entre 5.01 y 10.

▶ 10-15 -> 3 : todos los casos entre 10.01 y 15.

El primer paso para reclasificar consiste en crear una matriz de datos con los límites inferior y superior de cada categoría, y con el valor que sustituirá a los valores originales.

matriz_clas <- matrix(c(-Inf, -0.1, 1,
                         -0.1, 0.4, 2,
                          0.4, Inf, 3),
                    ncol = 3, 
                    byrow = TRUE)
matriz_clas
##      [,1] [,2] [,3]
## [1,] -Inf -0.1    1
## [2,] -0.1  0.4    2
## [3,]  0.4  Inf    3

A continuación se aplica la función terra::classify() al ráster original, aplicando la matriz anterior.

L8_clases <- classify(L8, 
                      rcl = matriz_clas, 
                      right=TRUE, 
                      include.lowest=TRUE)
names(L8_clases) <- "usos suelo"

plot(L8_clases, 
     type = "classes",
     levels = c("Agua","Matorral","Cultivado"),
     col = c("blue", "brown", "green"))

Operaciones combinando varias capas

Es posible calcular un valor de celda a partir de diferentes valores almacenados en varias capas de un objeto SpatRaster. El ejemplo más común es probablemente el cálculo de un índice multiespectral, por ejemplo el NDVI. La imagen que sirve de ejemplo, obtenida por el satélite Sentinel-2 y extraida de Copernicus Open Access Hub, contiene dos bandas espectrales, con una resolución de 10m: la banda roja y el infrarrojo cercano. El índice NDVI se calculará en todas y cada una de los píxeles, a partir de las capas de datos ráster de una imagen de satélite multiespectral.

plot(sentinel2a)

names(sentinel2a) <- c("B04", "B08")

A partir de ambas bandas, se puede obtener otro objeto ráster cuyos valores provienen de la aplicación de una fórmula a los datos originales:

sentinel2_NDVI <- (sentinel2a[["B08"]] - sentinel2a[["B04"]]) / (sentinel2a[["B08"]] + sentinel2a[["B04"]])
plot(sentinel2_NDVI)

📝 ACTIVIDAD DE EVALUACIÓN OPERACIONES LOCALES:

ACTIVIDAD DE EVALUACIÓN 1

► En esta actividad se utilizará una imagen de satélite denominada RT_T36MZE_20160108T075312.tif.

  1. Es una imagen del satélite Sentinel 2 con todas las bandas que poseen una resolución de 20 m. Transfórmalo en un objeto ráster multibanda denominado raster y cambia el nombre de todas las bandas.

  2. Elabora dos objetos ráster, denominados indice1 e indice2, resultantes de aplicar las siguientes fórmulas al objeto ráster. No es necesario que se extraigan la bandas.

\[indice1 = (B8 – B4) / (B8 + B4)\].

\[indice2 = (B3 - B8) / (B3 + B8)\]

  1. Representa gráficamente ambos ráster en el mismo lienzo.

  2. A partir del objeto ráster indice1 y aplicando las categorías y umbrales siguientes.

Intervalos Etiqueta Uso de suelo
< a 0.20 1 Agua
0.20 a 0.30 2 Urbano
0.30 a 0.40 3 Cultivos
> 0.40 4 Bosque

  • Clasifica los usos de suelo

  • Representa gráficamente esas categorías

  • Calcula el porcentaje de la superficie que representa cada categoría en forma de tabla. Verifica primero que es ráster categórico. Si no lo fuera, conviértelo.

  • Elabora los mapas que reproducen las diferentes categorías con la función terra::segregate().

  • Aplica la función terra::clamp() para resaltar el lago utilizando el valor de -0.05 como límite superior.

  • Finalmente, representa de manera conjunta sendos mapas con el objeto ráster indice1 y con el objeto ráster lago. Añade título y etiquetas en los ejes X e Y a conveniencia.

ACTIVIDAD DE EVALUACIÓN 2

► En esta actividad se utilizará una imagen con la temperatura de la superficie terrestre denominada LC09_L2SP_201032_20220718_20220720_02_T1_ST_B10.tif, así como un fichero vectorial denominado limite_madrid.json, ambos incluidos en el fichero que se puede descargar desde internet.

  • Crea, a partir del fichero ráster, un objeto SpatRaster con el nombre lst y un objeto SpatVector con el nombre limite.

  • Utiliza el objeto SpatVector para recortar el objeto SpatRaster.

  • Calcula la temperatura de la superficie terrestre en ºC aplicando la siguiente fórmula:

\[LST = (DN * 0.00341802 + 149.0) - 273.15\]

  • Representa gráficamente la temperatura de superficie utilizando como paleta de colores la función colorRampPalett(). Construye esta paleta con los colores “lightblue”, “lightgreen”, “yellow”, “orange”, “red”, y 100 intervalos de colores. El título del mapa deberá ser Temperatura de la superficie terrestre.

  • Graba el objeto ráster con la temperatura de la superficie terrestre

  • Elabora un histograma con la temperatura de superficie

  • Calcula los valores promedio, máximo, mínimo y los cuartiles 1º y 3º de la temperatura de superficie.

  • Calcula la diferencia entre el valor de cada pixel y la temperatura media de todo el conjunto. Represéntala gráficamente utilizando una paleta con los colores “darkblue”, “lightblue”, “white”, “orange”, “red”, y 100 intervalos de colores

  • A continuación, identifica zonas muy cálidas o muy frescas (“Hot/Cold spots”) transformando los valores originales en puntuaciones z y seleccionando sólo aquellos puntos cuyo valor sea <-1.96 (áreas muy cálidas) o > +1.96 (áreas muy frías). Esto último tendrás que hacerlo reclasificando el objeto ráster. Finamente, crea un gráfico sólo con ellas. Recuerda que las puntuaciones Z se calculan de la siguiente manera:

\[Z = \frac{X – μ}{σ}\] siendo

\(X\) el valor de cada pixel

\(μ\) el promedio de los valores de LST en todo el objeto ráster.

\(σ\) la desviación típica o estandar de los valores de LST en todo el objeto ráster.

  • Finalmente, graba en formato *.tif tanto el objeto SpatRaster con los puntos cálidos y fríos.