LABORATORIO 3. PREPROCESAMIENTO DE IMÁGENES DE SATÉLITE: RESOLUCIÓN DE LAS IMÁGENES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla
2025-03-10
💡 OBJETIVOS DE LA ACTIVIDAD:
► Diferenciar los conceptos de resolución espacial y radiométrica.
► Familiarizarse con los procedimientos necesarios para modificar la resolución de las imágenes de satélite, conociendo sus ventajas y sus inconvenientes.
💡 MATERIALES PARA LA ACTIVIDAD:
Los materiales para el desarrollo de esta actividad son los mismos de la actividad pasada.
► Imágenes originales. Contiene los siguientres ficheros:
INTRODUCCIÓN
Una imagen de satélite se caracteriza por diferentes tipos de resolución:
Resolución espacial.
Resolución espectral.
Resolución radiométrica.
Resolución temporal.
En este apartado se trabajará fundamentalmente con la resolución espacial y la resolución radiométrica.
Resolución espacial: este concepto designa al objeto más pequeño que se puede distinguir en una imagen. Está determinada por el tamaño del píxel (normalmente medido en unidades de longitud -metros- sobre el terreno), la altura del sensor con respecto a la Tierra, el ángulo de visión, la velocidad de escaneado y las características ópticas del sensor.
Por ejemplo las imágenes Landsat TM, tienen una resolución espacial de 30x30 m en las bandas 1,2,3,4, 5 y 7 y de 120x120m en la 6 (térmica).
El sensor SPOT - HRV tiene una resolución de 10x10m.
Satélites meteorológicos como METEOSAT tienen un píxel con una tamaño de 5 km.
Resolución radiométrica: se la conoce también por resolución dinámica, y se refiere a la cantidad de niveles de gris en que se divide la radiación recibida para ser almacenada y procesada posteriormente. Esto depende del conversor analógico digital incluido en el sensor. Por ejemplo, en Landsat TM es de 256. Una mayor resolución facilita el distinguir mejor las pequeñas diferencias de radiación.
El remuestreo (resampling) es el proceso para modificar la resolución de los píxeles de una imagen, siendo uno de las fases más comunes en cualquier análisis que implique el uso de imágenes de satélite.
El remuestreo espacial se realiza bien aumentando la resolución espacial (reducir el tamaño de píxel), lo que se conoce como ampliación, o bien disminuyendo la resolución espacial (aumentando el tamaño de píxel), conocido como reducción. El proceso de remuestreo calcula los nuevos valores de los píxeles a partir de los valores de la imagen no corregida, siendo tres los métodos más comunes:
El vecino más cercano (nearest neighbour) utiliza el valor del píxel de la imagen original que está más cerca de la nueva ubicación del píxel en la imagen corregida. Este es el método más simple y no altera los valores originales, pero puede provocar que algunos valores de píxeles se dupliquen mientras que otros se pierden, y tiende a dar como resultado una apariencia de imagen inconexa o en bloques.
La interpolación bilineal pondera los cuatro píxeles en la imagen original más cercanos a la nueva ubicación del píxel de salida, creando valores digitales completamente nuevos. Esto puede ser indeseable si se va a realizar un procesamiento y análisis adicionales, como la clasificación basada en la respuesta espectral.
La convolución cúbica va aún más lejos calculando un promedio ponderado según la distancia de un bloque de dieciséis píxeles de la imagen original que rodea la nueva ubicación del pixel de salida.
El remuestreo radiométrico se realiza para aumentar o disminuir la resolución radiométrica de la imagen.
REMUESTREO ESPACIAL
Uno de los inconvenientes de los objetos ráster es que suelen ofrecer resoluciones espaciales diferentes. En el caso de las imágenes de satélite, algunas bandas disponen de mayor resolución que otras. Esta situación supone una reducción en la nitidez de los elementos que aparecen en la imagen, lo que dificulta su interpretación.
El paquete terra contiene varias funciones para realizar tanto el remuestreo espacial como el radiométrico. Comenzaremos la actividad activando dicho paquete.
library(terra)## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.21setwd("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/resolucion_imagenes/")En primer lugar, debe importarse la imagen sobre la que se realizarán las operaciones de remuestreo.
imagen <- rast("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/resolucion_imagenes/fichero_remuestreo.tif")A continuación, se puede preguntar a R acerca de la resolución de la imagen y del número de celdas.
res(imagen) # Resolución del fichero## [1] 30 30ncell(imagen) ## [1] 59763461Aumento de la resolución espacial
Este remuestreo aumenta la resolución espacial, es decir, disminuye el tamaño de píxel del objeto ráster.
En primer lugar, necesitamos crear un ráster vacío con la resolución espacial del nuevo ráster (10 m en este caso). Este nuevo ráster actuará como ráster de salida (o final), y cuyos valores de píxel se calcularán a partir del ráster de entrada, tras aplicar alguno de los métodos de remuestreo.
IMPORTANTE. Este nuevo ráster debe poseer el mismo sistema de coordenadas de referencia del raster de entrada, así como sus mismas dimensiones.
Para crear un nuevo ráster que cumpla ambas condiciones la sintaxis será:
imagen10 <- rast(crs = crs(imagen), # Mismo sistema de coordenadas
resolution = 10, # Resolución del nuevo fichero
xmin = xmin(imagen), xmax = xmax(imagen), # Dimensiones en el eje x, las mismas del fichero original
ymin = ymin(imagen), ymax = ymax(imagen)) # Dimensiones en el eje y, las mismas del fichero original
imagen10## class : SpatRaster
## dimensions : 23403, 22983, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 308985, 538815, 2757585, 2991615 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 44N (EPSG:32644)res(imagen10) # Resolución del nuevo fichero## [1] 10 10ncell(imagen10) ## [1] 537871149Método del vecino mas cercano.
En el método del vecino más cercano, el valor de píxel de la imagen ráster de salida se calcula a partir del píxel más cercano de la imagen de entrada. Es un método muy simple que no cambia el valor de la celda de entrada en el ráster de salida. Por ello, es el más recomendable para su aplicación sobre datos discretos (por ejemplo, un raster con usos del suelo -pe Corine Land cover-).
El remuestreo por el método del vecino más cercano se puede realizar
usando el argumento method = "near" en la función
terra::resample.
imagen10_ngb <- resample(imagen,
imagen10,
method = "near")## |---------|---------|---------|---------|========================================= res(imagen10_ngb) # Resolución del nuevo fichero## [1] 10 10ncell(imagen10_ngb) # Número de celdas del nuevo fichero## [1] 537871149Puede observarse que el número de celdas es mucho mayor que el de la imagen original, ya que la resolución espacial es ahora de 10 m.
Método “bilineal”.
En el método de interpolación bilineal, el valor de los píxeles en el raster de salida se calcula como el promedio ponderado de la distancia de los cuatro píxeles más cercanos, lo que proporciona una apariencia más suave a ese ráster. Es más apropiado para el remuestreo con datos continuos (por ejemplo, un MDT o una imagen de satélite).
El remuestreo por este método se puede realizar usando el argumento
method = "bilinear".
imagen10_bl <- resample(imagen,
imagen10,
method = "bilinear")## |---------|---------|---------|---------|========================================= res(imagen10_bl) # Resolución del nuevo fichero## [1] 10 10ncell(imagen10_bl) # Número de celdas del nuevo fichero## [1] 537871149A continuación se representan ambas imágenes para observar el resultado de las mismas
# Una fila, dos columnas
par(mfrow = c(1, 2))
# Gráficos a combinar
plot(imagen10_ngb,
main = "Vecino más próximo",
col = gray(0:100/100))
plot(imagen10_bl,
main = "Método bilineal",
col = gray(0:100/100))
par(mfrow = c(1, 1)) # Volvemos al estado originalReducción de la resolución espacial
La reducción de resolución espacial aumenta el tamaño de píxel del conjunto de datos ráster. Los métodos de reducción son similares a los empleados en el caso de aumento de la resolución espacial.
Métodos del vecino más próximo y bilineal
El método del vecino más cercano produce un efecto granulado y grueso en el ráster de salida. Para aplicar el vecino más cercano y el método de interpolación bilineal para el remuestreo en R, primero tenemos que crear un ráster vacío de menor resolución espacial, digamos 90 m en este caso.
imagen90 <- rast(crs = crs(imagen),
resolution = 90, # Resolución del nuevo fichero
xmin = xmin(imagen), xmax = xmax(imagen), # Dimensiones en el eje x, las mismas del fichero original
ymin = ymin(imagen), ymax = ymax(imagen)) # Dimensiones en el eje y, las mismas del fichero original► Método vecino más cercano
imagen90_ngb <- resample(imagen,
imagen90,
method = "near")
res(imagen90_ngb) # Resolución del nuevo fichero## [1] 90 90ncell(imagen90_ngb) # Número de celdas del nuevo fichero## [1] 6640400► Método bilinear
imagen90_bl <- resample(imagen, imagen90, method = "bilinear")
res(imagen90_bl) # Resolución del nuevo fichero## [1] 90 90ncell(imagen90_bl) # Número de celdas del nuevo fichero## [1] 6640400Agregación espacial
A veces es necesario cambiar la resolución espacial de los datos ráster para investigar los efectos de escala, puesto que crear una secuencia de resoluciones espaciales e investigar cual es su influencia en los resultados de una investigación puede ayudar a identificar qué resolución espacial es la más adecuada para cartografiar la distribución espacial de un determinado proceso.
También puede reducir la escala de un ráster subdividiendo los
píxeles, si bien esto no tiene ninguna influencia en los valores de
píxeles; simplemente aumenta el número de píxeles. La disminución de la
resolución espacial de un ráster se logra agregando los valores de los
píxeles vecinos y asignándolos a un nuevo píxel con una menor resolución
espacial. Para valores continuos, esto generalmente significa promediar
los valores de píxeles. Para hacer esto, usamos
raster::aggregate(), modificando la escena para que la
resolución original de 30 m pase a 60 m (factor 2).
Sin embargo, la reducción de escala puede ser útil para hacer que las resoluciones de dos imágenes sean compatibles.
La agregación espacial es un procedimiento similar al del remuestreo, pero en este caso se mantiene la estructura de píxeles del ráster de entrada, calculándose los nuevos valores del raster de salida a partir de los resultados de una función específica proporcionada por el usuario (por ejemplo la media, modo, mínimo, máximo, etc).
En R, la agregación espacial de las imágenes ráster se realiza
mediante la función terra::aggregate[https://rspatial.github.io/terra/reference/aggregate.html].
Las funciones utilizadas para agregar el valor al ráster de salida se
especifican en el argumento divertido fun. El número de
celdas para agregar tanto en dirección horizontal como vertical se
menciona en el factor de agregación con el argumento fact,
que es en última instancia quién determina la resolución espacial del
ráster de salida.
En las secciones posteriores se muestran diferentes vías de agregación espacial para obtener un ráster de salida con una resolución espacial de 90 m desde un ráster de entrada de resolución de 30 m mediante las funciones media, modo, mínimo y máximo, si bien esta función puede ser modificada en función del interés del usuario.
Por último, téngase en cuenta que si tenemos valores discretos en nuestros píxeles (como en una clasificación), debemos evitar promediar porque el promedio numérico de la clase “1” y la clase “2” daría como resultado una clase indefinida “1.5”. Por lo tanto, para los datos categóricos normalmente usaríamos el valor más frecuente (fun = modal), que, por ejemplo, devolvería “bosque” para el conjunto de valores c(“Bosque”, “Bosque”, “Agua”, “Urbano”). Este procedimiento se conoce como filtro mayoritario.
Con el valor promedio.
Esta función calcula la media de los valores de los píxeles situados dentro de la ventana especificada por el usuario.
Para realizar este se asigna el argumento fun = mean;
además, se asigna el argumento fact = 3” para crear un
ráster con una resolución de píxeles de 90m, al multiplicar la
resolución espacial de 30 m del ráster de entrada por 3 (3x30=90). Si
desea eliminar el NA o los datos faltantes del cálculo para la
agregación espacial, se utiliza el argumento na.rm = TRUE.
imagen90_media <- aggregate(imagen,
fact = 3,
fun = "mean",
na.rm = TRUE)## |---------|---------|---------|---------|========================================= res(imagen90_media) ## [1] 90 90ncell(imagen90_media)## [1] 6642954Con el valor modal.
Los valores de píxel con mayor frecuencia, la moda de los valores de píxeles en el ráster de entrada se asigna al ráster de salida como se muestra en la siguiente figura.
En este caso se usa el argumento fun = modal.
imagen90_moda <- aggregate(imagen,
fact = 3,
fun = "modal",
na.rm = TRUE)## |---------|---------|---------|---------|========================================= res(imagen90_moda) ## [1] 90 90ncell(imagen90_moda)## [1] 6642954Con el valor mínimo.
Esta función usa el valor mínimo de los píxeles del ráster de entrada correspondientes a la ventana espacial especificada por el usuario.
imagen90_minimo <- aggregate(imagen,
fact = 3,
fun = "min",
na.rm = TRUE)## |---------|---------|---------|---------|========================================= res(imagen90_minimo) ## [1] 90 90ncell(imagen90_minimo)## [1] 6642954Con el valor máximo.
El valor máximo de los píxeles del raster de entrada en la ventana especificada por el usuario es asignado como valor de los píxeles en el ráster de salida.
En R, el argumento fun = maxse usa en la función
aggregate()
imagen90_maximo <- aggregate(imagen,
fact = 3,
fun = "max",
na.rm = TRUE)## |---------|---------|---------|---------|========================================= res(imagen90_maximo) ## [1] 90 90ncell(imagen90_maximo)## [1] 6642954A continuación se representan ambas imágenes para observar el resultado de las mismas
# Una fila, dos columnas
par(mfrow = c(2, 2))
# Gráficos a combinar
plot(imagen90_media,
main = "Valor promedio",
col = gray(0:100/100))
plot(imagen90_moda,
main = "Valor modal",
col = gray(0:100/100))
plot(imagen90_minimo,
main = "Valor mínimo",
col = gray(0:100/100))
plot(imagen90_maximo,
main = "Valor máximo",
col = gray(0:100/100))
par(mfrow = c(1, 1)) # Volvemos al estado originalEl complemento de aggregate() es
terra:::disagg()[https://rspatial.github.io/terra/reference/disaggregate.html],
pero a diferencia del primero no requiere una función, ya que solo
asigna sus valores a píxeles más pequeños. Revisando los resultados se
observa que la resolución ha vuelto a 30 m, pero la imagen es
exactamente la misma de la versión a 60 m.
imagen_desagregada <- disagg(imagen,
fact = 2) ## |---------|---------|---------|---------|============== res(imagen_desagregada)## [1] 15 15REMUESTREO RADIOMÉTRICO
El concepto “resolución radiométrica” hace referencia a la cantidad de niveles digitales en los que es posible guardar la información digital contenida en una imagen. Cada imagen se guarda en bytes y cada byte está compuesto por 8 bits. Un valor, computacionalmente hablando, se puede guardar como un «0» o un «1», decir que algo tiene 8 bits de información hace referencia a la posibilidad de combinar esos ceros y unos en 8 sub espacios diferentes, es decir, 2^8 que es igual a 256 posibilidades.
Pues bien la resolución radiométrica hace referencia a cuanta información distingue el sensor: a mayor resolución radiométrica, mayor posibilidad de que la información sea más variada.
El remuestreo radiométrico aumenta o disminuye la resolución
radiométrica de la imagen ráster. Para ello, es posible utilizar dos
procedimientos: la función RStoobox::rescale() del paquete
RStoolbox, o una función creada por el propio
usuario.
En primer lugar, se importará la imagen.
imagen2 <- rast("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/resolucion_imagenes/LC08_L1TP_201032_20220507_20220512_02_T1_B5.tif")
imagen2## class : SpatRaster
## dimensions : 7911, 7791, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 314385, 548115, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source : LC08_L1TP_201032_20220507_20220512_02_T1_B5.tif
## name : LC08_L1TP_201032_20220507_20220512_02_T1_B5A continuación, con la función terra::datatype[https://rdrr.io/cran/terra/man/datatype.html]
preguntamos a R que tipo de datos contiene la imagen
datatype(imagen2) ## [1] "INT2U"El remuestreo puede realizarse de dos maneras:
► OPCIÓN 1: a través de la función
RStoolbox::rescaleImage, para lo que habrá que activar el
paquete Rstoolbox.
library(RStoolbox)## Warning: package 'RStoolbox' was built under R version 4.3.3## This is version 1.0.0 of RStoolboxA continuación, la función rescaleImage adapta los
valores del primer formato al segundo.
img8bit <- rescaleImage(imagen2,
xmin = 0,
xmax = 65535,
ymin= 0,
ymax= 255)
datatype(img8bit)## [1] ""img8bit## class : SpatRaster
## dimensions : 7911, 7791, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 314385, 548115, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source(s) : memory
## name : LC08_L1TP_201032_20220507_20220512_02_T1_B5
## min value : 16.92607
## max value : 255.00000► OPCIÓN 2. A través de una fórmula ad-hoc creada por el usuario para convertir los valores de los píxeles de 8 bit a 255 bit
img8bit_bis <- app(imagen,
fun = function(x){((x-min(x))*255)/
(max(x)-min(x)) + 0})A continuación se representan ambas imágenes para observar el resultado de las mismas
# Una fila, dos columnas
par(mfrow = c(1, 2))
# Gráficos a combinar
plot(imagen,
main = "Imagen original",
col = gray(0:100/100))
plot(img8bit,
main = "Imagen 8 bit",
col = gray(0:100/100))
par(mfrow = c(1, 1)) # Volvemos al estado originalCualquiera de estas dos imágenes puede salvarse en el disco duro como fichero GeoTiff, eso sí, especificando ahora
writeRaster(img8bit,
"D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/resolucion_imagenes/img8bit.tif",
datatype = "INT1U",
overwrite = TRUE)A continuación se importará de nuevo el ráster para verificar sus datos.
img8bit_tris <- rast("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/resolucion_imagenes/img8bit.tif") # Abre de nuevo el ráster
datatype(img8bit_tris) # Verifica el tipo de datos del nuevo fichero raster## [1] "INT1U"Comparación de la media de imágenes ráster remuestreadas con los métodos discutidos anteriormente
Con el fin de comparar los datos ráster de entrada con los resultados
obtenidos del remuestreo, se procede a continuación al cáluclo del
promedio de los valores de cada imagen ráster remuestreada por los
métodos discutidos anteriormente. Para ello se utiliza la función
terra:global, que devuelve un único valor calculado a
partir de todos los valores de una capa.
a <- global(imagen,
stat = "mean",
na.rm = TRUE) # Ignora los valores NODATA en los cálculos► Promedio del raster remuestreado a una resolución de 10m por el método del vecino más próximo
b <- global(imagen10_ngb, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 10m por el método de interpolación bilineal.
c <- global(imagen10_bl, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método del vecino más próximo.
d <- global(imagen90_ngb, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método de interpolación bilineal.
e <- global(imagen90_bl, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método de agragación por la media.
f <- global(imagen90_media, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método de agragación por la moda.
g <- global(imagen90_moda, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método de agragación por la mínima.
h <- global(imagen90_minimo, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado a una resolución de 90m por el método de agragación por el máximo.
i <- global(imagen90_maximo, stat = "mean", na.rm = TRUE)► Promedio del raster remuestreado radiométricamente a una imagen de 8 bit.
j <- global(img8bit, stat = "mean", na.rm = TRUE)► Comparación de la media de todas las imágenes
promedio_resample <- c(a, b, c, d, e, f, g, h, i, j)
promedio_resample## $mean
## [1] 11556.87
##
## $mean
## [1] 11556.87
##
## $mean
## [1] 11556.87
##
## $mean
## [1] 11556.91
##
## $mean
## [1] 11556.85
##
## $mean
## [1] 11552.4
##
## $mean
## [1] 11056.35
##
## $mean
## [1] 11046.71
##
## $mean
## [1] 12058.88
##
## $mean
## [1] 70.78195A partir de estos datos, se deduce que la media de los valores de píxel en los remuestreos es bastante similar, mientras que en el caso de los procedimientos de agregación, el valor medio de las imágenes varía según el modelo utilizado.
rm(list=ls())