C/ Mariano Esquillor, s/n. 50018 - Zaragoza
+34 619 12 46 15
info@remot-technologies.com

Webinar R

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.

Página web https://scihub.copernicus.eu

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.

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)

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:

Video del webinar:

Aquí os dejamos el vídeo completo del webinar.

A %d blogueros les gusta esto: