Introducción a R en el tratamiento de imágenes satélite

RemOT-Technologies lenguajes de programación R 10 junio, 2020 | 0

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.

Página web de Copernicus Open Access Hub. — Página web https://scihub.copernicus.eu

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 = 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 &lt;- as(extent(sen2), 'SpatialPolygons')
pts &lt;- 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.

lenguajes de programación r Sentinel-2

Introducción a R en el tratamiento de imágenes satélite