Como bien hemos visto en anteriores artículos, R es un software con una enorme capacidad para el tratamiento de información espacial y es por eso que en los últimos años ha conseguido llegar a un elevado número de usuarios en este ámbito. Sin embargo, esta popularidad creciente va de la mano de la liberalización de la información espacial, favorecida por la puesta en marcha de misiones de gran calibre como, por ejemplo, Sentinel 2.
La misión de la Agencia Espacial Europea (ESA) ha permitido acceder de forma libre a imágenes multiespectrales de alta resolución espacial para toda la Tierra. Además, al estar compuesto por dos satélites hermanos (Sentinel-2A y Sentinel-2B), permite obtener imágenes con una alta resolución temporal, con imágenes disponibles cada 5 días de media.
Debido a esta mayor accesibilidad a los datos y a la mayor calidad de los mismos, han aparecido nuevas posibilidades y necesidades de explotación. Uno de estos ejemplos son las firmas espectrales.
¿Qué son las firmas espectrales?
Las firmas espectrales o también llamadas signaturas espectrales son una caracterización de la reflectividad de una superficie u objeto para unas longitudes de onda determinadas. Esta característica, única para cada tipo de superficie, permite diferenciar e identificar distintos elementos de la superficie terrestre, hasta el punto de poder clasificarlos por proximidad espectral.
El comportamiento, depende de las características espectrales de cada superficie y, en ciertas superficies, de la variabilidad temporal de esta. Estas diferencias son especialmente perceptibles en superficies vegetales, en las que existe esa alta variabilidad debido a las moderadas diferencias de reflectividad según el estado fotosintético de la planta.
Dependiendo de los valores de reflectividad obtenidos a través del sensor, podemos ser incluso capaces de diferenciar el estado vegetativo de la planta y discriminar entre vegetación sana o enferma.
Acceder a la información de Sentinel 2
Como ya hemos visto en anteriores artículos, podemos acceder a información ráster a través de R con la librería raster, ya instalada en el núcleo de R. Para ello vamos a cargar la librería y vamos a listar todos los archivos de nuestra carpeta con la función list.files(). Vamos a utilizar el parámetro full.names para devolver la ruta completa a cada uno de los archivos.
#Cargar librería raster
library(raster)
path_folder = "E:/RUTA_CARPETA/S2B_MSIL1C_20190514_bandas/"
jp2_files <- list.files(path_folder, full.names = T)A continuación, vamos a filtrar los archivos ráster seleccionados de nuestra carpeta por nombre o patrón de texto.
pattern <- c("B02", "B03", "B04", "B05", "B06", "B07","B08", "B11", "B12")
stack_files <- unique(grep(paste(pattern,collapse="|"), jp2_files, value=TRUE))Seleccionados todos los archivos en los que aparecen los nombres indicados, vamos a crear un objeto RasterStack con el conjunto de rásters seleccionados.
sen2 <- stack(stack_files)Seleccionar coordenadas utilizando la función click() del paquete raster
En este punto, vamos a crear un data.frame vacío, en el que vamos a introducir distintas coordenadas haciendo clic sobre el mapa, con la ayuda de la función click() del paquete raster. Cada uno de estos clics (hemos definido 10 en este caso), identificará las coordenadas en el mapa y serán introducidas en nuestra tabla de datos.
df_clic <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(df_clic) <- c("x","y")
plotRGB(brick(jp2_files[14]), r=1, g=2, b=3)
for (c in 1:10){
df_clic <- rbind(df_clic,as.data.frame(click()))
}Ahora vamos a convertir la tabla con coordenadas (x,y) a formato SpatialPointsDataFrame,con el que podremos representar nuestros puntos de forma gráfica.
spdf <- SpatialPointsDataFrame(coords = df_clic[,c(1,2)], data = df_clic, proj4string = crs(sen2))
plotRGB(brick(jp2_files[14]), r=1, g=2, b=3)
text(spdf,1:10,cex=1.5)
Extracción de valores ráster por coordenadas
Tras obtener los puntos, vamos a realizar la extracción de los valores ráster para cada una de las localizaciones seleccionadas. Para ello vamos a utilizar la función extract() del paquete raster y vamos a convertir el resultado a un nuevo data.frame, al que le vamos a añadir las longitudes de onda para cada una de las bandas. Además, vamos a modificar los nombres de cabecera de las columnas y filas con colnames() y rownames().
ext_sen <- as.data.frame(t(extract(sen2, spdf)))
print(ext_sen)
ext_sen <- cbind(ext_sen,c(492,560,665,705,740,783,842,1610,2190))
rownames(ext_sen) <- c("blue","green","red","RE1","RE2","RE3","NIR","SWIR1","SWIR2")
colnames(ext_sen) <- c("1","2","3","4","5","6","7","8","9","10","bands")Modificar data.frame a formato largo con melt (librería reshape)
Vamos a convertir nuestra tabla de datos a formato largo utilizando al función melt() de la librería reshape. Este formato se utiliza para que cada fila contenga un dato y a través del resto de columnas se contextualiza el mismo. En este caso, queremos obtener un valor de reflectividad para cada coordenada y banda. También vamos a escalar nuestros valores al intervalo 0-100.
library(reshape)
ext_melt <- melt(ext_sen,id=c("bands"))
ext_melt$value <- ext_melt$value/10000*100Representar firmas espectrales con ggplot2
Ya tenemos todos nuestros datos listos para ser representados en un nuevo gráfico con ggplot2. Antes, vamos a crear una paleta de 10 colores aleatorios con la función distinctColorPalette() de la librería randomcoloR.
library(randomcoloR)
palette <- distinctColorPalette(10)
library(ggplot2)
ggplot(data = ext_melt, aes(x=bands, y=value, color=variable)) +
geom_line(size=1.3, alpha=0.6) + geom_point(shape=19, size=2) + scale_color_manual(values = palette) + scale_x_sqrt(breaks=c(492,560,665,705,740,783,842,1610,2190)) + ylim(0, 100) + annotate("text", x = c(492,560,665,705,740,783,842,1610,2190), y=0, label = c("blue","green","red","RE1","RE2","RE3","NIR","SWIR1","SWIR2"), size=2.5) + theme_bw(base_size = 10) + labs(title="Firmas espectrales Sentinel-2", x="Longitud de onda (μm)", y="Reflectividad (%)")
¿Para qué podremos utilizar firmas espectrales en teledetección?
- Las firmas espectrales nos permiten identificar coberturas u objetos de la superficie terrestre con características espectrales similares.
- Pueden ser utilizadas para calibrar imágenes satélites sin corrección radiométrica.
- Permiten distinguir estados fenológicos de los cultivos o vegetación en general.
Aquí el código completo:
library(raster)
path_folder = "E:/ARTICULOS/R/4.FIRMAS_ESPECTRALES/S2B_MSIL1C_20190514_bandas/"
jp2_files <- list.files(path_folder, full.names = T)
pattern <- c("B02", "B03", "B04", "B05", "B06", "B07","B08", "B11", "B12")
stack_files <- unique(grep(paste(pattern,collapse="|"), jp2_files, value=TRUE))
sen2 <- stack(stack_files)
df_clic <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(df_clic) <- c("x","y")
plotRGB(brick(jp2_files[14]), r=1, g=2, b=3)
for (c in 1:10){
df_clic <- rbind(df_clic,as.data.frame(click()))
}
spdf <- SpatialPointsDataFrame(coords = df_clic[,c(1,2)], data = df_clic, proj4string = crs(sen2))
text(spdf,1:10,cex=1.5)
ext_sen <- as.data.frame(t(extract(sen2, spdf)))
print(ext_sen)
ext_sen <- cbind(ext_sen,c(492,560,665,705,740,783,842,1610,2190))
rownames(ext_sen) <- c("blue","green","red","RE1","RE2","RE3","NIR","SWIR1","SWIR2")
colnames(ext_sen) <- c("1","2","3","4","5","6","7","8","9","10","bands")
library(reshape)
ext_melt <- melt(ext_sen,id=c("bands"))
ext_melt$value <- ext_melt$value/10000*100
library(randomcoloR)
palette <- distinctColorPalette(10)
library(ggplot2)
ggplot(data = ext_melt, aes(x=bands, y=value, color=variable)) +
geom_line(size=1.3, alpha=0.6) +
geom_point(shape=19, size=2) +
scale_color_manual(values = palette) +
scale_x_sqrt(breaks=c(492,560,665,705,740,783,842,1610,2190)) +
ylim(0, 100) +
annotate("text", x = c(492,560,665,705,740,783,842,1610,2190), y=0,
label = c("blue","green","red","RE1","RE2","RE3","NIR","SWIR1","SWIR2"), size=2.5) +
theme_bw(base_size = 10) +
labs(title="Firmas espectrales Sentinel-2", x="Longitud de onda (μm)", y="Reflectividad (%)")
¿Quieres comentarnos algo? Adelante!