PREPROCESAMIENTO DE IMÁGENES DE SATÉLITE: RESOLUCIÓN DE LAS IMÁGENES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla y Francisco Conde
2023-04-01
1 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 la imagen. Está determinada por el tamaño del píxel, medido en metros sobre el terreno, esto depende de 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, mientras que los satélites meteorológicos como NOAA, el píxel representa un tamaño desde 500 a 1100m de lado.
Resolución radiométrica: Se la llama a veces también 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 usado. Así por ejemplo Landsat TM es de 28 = 256. Esto significa que tenemos una mejor resolución dinámica en el TM y podemos distinguir mejor las pequeñas diferencias de radiación.
Si en el caso de los mosaicos es condición básica disponer de imágenes con idéntica resolución, en otras situaciones los usuarios deben trabajar con datos ráster que poseen diferentes resoluciones espaciales y radiométricas. Para poder trabajor con todos ellos, debe modificarse la resolución de alguno de los datos. El remuestreo 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 (que supone disminuir el tamaño de píxel), lo que se conoce como ampliación, o bien disminuyendo la resolución espacial (lo que supone aumentar 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 remuestreo según 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.
El remuestreo según una 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.
El remuestreo según 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.
2 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.
Diferentes resoluciones espaciales de una imagen
El paquete raster contiene varias funciones para realizar tanto el remuestreo espacial como el radiométrico. Comenzaremos la actividad activando dicho paquete.
library(raster)## Loading required package: spEn primer lugar, debe importarse la imagen sobre la que se realizarán
las operaciones de remuestreo. Dado que el fichero
mosaico1es un fichero multibanda se importará como
RasterBrick.
imagen <- brick("./datos/mosaico/originales/mosaico1.tif")Sin embargo, el remuestreo espacial debe realizarse banda a banda. Por lo tanto, una vez importado el fichero raster y convertido en objeto se extraerá la banda 5 para mostrar el procedimiento.
imagen <- imagen[[5]] # Seleccionamos la banda 5
res(imagen) # Resolución espacial de la imagen## [1] 30 30ncell(imagen) # Número de celdas de la imagen importada## [1] 1699296imagen## class : RasterLayer
## band : 5 (of 7 bands)
## dimensions : 1488, 1142, 1699296 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : mosaico1.tif
## names : mosaico1.5
## values : 4546, 65535 (min, max)2.1 Aumento 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 <- raster(crs = crs(imagen), # Mismo sistema de coordenadas
resolution = 10, # Resolución del nuevo fichero
xmn = xmin(imagen), xmx = xmax(imagen), # Dimensiones en el eje x, las mismas del fichero original
ymn = ymin(imagen), ymx = ymax(imagen)) # Dimensiones en el eje y, las mismas del fichero original
imagen10## class : RasterLayer
## dimensions : 4464, 3426, 15293664 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs#values(imagen10)2.1.1 Mé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-).
Remuestreo según el vecino más próximo
El remuestreo por el método del vecino más cercano se puede realizar
usando el argumento method = "ngb" en la función
raster::resample del paquete ráster.
imagen10_ngb <- resample(imagen,
imagen10,
method = "ngb")
res(imagen10_ngb) # Resolución del nuevo fichero## [1] 10 10ncell(imagen10_ngb) # Número de celdas del nuevo fichero## [1] 15293664Puede 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.
2.1.2 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, los de un MDT o los de una imagen de satélite).
Remuestreo según el método bilineal
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] 15293664A 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 original2.2 Reducció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.
2.2.1 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 <- raster(crs = crs(imagen),
resolution = 90, # Resolución del nuevo fichero
xmn = xmin(imagen), xmx = xmax(imagen), # Dimensiones en el eje x, las mismas del fichero original
ymn = ymin(imagen), ymx = ymax(imagen)) # Dimensiones en el eje y, las mismas del fichero original- Método vecino más cercano
imagen90_ngb <- resample(imagen,
imagen90,
method = "ngb")
res(imagen90_ngb) # Resolución del nuevo fichero## [1] 90 90ncell(imagen90_ngb) # Número de celdas del nuevo fichero## [1] 188976- 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] 1889762.3 Agregació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 raster::aggregatedel paquete
ráster. 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.
2.3.1 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.
Agregación mediante la función media
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] 1889762.3.2 Con 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.
Agregación mediante la función moda
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] 1889762.3.3 Con 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.
Agregación mediante la función mínimo
En R, el argumento fun = minse usa en la función
aggregate()
imagen90_minimo <- aggregate(imagen,
fact = 3,
fun = min,
na.rm = TRUE)
res(imagen90_minimo) ## [1] 90 90ncell(imagen90_minimo)## [1] 1889762.3.4 Con 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.
Agregación mediante la función maximo
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] 188976A 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
raster:::disaggregate(), pero a diferencia del primero no
requiere una función, ya que solo asigna sus valores a píxeles más
pequeños. Si revisa los resultados, verá 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 <- disaggregate(imagen,
fact = 2)
res(imagen_desagregada)## [1] 15 153 REMUESTREO 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.
Remuestreo según el vecino más próximo
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.
- La sintaxis de la función
RStoolbox::rescalees la siguiente:
rescaleImage(x, y, xmin, xmax, ymin, ymax, forceMinMax = FALSE)
Siendo los argumentos:
x Objeto ráster a normalizar.
y Objeto ráster que se comporta como la imagen
de referencia (aunque este argumento es opcional). Se usa para extraer
los valores min y max si se omiten.
xmin Valor minimo de x. O un solo valor o un
valor por capa en x. Si no se proporciona xmin, se extraerá de x.
xmax Valor máximo de x. O un solo valor o un
valor por capa en x. Si no se proporciona xmin, se extraerá de x.
ymin Valor mínimo de y. O un solo valor o un
valor por capa en y. Si no se proporciona xmin, se extraerá de y.
ymax Valor máximo de y. O un solo valor o un
valor por capa en y. Si no se proporciona xmin, se extraerá de y.
forceMinMax Valor lógico. Obliga a una
actualización de los valores mínimo y máximo en x o en y.
En el siguiente ejemplo, la imagen original, cuyo tipo es INT2U, con un tamaño de 16 bits y valores radiométricos de 0 a 65535, se transformará a una imagen de 8 bits con valores de píxeles de 0 a 255 y tipo de datos INT1U.
En primer lugar, se importará la imagen.
imagen <- raster("./datos/remuestreo/originales/LC08_L1TP_201032_20220507_20220512_02_T1_B5.tif")
imagen## class : RasterLayer
## dimensions : 7911, 7791, 61634601 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 314385, 548115, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : LC08_L1TP_201032_20220507_20220512_02_T1_B5.TIF
## names : LC08_L1TP_201032_20220507_20220512_02_T1_B5
## values : 0, 65535 (min, max)A continuación, con la función raster::dataType
preguntamos a R que tipo de datos contiene la imagen
dataType(imagen) ## [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)A continuación, la función rescaleImage adapta los
valores del primer formato al segundo.
img8bit <- rescaleImage(imagen,
xmin = 0,
xmax = 65535,
ymin= 0,
ymax= 255)
dataType(img8bit)## [1] "FLT4S"img8bit## class : RasterLayer
## dimensions : 7911, 7791, 61634601 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 314385, 548115, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 16.92607, 255 (min, max)- 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 <- calc(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_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/remuestreo/resultados/img8bit.tif",
datatype = "INT1U",
overwrite = TRUE)A continuación se importará de nuevo el ráster para verificar sus datos.
img8bit_tris <- raster("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/remuestreo/resultados/img8bit.tif") # Abre de nuevo el ráster
dataType(img8bit_tris) # Verifica el tipo de datos del nuevo fichero raster## [1] "INT1U"3.0.1 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
raster:cellStats, que devuelve un único valor calculado a
partir de todos los valores de una capa.
a <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(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 <- cellStats(imagen90_maximo, stat = "mean", na.rm = TRUE)- Promedio del raster remuestreado radiométricamente a una imagen de 8 bit.
j <- cellStats(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## [1] 18190.96078 18365.23541 18365.25908 18365.69953 18366.10569 18365.02374
## [7] 18363.86504 16562.59332 20286.85298 70.78195Crea un dataframe con esos valores medios
nombre_procedimiento <- c("Raster original", "Incremento resolución según el vecino más próximo", "Incremento resolución según interpolación bilineal", "Reducción resolución según el vecino más próximo", "Reducción resolución según la interpolación bilineal", "Agregación según la media", "Agregación según la moda","Agregación según el mínimo", "Agregación según el máximo", "Remuestreo radiométrico a 8 bit")
data.frame(promedio_resample, row.names = nombre_procedimiento)## promedio_resample
## Raster original 18190.96078
## Incremento resolución según el vecino más próximo 18365.23541
## Incremento resolución según interpolación bilineal 18365.25908
## Reducción resolución según el vecino más próximo 18365.69953
## Reducción resolución según la interpolación bilineal 18366.10569
## Agregación según la media 18365.02374
## Agregación según la moda 18363.86504
## Agregación según el mínimo 16562.59332
## Agregación según el máximo 20286.85298
## Remuestreo radiométrico a 8 bit 70.78195A partir de este dataframe, 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())