TIDYVERSE 2

Danilo Verdugo (dynageo@gmail.com)

RASTER: Sentinel-2

CDSE

API del ecosistema de datos de Copernicus

Este paquete proporciona la interfaz para la API del Ecosistema Espacial de Datos de Copernicus, principalmente para buscar en el catálogo de datos disponibles de las misiones Sentinel de Copernicus y obtener imágenes del área de interés según las bandas espectrales seleccionadas.

https://zivankaraman.github.io/CDSE/index.html

https://shapps.dataspace.copernicus.eu/dashboard/#/account/settings

Existe el paquete CDSE para trabajar con imágenes.

library(CDSE)
GetCollections(as_data_frame = TRUE)
                     id                   title
1        sentinel-2-l1c          Sentinel 2 L1C
2    sentinel-3-olci-l2      Sentinel 3 OLCI L2
3         landsat-ot-l1 Landsat 8-9 OLI-TIRS L1
4       sentinel-3-olci         Sentinel 3 OLCI
5      sentinel-3-slstr        Sentinel 3 SLSTR
6 sentinel-3-synergy-l2   Sentinel 3 Synergy L2
7        sentinel-1-grd          Sentinel 1 GRD
8        sentinel-2-l2a          Sentinel 2 L2A
9        sentinel-5p-l2    Sentinel 5 Precursor
                                                                   description
1                                     Sentinel 2 imagery processed to level 1C
2                 Sentinel 3 data derived from imagery captured by OLCI sensor
3                        Landsat 8-9 Collection 2 imagery processed to level 1
4                                   Sentinel 3 imagery captured by OLCI sensor
5                                  Sentinel 3 imagery captured by SLSTR sensor
6 Sentinel 3 data derived from imagery captured by both OLCI and SLSTR sensors
7                                     Sentinel 1 Ground Range Detected Imagery
8                                     Sentinel 2 imagery processed to level 2A
9                      Sentinel 5 Precursor imagery captured by TROPOMI sensor
                 since            instrument  gsd bands constellation long.min
1 2015-11-01T00:00:00Z                   msi   10    13    sentinel-2     -180
2 2016-04-17T11:33:13Z                  olci  300    NA          <NA>     -180
3 2013-03-08T00:00:00Z oli/tirs/oli-2/tirs-2   30    17          <NA>     -180
4 2016-04-17T11:33:13Z                  olci  300    21          <NA>     -180
5 2016-04-17T11:33:13Z                 slstr 1000    11          <NA>     -180
6 2016-04-17T11:33:13Z            olci/slstr  300    NA          <NA>     -180
7 2014-10-03T00:00:00Z                 c-sar   NA    NA    sentinel-1     -180
8 2016-11-01T00:00:00Z                   msi   10    12    sentinel-2     -180
9 2018-04-30T00:18:50Z               tropomi 5500    NA          <NA>     -180
  lat.min long.max lat.max
1     -56      180      83
2     -85      180      85
3     -85      180      85
4     -85      180      85
5     -85      180      85
6     -85      180      85
7     -85      180      85
8     -56      180      83
9     -85      180      85

Buscar en la página epsg, la versión wkt de la proyección UTM H19S para evitar algunos mensajes de alerta por versión de cambio de librería PROJ.

coor = "+proj=utm +zone=19 +south +ellps=GRS80 +units=m +no_defs +type=crs"

Pasos para descargar límites, seleccionar y proyectar mi polígono de trabajo.

mapa_crudo <- gadm(country = "CHL", 
                   level = 3, 
                   path = tempdir()) 

olivar = st_as_sf(mapa_crudo) |>
  filter(NAME_3 == "Olivar") |>
  st_transform(crs=coor)

plot(olivar$geometry)
#Extracción DTM
mapa_dtm = get_elev_raster(locations = olivar, 
                           z = 11, 
                           src = "aws", 
                           clip = "locations", 
                           tmp_dir = tempdir(), 
                           ncpu = ifelse(future::availableCores() > 2, 2, 1)) |>
  rast()

plot(mapa_dtm)
oclient = GetOAuthClient(id = 'sh-0cb438ad-5ed4-4419-b3ef-c86e896049ed',
                         secret = 'wyovOi29kgwKiCiNYuiTdYFkTqNb9vum')
#'sh-d484c1de-042d-4910-b7d3-1b2afd85c537'
#'NYiVMEkLdOHXkairtDhqZs2L1yKIDkSq'

Datos que vamos a utilizar:

https://drive.google.com/drive/folders/1wv5MpMv8BEcBy80ue80C2eSA29kBBgQw?usp=drive_link

raw_script_text <- paste(readLines("data/UCIBands.js"), collapse = "\n")
images <- SearchCatalog(aoi = olivar, from = "2025-12-01", to = "2025-12-15",
                        collection = "sentinel-2-l2a", with_geometry = TRUE,
                        client = oclient)

Algunos métodos para explorar los resultados, el objetivo es encontrar la fecha de la mejor imagen:

range(images$areaCoverage)
[1] 100 100
range(images$tileCloudCover)
[1]  0.09 58.90
img = images[images$tileCloudCover<1,]

day = img$acquisitionDate[1]

Paso 4: Ahora con la fecha y el área descargo mi imagen a la resolución de 10m.

ras <- GetImage(aoi = olivar, time_range = day, script = raw_script_text,
                       collection = "sentinel-2-l2a", format = "image/tiff",
                       mosaicking_order = "mostRecent", resolution = 10,
                       mask = T, buffer = 1, client = oclient )

plot(ras)
#NO OLVIDAR QUE PUEDES EXPORTAR
#writeRaster(ras, "data/mapa_tmp.tif", overwrite = T, NAflag=0)

SCL

Puede ser útil descargar la clasificación de uso de suelo automática, permite enmascarar valores anómalos.

#script para descargar.
raw_scl_text <- paste(readLines("data/SCLBand.js"), collapse = "\n")

#usamos la misma función.
scl <- GetImage(aoi = olivar, time_range = day, script = raw_scl_text,
                       collection = "sentinel-2-l2a", format = "image/tiff",
                       mosaicking_order = "mostRecent", resolution = 10,
                       mask = T, buffer = 1, client = oclient)

#forma de crear un mapa clasificado de acuerdo a los valores deseados.
mask = ifel(scl %in% c(0,7,8,9,10,11), NA,1)

plot(mask)

#NO OLVIDAR QUE PUEDES EXPORTAR
#writeRaster(scl, "data/mapa_scl.tif", overwrite = T, NAflag=0)

Operaciones

Preparando nuestra imagen.

Renombrar las bandas.

#Pre ordenamiento- colocando nombres
nombres = c("Aero","Azul","Verde","Rojo","RedEdge1","RedEdge2","RedEdge3","NIR","RedEdge4","H2O","SWIR1","SWIR2")
names(ras) = nombres

Combinación de bandas.

Color verdadero RGB.

#plotear las bandas
plotRGB(ras, r="Rojo", g="Verde", b="Azul", stretch="lin")

Falso color IR-R-G

plotRGB(ras,r="NIR",g="Rojo",b="Verde",stretch="lin")

Falso color para uso suelos.

plotRGB(ras,r="SWIR2",g="SWIR1",b="NIR",stretch="lin")

Operación de Bandas

Razón de Bandas

#Razón de Sultán
r = ras$SWIR2 / ras$Azul
g = ras$SWIR2 / ras$SWIR1
b = (ras$SWIR2 / ras$NIR) * (ras$Rojo / ras$NIR)

sultan= c(r,g,b)
names(sultan) = c("R", "G","B")
plotRGB(sultan,r=1,g=2,b=3,stretch="lin")

Índices

NDVI

ndvi = (ras$NIR - ras$Rojo)/(ras$NIR+ras$Rojo)
plot(ndvi)

NDSI

ndsi = (ras$Verde-ras$SWIR1)/(ras$Verde+ras$SWIR1)
plot(ndsi)

Filtrado

#Filtrar
vegetacion = ifel(ndvi > .75, 1, NA)
plot(vegetacion)

Conversión

#Vectorizar y exportar
tmp_pol=as.polygons(vegetacion)
plot(tmp_pol)

Escritura

note que debemos convertir a objeto ‘sf’ con st_as_sf().

st_write(st_as_sf(tmp_pol),"data/test1.shp", append=F)