Introducción a R en el tratamiento de imágenes satélite
Tras el éxito del webinar (link al video al final del post) que realizamos con nuestros amigos de Geoinnova, y el aluvión de emails que hemos recibido para preguntarnos por el script, os dejamos esta entrada del blog en el que se recopilan, uno a uno, los pasos realizados para iniciarte en R y usarlo como introducción a R en el tratamiento de imágenes de teledetección.
¿Qué vamos a aprender?
- Instalación de paquetes necesarios para el tratamiento de imágenes en R.
- Descarga de imágenes de satélite
- Cálculo de estadísticas de imágenes
- Cálculo de índices aplicados a las imágenes
También podrás encontrar más abajo los comandos que se utilizarán durante el webinar. De este modo solamente tendrás que copiar y pegar en RStudio y podrás replicar todo lo visto de manera más cómoda, a la vez que transcurre el webinar o posteriormente.
Antes de empezar
Tendremos que hacer dos cosas, la primera será instalar RStudio y la segunda registrarnos en la plataforma para descargas de imágenes de satélite Sentinel de la Agencia Espacial Europea (ESA).
Instalación de RStudio
Utilizaremos RStudio durante el desarrollo del webinar.
RStudio es un IDE de código abierto desarrollado para facilitar el trabajo con el lenguaje R que nos permite hacer scripts para análisis de datos, aplicaciones web interactivas, reportes, gráficos y, lo que más nos interesa a nosotros, añadir paquetes que nos permitirán hacer teledetección.
Existe una versión libre y gratuita, con licencia AGPL v3, que puedes descargar desde aquí: https://rstudio.com/products/rstudio/download/
Registro en Copernicus Open Access Hub
El primer paso que tendremos que realizar será acceder a Copernicus Open Access Hub y registrarnos en https://scihub.copernicus.eu/ para poder descargar imágenes Sentinel.
A partir de aquí, encontrarás los comandos para que copies y pegues en RStudio y generes tu propio script en lenguaje R.
Instalar paquete sen2r
install.packages("sen2r")Cargar librería sen2r
library(sen2r)Login
write_scihub_login(username = getPass::getPass(), password = getPass::getPass())Definir ventana temporal
time_window <- as.Date(c("2020-05-01", "2020-05-05"))Tile – Google Earth Pro
#https://sentinel.esa.int/documents/247904/1955685/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml
tile = "30TXM"Limitar el % de nubes de la imagen
max_cloud = 20Comprobar productos disponibles
list_safe <- s2_list(tile= tile, time_interval = time_window, max_cloud = max_cloud)Metadatos
safe_getMetadata(list_safe, "sensing_datetime")
safe_getMetadata(list_safe, "nameinfo")Crear directorio de descarga
path <- "C:/WEBINAR_R/" #Indicar la ruta a la carpeta donde está guardada la imagen
path_save <- paste0(path,"S2A")
dir.create(path_save)Descargar
s2_download(list_safe, outdir = path_save)Cargar imágenes raster en R
En teledetección es muy importante conocer las características de las imágenes con las que trabajamos. Aquí tenéis un par de links que dan acceso a la documentación de Sentinel-2.
- https://sentinel.esa.int/web/sentinel/missions/sentinel-2
- https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/spatial
Comprobar directorios
dirs <- list.dirs(path = path_save)
dirs <- dirs[grep("IMG_DATA", dirs)]
v_files <- list.files(dirs[1], full.names = T, recursive = T)Acceder a un objeto
print(v_files[2])Filtrar por patrón de texto
v_files <- list.files(path = path_save, pattern = "20m.jp2", full.names = TRUE, recursive = T)Instalar y cargar la librería raster de R
install.packages("raster") #Instalación de la librería
library(raster) #Cargamos la libreríaColorear imagen en color verdadero
plotRGB(brick(v_files[12]), r=1, g=2, b=3)Cargar banda 2 (azul) desde archivo – ruta
b02 <- raster(v_files[2])Explorar información de la imagen
crs(b02) #CRS
res(b02) #Resolución X - y
xres(b02) #Resolución X
yres(b02) #Resolución y
ncell(b2) #Número de celdas
dim(b02) #Dimensiones
extent(b02) #Extensión x-yGuardar valores del raster como un objeto
b02_hist <- getValues(b02)Crear el histograma
hist(b02_hist, breaks=5000, xlim=c(0,3000)) #HistogramaRecortar la imagen (“Clip”)
xy_ext <- extent(xmin(b02),xmax(b02),ymin(b02),ymax(b02))
xy_ext <- extent(255000, 260000, -375000, -370000)Clip cuadrado central de la imagen (1/4)
xmax <- xmax(b02) - ((abs(xmax(b02)-xmin(b02))/2)/2)
xmin <- xmin(b02) + ((abs(xmax(b02)-xmin(b02))/2)/2)
ymax <- ymax(b02) - ((abs(ymax(b02)-ymin(b02))/2)/2)
ymin <- ymin(b02) + ((abs(ymax(b02)-ymin(b02))/2)/2)Cuadro delimitador
xy_ext <- extent(xmin,xmax,ymin,ymax)
plot(b02)
plot(xy_ext,add=T)Crop
b02_crop <- crop(b02,xy_ext)
plot(b02_crop)Filtrar archivos por nombre
patterns <- c("B02", "B03", "B04", "B05", "B06", "B07", "B11", "B12", "B8A")
brick_files <- unique(grep(paste(patterns,collapse="|"), v_files, value=TRUE))Crear stack de bandas
sen2 <- stack(brick_files)Crop stack
sen2 <- crop(sen2,xy_ext)
plot(sen2)Modificar nombres de las bandas
names(sen2) <- c('blue','green','red','RE1','RE2','RE3','SWIR1','SWIR2','NIR')Composiciones de color
plotRGB(sen2, r="red", g="green", b="blue", axes=TRUE, stretch ="lin", main ="True Color") # 4 3 2
plotRGB(sen2, r="NIR", g="red", b="green", axes=TRUE, stretch ="lin", main ="Falso Color Infrarojo") # 8 4 3
plotRGB(sen2, r="SWIR2", g="NIR", b="red", axes=TRUE, stretch ="lin", main ="SWIR") # 12 8 4
plotRGB(sen2, r="NIR", g="RE3", b="blue", axes=TRUE, stretch ="lin", main ="Vegetación") # 8 5 2(4)
plotRGB(sen2, r="SWIR2", g="SWIR1", b="blue", axes=TRUE, stretch ="lin", main ="Masas de Agua") # 12 11 2Interactuar con las imágenes desde R
En R se pueden desarrollar funciones que nos permitan interactuar con el mapa o imagen sin salir del entorno de programación. Para ello:
plot(sen2$red)Zoom
zoom(sen2$red) #Selección en mapaSeleccionar coordenadas
click() #Click mapaSeleccionar valores
click(sen2$red) #Click mapaSeleccionar extensión
drawExtent() #Selección en mapa
#plot(crop(sen2_crop$red,drawExtent()))Guardar resultados
En este paso aprenderemos cómo guardar las imágenes.
Crear directorio de resultados
path_results <- paste0(path_save,"RES/")
dir.create(path_results)Guardar archivo como TIF
writeRaster(subset(sen2,c("NIR","red","green")), filename=paste0(path_results,"FALSECOLOR.tif",sep=""),format="GTiff",datatype="FLT4S", overwrite=TRUE)Crear índices
Uno de los mayores potenciales de la teledetección es la posibilidad de crear índices a partir de imágenes que nos permitan analizar qué está ocurriendo sobre la superficie. A partir de aquí vamos a ver cómo se aplican algunos de los más utilizados.
NDVI (Normalized Difference Vegetation Index)
NDVI <- (sen2$NIR - sen2$red) / (sen2$NIR + sen2$red) # (B05 - B04) / (B05 + B04)
plot(NDVI)NDWI (Normalized Difference Water Index)
NDWI <- (sen2$green - sen2$NIR) / (sen2$green + sen2$NIR) # (B03 - B08) / (B03 + B08) McFeeters
plot(NDWI)Obtener valores superiores a un determinado valor (Filtrado de datos)
NDWI_filt <- NDWI > 0.3 #Detectar masas de agua
plot(NDWI_filt)MSI (Moisture Soil Index)
MSI <- sen2$SWIR2 / sen2$red # B11 / B08
plot(MSI)GNDVI (Green Normalized Difference Vegetation Index)
GNDVI <- (sen2$NIR - sen2$green) / (sen2$NIR + sen2$green) # (B08 - B03) / (B08 + B03)
plot(GNDVI)CVI – Chlorophyll Vegetation Index
CVI <- sen2$NIR * (sen2$red / (sen2$red^2)) # CVI
plot(CVI)Extracción de datos del raster
plotRGB(sen2, r="red", g="green", b="blue", axes=TRUE, stretch ="lin", main ="True Color") # 4 3 2Selección de coordenadas – AOI
c1 <- click() #Click mapa
plot(NDVI)
c2 <- click() #Click mapaJuntar coordenadas en objeto
xy <- rbind(c1,c2)Extracción de valores raster a partir de coordenadas
ext_sen <- extract(sen2, xy)Reordenar tabla por nombres
colnames(ext_sen) #Nombres de tabla
order_col <- c("blue","green","red","RE1","RE2","RE3","NIR","SWIR1","SWIR2") #Vector de nombres
ext_sen <- ext_sen[,order_col] #Reordenar
print(ext_sen)Crear gráfico de firmas espectrales
Crear gráfico vacío
plot(0, ylim=c(0, 10000), xlim = c(1,9), xlab = "Bandas", ylab = "Reflectancia", xaxt='n', main="Gráfico de firmas espectrales")Fijar nombres del eje X
axis(1, at=1:9, lab=order_col)Añadir información espectral
col = 1
for (i in 1:nrow(ext_sen)){
lines(ext_sen[i,], type = "o", lwd = 3, lty = 1, col=col, pch=1)
col = 1 + col
}Selección de coordenadas aleatorias
pol_ext <- as(extent(sen2), 'SpatialPolygons')
pts <- spsample(pol_ext[1,], 20, type = 'random')
plot(sen2$blue)<br>points(pts)Y esto ha sido todo. Puedes descargar el script completo desde aquí: https://bit.ly/3dSMkBx
Si quieres seguir profundizando en los conceptos de teledetección con R:
- https://custom-scripts.sentinel-hub.com/sentinel-2/composites/
- https://www.indexdatabase.de/db/s-single.php?id=96
Video del webinar:
Por último, os dejamos el vídeo completo del webinar.
