Bienvenido al webinar de introducción a R organizado en conjunto con nuestros amigos de Geoinnova.
En este webinar aprenderás los conceptos básicos para iniciarte en R y usarlo para el procesado de imágenes de teledetección.
Qué veremos en el webinar
- 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
Debajo tienes el link (expirará pasados unos días) con la imagen que usaremos durante el webinar.
Imagen Sentinel-2: https://we.tl/t-HY8CinTSsv
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 tu propio script de R y puedas seguir el webinar.
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 = 20
Comprobar 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ía
Colorear 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-y
Guardar valores del raster como un objeto
b02_hist <- getValues(b02)
Crear el histograma
hist(b02_hist, breaks=5000, xlim=c(0,3000)) #Histograma
Recortar 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 2
Interactuar 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 mapa
Seleccionar coordenadas
click() #Click mapa
Seleccionar valores
click(sen2$red) #Click mapa
Seleccionar 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 2
Selección de coordenadas – AOI
c1 <- click() #Click mapa
plot(NDVI)
c2 <- click() #Click mapa
Juntar 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)
Llegamos al final del webinar. ¡Esperamos que os haya gustado tanto como a nosotros prepararlo!
Puedes descargar el script completo desde aquí: https://we.tl/t-Cnk7vIXWEL (El link expirará en unos días)
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:
Aquí os dejamos el vídeo completo del webinar.