💡 OBJETIVOS:

💡 MATERIALES PARA LA ACTIVIDAD:

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

Imágenes satélite. Es la misma imagen multibanda utilizada en la actividad previa.

script_laboratorio_6_density_slicing

INTRODUCCIÓN

Bajo el término “Density slicing” se engloba un procedimiento para establecer unos umbrales sobre los valores de un índice espectral. En este procedimiento se establecen una serie de valores-umbral especificados por el usuario. Esto resulta en la división de los valores de una imagen en una serie de categorías excluyentes, que pueden asociarse a un tipo de uso de suelo, o superficie terrestre (por ejemplo, los valores de NDVI pueden usarse para definir áreas sin vegetación, pradera, matorral o bosque).

Umbralizar un objeto ráster es un tipo especial de cuantificación que segmenta los valores de los píxeles en dos clases, dependiendo de un valor-umbral muy concreto. La finalidad es identificar qué píxeles conforman los fenómenos espaciales buscados, asignándolos a cierto grupo, llamado comúnmente “segmento”, y qué píxeles son sólo el entorno de estos objetos. El objeto a segmentar está compuesta por valores numéricos, por lo que la pertenencia de un píxel a un cierto segmento se decide mediante la comparación de su valor con el valor umbral.

En ambos casos, se pasa de una imagen con valores cuantitaticos continuos a una imagen con valores cualitativos discretos que se puede convertir en un ráster categórico. Este último puede ser luego usado para un análisis detallado de cada una de las clases de manera separada. Por ejemplo, si se quisiera calcular la superficie total de un embalse en dos momentos del tiempo diferentes para analizar su evolución, esto se podría estimar mejor si los píxeles que no son agua desaparecieses, quedando sólo aquellos que representan el agua.

La clave de este procedimiento es la elección de los valores-umbral, que puede realizarse de manera manual, o apoyándose en ciertas herramientas estadísticas. Siguiendo con el ejemplo, un valor del índice NDWI mayor que 0 representaría el valor umbral; a todos los píxeles con un valor superior a ese umbral se les puede asignar el valor 1, mientras que los inferiores pueden pasar a NA. A partir de las imágenes de salida se podría estimar fácilmente la superficie de las masas de agua y su variación en el tiempo.

PREPARACIÓN DE LOS DATOS

Los paquetes básicos para realizar este ejercicio son los siguientes:

library(terra)                             # Importación y transformación de imágenes ráster
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21
library(ggplot2)                           # Visualización de datos ráster con RStoolbox
## Warning: package 'ggplot2' was built under R version 4.3.3
library(sf)                                # Manejo de ficheros en formato vectorial.
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(RStoolbox)                        # Transformación y visualización de datos ráster.
## Warning: package 'RStoolbox' was built under R version 4.3.3
## This is version 1.0.0 of RStoolbox
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
## Warning: package 'readr' was built under R version 4.3.2
## Warning: package 'purrr' was built under R version 4.3.2
## Warning: package 'dplyr' was built under R version 4.3.2
## Warning: package 'stringr' was built under R version 4.3.2
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

En el presente apartado trabajaremos con el objeto ndvi, que contiene los valores del índice homónimo, y el objeto imagen, con el que se ha trabajado en unidades anteriores. Para ello, se establece la carpeta de trabajo

setwd("D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/datos/")      

A continuación se importan ambos objetos ráster

imagen <- rast("D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/datos/imagen.tif")
names(imagen) <- c("coastal aerosol", "blue", "green", "red", "NIR", "SWIR1", "SWIR2")

A continuación, se calculará el índice NDVI

ndvi <- (imagen[["NIR"]] - imagen[["red"]]) / (imagen[["NIR"]] + imagen[["red"]])
ndvi
## class       : SpatRaster 
## dimensions  : 1245, 1497, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source(s)   : memory
## varname     : imagen 
## name        :       NIR 
## min value   : -0.959918 
## max value   :  0.883263

Se realiza una primera representación gráfica del índice NDVI y se compara con la imagen en falso color.

par(mfrow = c(1,2))
#imagenFCC <- imagen[[c(5,4,3)]]
plotRGB(imagen, 
        r=5, g=4, b=3, 
        stretch = "lin",
        main="Imagen en falso color", 
        mar=c(3.5, 3.5, 2.5, 2.5))

plot(ndvi, main = "Índice NDVI")

par(mfrow = c(1,1))

CÁLCULO DE UMBRALES

Método manual

Un enfoque tradicional es utiliza el valor 0 como umbral. Esto significa que los pixeles que tengan un valor superior a 0 se clasificarán como perteneciente a una determinada clase, por ejemplo, vegetación o agua, asignando NA a los píxeles con valores inferiores.

En la función terra:: classify()[https://rspatial.github.io/terra/reference/classify.html] los argumentos principales son la imagen ráster de entrada y la matriz que sirve de guía para la reclasificación. Esta matriz contiene tres columnas.

  • La primera y la segunda columnas contienen los umbrales (límites inferior y superior de un rango).

  • La tercera columna contiene el valor que se asignará a ese rango.

A continuación se muestra un ejemplo:

umbral_veg1 <- c(-Inf, 0, NA,  
                 0, Inf, 1)

matriz_clas_veg1 <- matrix(umbral_veg1, 
                           ncol=3, 
                           byrow=TRUE)
matriz_clas_veg1
##      [,1] [,2] [,3]
## [1,] -Inf    0   NA
## [2,]    0  Inf    1

A continuación se clasifica el objeto ndvi en las categorías señaladas, y se representa gráficamente.

imagen_veg1 <- classify(ndvi, 
                        matriz_clas_veg1)
plot(imagen_veg1)

A pesar de su sencillez, los resultados puede introducir errores debido a cambios espacio temporales en el brillo y contraste de las imágenes.

📝 ACTIVIDAD DE EVALUACIÓN 1:

A continuación se realizará un ejercicio en el que se va a utilizar el objeto ráster ndvi generado anteriormente. Si consideramos que los píxeles con valores superiores a 0.4 probablemente correspondan a vegetación, utilícese este umbral para determinar la superficies con algún tipo de vegetación.

► Crea la matriz con las clases por debajo y por encima.

► Clasifica el ráster según los umbrales definidos

Hasta este momento se ha clasificado el objeto ráster en 2 categorías, por encima y por debajo de un determinado valor. Sin embargo, también se pueden obtener categorías intermedias. Por ejemplo, si desea retener píxeles con valores NDVI que van desde 0 a 0,4, probablemente representando suelos con escasa vegetación se debe asignar el valor ‘1’ a todos los píxeles dentro de ese rango de valores y NA a los píxeles por encima y por debajo de esos límites

umbral_suelo <- c(-Inf, 0, NA,
                0, 0.4, 1, 
                0.4, Inf, NA)

matriz_clas_suelo <- matrix(umbral_suelo, 
                            ncol=3, 
                            byrow=TRUE)

imagen_suelo <- classify(ndvi, 
                         matriz_clas_suelo)
plot(imagen_suelo, col = "tan")

Esta última imagen puede superponerse a una imagen en color natural para mejorar la interpretación del tipo de superficie

plotRGB(imagen, r=4, g=3, b=2, axes=TRUE, stretch="lin")
plot(imagen_suelo, main = 'Mapa suelo desnudo', col = "tan", add=TRUE, legend=FALSE)

Umbrales según histograma

Una manera de definir uno (o varios) umbral(es) cuando se carece de valores predeterminados es explorar la distribución de esos valores en el objeto ráster mediante la función hist(), que produce un histograma a partir del cual se pueden identificar tanto valores erróneos, valores extremos pero también cambios en la distribución de frecuencias que pudieran sugerir discontinuidades (umbrales).

hist(ndvi,
     main = "Distribución de valores del NDVI según píxeles",
     xlab = "NDVI",
     ylab= "Frecuencia",
     col = "green",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
## Warning: [hist] a sample of 54% of the cells was used
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

En este caso, obsérvese que hay un corte que coincide con el valor 0.15, que podría representar el tránsito entre superficies no vegetadas y las vegetadas. A su vez, a partir de -0.1 se observa otra discontinuidad, que podría corrresponder a la aparición de agua. Por la derecha de la distribución la frecuencia de valores cae rápidamente por encima de 0.35, que podría representar la aparición de bosques.

Método de Otsu

Este procedimiento manual puede ser refinado utilizando el método de discretización de clases de Otsu, llamado así en honor a Nobuyuki Otsu en 1979 (Otsu, 1979), utiliza técnicas estadísticas, para resolver el problema. En concreto, utiliza la varianza, una medida de la dispersión de valores. Su método calcula el valor umbral de forma que la dispersión dentro de cada segmento sea lo más pequeña posible, pero al mismo tiempo la dispersión sea lo más alta posible entre segmentos diferentes. Para ello se calcula el cociente entre ambas varianzas y se busca un valor umbral para el que este cociente sea máximo.

A continuación se mostrará una aplicación práctica de esta metodología, aplicándose a un raster en el que se ha calculado el índice NDWI: los valores negativos corresponden a terreno continental y los positivos a agua.

\[NDWI = (G - NIR)/(G + SWIR)\]

El primer procedimiento consiste en calcular el citado índice

ndwi <- (imagen[["green"]] - imagen[["NIR"]]) / (imagen[["green"]] + imagen[["NIR"]])
plot(ndwi, 
     col=rev(topo.colors(100)), main="NDWI")

Para destacar y extraer los cuerpos de agua de mayor tamaño utilizando el método de Otsu se recurre a una función del paquete EBImage. EBImage proporciona una funcionalidad de propósito general para el procesamiento y análisis de imágenes (Pau et al. [2010]) y es un paquete creado y mantenido por la iniciativa Bioconductor. Para la mantención e instalación de sus propios paquetes instalamos en primer lugar el paquete para gestión, BiocManager (Morgan [2021]).

#install.packages("BiocManager")
#library(BiocManager)

A continuación, se puede instalar el paquete deseado (por ejemplo EBImage) utilizando la función install() del paquete BiocManager.

#if (!require("BiocManager", quietly = TRUE))
#install.packages("BiocManager")

# BiocManager::install("EBImage", force = TRUE)

Y luego cargar en la sesión actual:

library(EBImage)
## 
## Attaching package: 'EBImage'
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:terra':
## 
##     flip, rotate, thresh, watershed

Las funciones del paquete requieren la información de nuestro modelo ráster en formato de array o matriz. La función terra::as.array() [https://rspatial.github.io/terra/reference/coerce.html] que toma un raster y devuelve como array.

array <- as.array(ndwi)
str(array)
##  num [1:1245, 1:1497, 1] -0.465 -0.469 -0.497 -0.545 -0.564 ...

La función EBIimage::otsu() devuelve un valor umbral o de corte para binarizar la imagen

corte <- otsu(array, 
              range = c(-1,1))                              # Número de niveles de gris
corte
## [1] -0.08984375

Aplicamos el valor de corte, pero antes duplicaremos la imagen original para trabajar sobre ella.

agua <- ndwi                                           
agua[agua >= corte] <- 1
agua[agua < corte] <- NA

El resultado del proceso anterior es una imagen binaria, en el que el valor 1 para los píxeles asignados como cuerpo de agua y el resto NA.

plot(agua, 
     col= "blue", 
     axes=F,
     box=T, 
     legend=F)

El método de Otsu es adecuado para la clasificación de índices espectrales sencillos, como NDVI (para separar las superficies con vegetación de las que no lo tienen), NDBI (extraer las zonas construidas) o el NDWI (extraer las superficies con agua). Como ventajas adicionales se encuentra el que es un procedimiento automático, sin intervención de la subjetividad del usuario, que funciona especialmente bien cuando los valores de esos índices muestran una distribución bimodal (histograma con dos clases dominantes), pero, por contra, no suele funcionar bien si la distribución es unimodal o presenta diversas clases dominantes.

POST PROCESAMIENTO DE UN RÁSTER BINARIO

A continuación vamos a revisar alguna operaciones adicionales sobre imágenes binarias de este tipo. Por ejemplo, este procedimiento facilita la identificación y delimitación espacial de algunos tipos de superficies (pe. una zona deforestada, un lago o pantano, un glaciar); comparando los resultados de este procedimiento en diferentes momentos de tiempo es posible seguir su evolución.

Parches

Cada grupo de píxeles con valor 1 adyacentes o en contacto forman un “parche”. La identifiación y agrupación de los píxeles individuales en conglomerados más grandes se realiza con la función terra::patches(). Esta función detecta grupos de píxeles conectados y proporciona una identificación única. Cuando todos los vecinos correspondientes al tipo pertenecen al mismo parche, el píxel central se asigna a ese id. NA y cero se utilizan como valores de fondo (es decir, estos valores se utilizan para separar grupos).

masas_agua <- patches(agua, 
                directions = 8,          # La definición de las celdas adjacentes se puede usar el kernel de vecinos tipo `reina` o `torre`. 
                allowGaps = FALSE,
                zeroAsNA=TRUE)           # Argumento `allowGaps=TRUE`: los ID no son consecutivos (por ejemplo, puede haber parches con los IDs 1, 2, 3 and 5, pero no el 4). Si se establece como FALSE, los números serán consecutivos.

plot(masas_agua, col = "blue")

A continuación se observa el gran número de agrupaciones o parques, debido a que muchos tienen una superficie que corresponde a un único pixel. El argumento range de la función terra::global puede servir para identificar el número de agrupaciones o parches.

global(masas_agua, 
       range, 
       na.rm = TRUE)
##         X1   X2
## patches  1 1818

Puede observarse que se han detectado 1818 “parches” con píxeles identificados como agua. Su numeración corresponde a un orden aleatorio de definición y no tiene relación con la ubicación geográfica, ni con su tamaño.

De la misma manera, cada uno de los elementos que conforman el objeto masas_agua puede dibujarse con una paleta secuencial que identifica cada uno de los individuos.

library(colorspace)
## Warning: package 'colorspace' was built under R version 4.3.3
## 
## Attaching package: 'colorspace'
## The following object is masked from 'package:terra':
## 
##     RGB
paleta <- qualitative_hcl(n = 11, 
                         palette = "Dark 3", 
                         rev = TRUE) 
plot(masas_agua, 
    col=colorRampPalette(paleta)(64), 
    axes = F, 
    box = T, 
    legend=F)

Si quisiéramos conocer la superficie cubierta por cada parche se utiliza la función terra::cellSice()[https://rspatial.github.io/terra/reference/cellSize.html]. Además, se puede realizar el cálculo de la superficie total de cada parche mediante la función terra::zonal[https://rspatial.github.io/terra/reference/zonal.html]. También se orde

superficie <- cellSize(masas_agua, unit="km") %>% 
  zonal(masas_agua, sum) %>%
  arrange(desc(area))

head(superficie) %>% round(2)
##   patches   area
## 1       6 163.30
## 2    1647   1.13
## 3    1805   1.05
## 4     411   0.98
## 5     752   0.72
## 6     383   0.59

Como se observa, hay un parche que ocupa la mayor parte de la superficie. Si quisiéramos eliminar todos los parches más pequeños, se procedería de la siguiente manera.

Primero se calcularía la frecuencia de píxeles por parche

tamaño_parches <- freq(masas_agua) %>% 
  arrange(desc(count))

A continuación, estableceríamos el umbral, por ejemplo, aquellas masas de agua con más de 1000 píxeles.

tamaño_minimo <- 1000
keep_ids <- tamaño_parches$value[tamaño_parches$count >= tamaño_minimo]

El procedimiento continua creando una máscara para grandes parches

# Set all other patch IDs to NA
masas_agua_filtradas <- classify(masas_agua, 
                             rcl = matrix(cbind(setdiff(tamaño_parches$value, keep_ids), NA), ncol=2))
plot(masas_agua_filtradas)

Opcional, también se puede enmascara los pequeños parches en el original.

filtered_raster <- mask(r, 
                        filtered_patches, 
                        maskvalues = NA)

A continuación se transformará el fichero ráster masas_agua en un fichero vectorial conformado por polígonos.

masas_agua_poligono <- as.polygons(masas_agua_filtradas) %>% 
  st_as_sf()

plot(st_geometry(masas_agua_poligono))

Se puede superponer ese objeto vectorial sobre la imagen en falso color.

plotRGB(imagen, r= 4, g = 3, b = 2, axes=TRUE, stretch="lin")
plot(st_geometry(masas_agua_poligono), add = TRUE, col = "blue")

Sobre el objeto vectorial se puede calcular la superficie total de cada uno de los polígonos, añadirlo como una nueva variable y representarlo gráficamente (se representa el logarítmo en base 10 para minimizar el escaso número de grandes polígonos).

masas_agua_poligono$superficie <- st_area(masas_agua_poligono)
masas_agua_poligono$perimetro <- st_perimeter(masas_agua_poligono)

inalmente, es posible exportar ese vector en formato shapefile para usar posteriormente como máscara, previa transformación en objeto sf:

st_write(masas_agua_poligono, 
         "D:/G174_OCW/Lab_5_transformaciones_matematicas/datos/masas_agua.shp", 
         append = TRUE)

O como en el punto 3.12 se define el uso del paquete mapview para crear sesiones interactivas de cartografía (Figura 9.43) a la cual agregamos los polígonos recién creados.

library(mapview)
## Warning: package 'mapview' was built under R version 4.3.2
mapview(masas_agua_poligono)
mapview(masas_agua_poligono, 
        map.types = c("OpenTopoMap"), 
        legend=F, 
        alpha.regions=0.35, 
        col.regions="lightblue", 
        color="darkblue",
        lwd=3)

📝 ACTIVIDAD DE EVALUACIÓN 2:

El fichero LC08_016028_20210824.tif es una imagen multibanda correspondiente a la ciudad de Ottawa, en Canadá, tomada por el satélite Landsat 8.

► Carga esta imagen en R y comprueba cuántas bandas tiene.

► Selecciona las 7 primeras bandas, que corresponderían a las bandas características de Landsat 8, y elimina todas las demás.

► Cambia los nombres de las bandas.

► Representa gráficamente una imagen en color real de la ciudad de Ottawa. Explica algunos de los elementos del territorio que pueden observarse en la imagen.

► Calcula el índice NDWI y aplica el procedimiento de clasifición de Otsu.

► Compara los resultados proporcionados por el objeto ráster imagen y los proporcionados por el fichero anterior.

EXTRACCIÓN DE VARIAS CATEGORÍAS MEDIANTE EL MÉTODO DEL DENSITY SLICING.

Un proceso algo más complejo es la definición de sucesivos umbrales, que conducen a la aparición de más clases, técnica que en el ámbito del análisis de imágenes se conoce con el nombre de Density slicing. Este procedimiento divide los valores de los píxeles en diferentes rangos y asigna un valor (y se imprime también un color específico) a cada rango en el ráster de salida.

Volviendo al índice NDVI, definido por valores que van de -1.0 a 1.0, la literatura señala que:

Mediante el “Density slicing” es posible generar diferentes categorías (que se interpretan como diferentes masas de vegetación) extendiendo el procedimiento anterior. Primero se crea la matriz de clases

ndvi_clases <- matrix(c(-1.0,  0.0, 1,  # Sin vegetación
                         0.0,  0.2, 2,   # Suelo desnudo/urbano
                         0.2,  0.4, 3,   # Vegetación dispersa
                         0.4,  0.6, 4,   # Matorral
                         0.6,  1.0, 5),  # Bosque
                        ncol = 3, 
                       byrow = TRUE)

A continuación se clasifica el objeto NDVI y se representa gráficamente

ndvi_clasificado <- classify(ndvi, 
                            rcl = ndvi_clases)

plot(ndvi_clasificado, 
     col = terrain.colors(5), 
     main = "NDVI clasificado")

El objeto ráster obtenido a transformado una variable continua (el índice ndvi) en una variable discreta (las 5 categorías de vegetación). Este ráster con valores discretos es un ráster numérico, que puede ser transformado en un ráster categórico. Para ello, primero es necesario definir los nombres de las categorías mediante un dataframe.

levels <- data.frame(value = 1:5,
                     label = c("Sin vegetación", "Suelo desnudo/urbano", "Vegetación dispersa", "Matorral", "Bosque"))

A continuación se aplicarán estos niveles al objeto ndvi_clasificado.

ndvi_factor <- as.factor(ndvi_clasificado)
levels(ndvi_factor) <- levels

Finalmente esto se puede representar gráficamente.

colores_ndvi<- c("#2c7fb8",  
                 "#fdae61",  
                 "#a6d96a",  
                 "#1a9641",  
                 "#004529")
plot(ndvi_factor, 
     col = colores_ndvi, 
     main = "Usos según NDVI")

Igualmente, se puede usar la función terra::segregate() se puede descomponer el objeto ráster en sus diferentes categorías y representarlas gráficamente

ndvi_factor_segun_categorias <- segregate(ndvi_factor, 
                                          keep = TRUE, 
                                          other = NA)
my_colors <- c("blue", "tan", "yellow", "palegreen", "darkgreen")
par(mfrow=c(1, 5))  # set plotting layout
for (i in 1:nlyr(ndvi_factor_segun_categorias)) {
  plot(ndvi_factor_segun_categorias[[i]], col=my_colors[i], main=paste("Layer", i))
}

En la siguiente imagen se superponen todas las categorías otra vez

plot(ndvi_factor_segun_categorias[[1]], col=my_colors[1], legend=FALSE, main="Overlay of Categories")
for (i in 2:nlyr(ndvi_factor_segun_categorias)) {
  plot(ndvi_factor_segun_categorias[[i]], col=my_colors[i], legend=FALSE, add=TRUE)
}

A continuación se eliminarán del *Global Environment` la mayoría de los objetos creados anteriormente.

rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1", "imagen_veg2", "levels", "matriz_clas_suelo", "matriz_clas_veg1", "matriz_clas_veg2", "my_colors", "ndvi_clases", "ndvi_clasificado", "ndvi_factor", "ndvi_factor_segun_categorias", "umbral_suelo", "umbral_veg1", "umbral_veg2")
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'imagen_veg2' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'matriz_clas_veg2' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'umbral_veg2' no encontrado

📝 ACTIVIDAD DE EVALUACIÓN 3:

► Repite la misma actividad con el fichero LC08_016028_20210824.tif que encontrarás en el fichero comprimido, imagen multibanda correspondiente a la ciudad de Ottawa, en Canadá, tomada por el satélite Landsat 8.

¿Se podría ajustar mejor los umbrales de esta clasificación?