2.1 Antes de usar GEE en R

El primer paso consiste en crear nuestro usuario en earthengine.google.com. Además, es necesario instalar CLI de gcloud (https://cloud.google.com/sdk/docs/install?hl=es-419), únicamente debes seguir las instrucciones en Google. Con respecto al lenguaje de GEE, muchas funciones que se aplican son similares a lo que se conoce de {tidyverse}. Se puede obtener más ayuda en https://r-spatial.github.io/rgee/reference/rgee-package.html y en la propia página de GEE.

Lo más esencial del lenguaje nativo Javascript de GEE es que se caracteriza por la forma de combinación de funciones y variables usando el punto, el que se sustuye por el $ en R. Todas las funciones GEE empiezan por el prefijo ee_* (ee_print(), ee_image_to_drive()).

Una vez que tenemos instalado gcloud y el paquete {rgee} podemos proceder a crear el entorno virtual de Python. La función ee_install() se encarga de instalar Anaconda3 y todos los paquetes necesarios. Para comprobar la correcta instalación de Python, y particularmente de los paquetes numpy y earthengine-api, podemos usar ee_check().

ee_install() # crear entorno virtual de Python
ee_check() # comprobar si todo está correcto

Antes de pasar a programar con la sintaxis propia de GEE se debe autenticar e inicializar GEE empleando la función ee_Initialize().

ee_Initialize(drive = TRUE) # autenticar e inicializar GEE

## ── rgee 1.1.5 ─────────────────────────────────────── earthengine-api 0.1.339 ── 
##  ✔ user: not_defined
##  ✔ Google Drive credentials:
 ✔ Google Drive credentials:  FOUND
##  ✔ Initializing Google Earth Engine:
 ✔ Initializing Google Earth Engine:  DONE!
## 
 ✔ Earth Engine account: users/dominicroye 
## ────────────────────────────────────────────────────────────────────────────────

2.2 Acceder a la predicción del GFS

En GEE se denomina ImageCollection una serie temporal de imágenes o datos multidimensionales. Cada dataset tiene asignado una ID y podemos acceder a ella haciendo la llamada ee$ImageCollection('ID_IMAGECOLLECTION'). Hay funciones auxilares que permiten la conversión de clases puramente de R a Javascript, por ejemplo, la fecha rdate_to_eedate(). Lo primero lo que hacemos es filtrar a la fecha más reciente con las últimas ejecuciones del modelo GFS.

Debemos saber que, a diferencia de R, únicamente cuando se envían tareas GEE se ejecuta el cálculo en los servidores, enviando todos los objetos creados. En la mayoría de los pasos se crean objetos EarthEngine.

## GFS forcast
dataset <- ee$ImageCollection('NOAA/GFS0P25')$filter(ee$Filter$date(rdate_to_eedate(today()-days(1)),
                                                                    rdate_to_eedate(today()+days(1))))
dataset

## EarthEngine Object: ImageCollection

Las ejecuciones del modelo se realizan cada 6 horas (0, 6, 12, 18), así que la función ee_get_date_ic() nos extrae las fechas para poder elegir la más reciente. Aquí es la primera vez que se ejecutan cálculos.

# vector de fechas únicas de ejecuciones
last_run <- ee_get_date_ic(dataset)$time_start |> unique()
last_run

## [1] "2023-02-22 00:00:00 GMT" "2023-02-22 06:00:00 GMT"
## [3] "2023-02-22 12:00:00 GMT" "2023-02-22 18:00:00 GMT"
## [5] "2023-02-23 00:00:00 GMT"

# seleccionamos la última
last_run <- max(last_run)

A continuación filtramos la fecha de última ejecución y seleccionamos la banda de la temperatura de aire a 2m. Las fechas de la predicción futura son atributos de cada ejecución con una predicción de hasta 336 horas (14 días) desde el día de la misma. Cuando queramos hacer cambios a cada imagen de una ImageCollection debemos hacer uso de la función map(), similar a la que conocemos del paquete {purrr}. En este caso redefinimos la fecha de cada imagen (system:time_start: fecha de ejecución) por la de la predicción (forecast_time). Es importante que la función a aplicar de R se encuentre dentro de ee_utils_pyfunc(), la que la traduce a Python. Después extraemos las fechas de la predicción.

# última ejecución y selección de variable 
temp <- dataset$filter(ee$Filter$date(rdate_to_eedate(last_run)))$select('temperature_2m_above_ground')

# definimos las fechas de predicción para cada hora
forcast_time <- temp$map(ee_utils_pyfunc(function(img)  {
  
 return(ee$Image(img)$set('system:time_start',ee$Image(img)$get("forecast_time")))

  })
)

# obtenemos las fechas de predicción
date_forcast <- ee_get_date_ic(forcast_time)
head(date_forcast)

##                            id          time_start
## 1 NOAA/GFS0P25/2023022300F000 2023-02-23 00:00:00
## 2 NOAA/GFS0P25/2023022300F001 2023-02-23 01:00:00
## 3 NOAA/GFS0P25/2023022300F002 2023-02-23 02:00:00
## 4 NOAA/GFS0P25/2023022300F003 2023-02-23 03:00:00
## 5 NOAA/GFS0P25/2023022300F004 2023-02-23 04:00:00
## 6 NOAA/GFS0P25/2023022300F005 2023-02-23 05:00:00

Aquí podríamos exportar los datos horarios de temperatura, pero también sería posible estimar la temperatura máxima o mínima diaria de los próximos 14 días. Para lograrlo debemos definir el inicio y final del periodo, y calcular el número de días. Lo que hacemos en términos simples es mapear sobre el número días para filtrar a cada día y aplicar la función max() o cualquier otra similar.

# definimos inicio final de las fechas
endDate <- rdate_to_eedate(round_date(max(date_forcast$time_start)-days(1), "day"))
startDate <- rdate_to_eedate(round_date(min(date_forcast$time_start), "day"))

# número de dias
numberOfDays <- endDate$difference(startDate, 'days')

# calculamos la máxima diaria 
daily <- ee$ImageCollection(
  ee$List$sequence(0, numberOfDays$subtract(1))$
  map(ee_utils_pyfunc(function (dayOffset) {
    start = startDate$advance(dayOffset, 'days')
    end = start$advance(1, 'days')
    return(forcast_time$
    filterDate(start, end)$
    max()$ # alternativa: min(), mean()
    set('system:time_start', start$millis()))
  }))
)

# fechas de la máxima diaria
head(ee_get_date_ic(daily))

##      id time_start
## 1 no_id 2023-02-23
## 2 no_id 2023-02-24
## 3 no_id 2023-02-25
## 4 no_id 2023-02-26
## 5 no_id 2023-02-27
## 6 no_id 2023-02-28

2.3 Mapa dinámico vía GEE

Dado que existe la posibilidad de añadir imágenes a un mapa dinámico en el editor de código de GEE, también podemos hacerlo usando desde R la función GEE Map.addLayer(). Seleccionamos simplemente el primer día con first(). En los otros argumento definimos el rango de los valores y la gama de colores.

Map$addLayer(
       eeObject = daily$first(),
       visParams = list(min = -45, max = 45,
                        palette = rev(RColorBrewer::brewer.pal(11, "RdBu"))),
       name = "GFS") + 
Map$addLegend(
  list(min = -45, max = 45, 
       palette = rev(RColorBrewer::brewer.pal(11, "RdBu"))), 
       name = "Temperatura máxima", 
       position = "bottomright", 
       bins = 10)

2.4 Exportar múltiples imagenes

El paquete {rgee} tiene una función muy útil para exportar una ImageCollection: ee_imagecollection_to_local(). Antes de usarla debemos fijar una región, la que se pretende exportar. En este caso, exportamos todo el globo con un rectángulo cubriendo la Tierra.

# extensión de la Tierra
geom <- ee$Geometry$Polygon(coords = list(
  c(-180, -90), 
  c(180, -90),
  c(180, 90),
  c(-180, 90),
  c(-180, -90)
),
proj = "EPSG:4326",
geodesic = FALSE)

geom #vemos que es un objeto EarthEngine de tipo geometría

## EarthEngine Object: Geometry

# carpeta temporal de descarga
tmp <- tempdir()

# descargamos 
ic_drive_files_2 <- ee_imagecollection_to_local(
  ic = daily$filter(ee$Filter$date(rdate_to_eedate(today()), rdate_to_eedate(today()+days(2)))), # elegimos los próximos 4 días
  region = geom,
  scale = 20000, 
  lazy = FALSE,
  dsn = path(tmp, "drive_"),
  add_metadata = TRUE
)

## ───────────────────────────────────── Downloading ImageCollection - via drive ──- region parameters
##  sfg      : POLYGON ((-180 -90, 180 -90, 180 90, -180 90, -180 -90)) 
##  CRS      : GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563, ..... 
##  geodesic : FALSE 
##  evenOdd  : TRUE 
## 
## Downloading: C:/Users/xeo19/AppData/Local/Temp/RtmpE9IUrn/drive_0.tif
## Downloading: C:/Users/xeo19/AppData/Local/Temp/RtmpE9IUrn/drive_1.tif
##  ────────────────────────────────────────────────────────────────────────────────

3.1 Datos

El primer paso consiste en importar los datos con ayuda de rast(). Además definimos el nombre de cada capa según su dimensión temporal correctamente.

# rutas a los datos descargados
forecast_world <- dir_ls(tmp, regexp = "tif")

# garantizamos el orden y importamos
file_ord <- str_extract(forecast_world, "_[0-9]{1,2}") |> parse_number()

forecast_rast <- rast(forecast_world[order(file_ord)])
forecast_rast

## class       : SpatRaster 
## dimensions  : 1002, 2004, 2  (nrow, ncol, nlyr)
## resolution  : 0.1796631, 0.1796631  (x, y)
## extent      : -180.0224, 180.0224, -90.01119, 90.01119  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : drive_0.tif  
##               drive_1.tif  
## names       : temperature_2m_above_ground, temperature_2m_above_ground

# definimos la dimensión temporal como nombre de cada capa
names(forecast_rast) <- seq(today(), today() + days(1), "day")
forecast_rast

## class       : SpatRaster 
## dimensions  : 1002, 2004, 2  (nrow, ncol, nlyr)
## resolution  : 0.1796631, 0.1796631  (x, y)
## extent      : -180.0224, 180.0224, -90.01119, 90.01119  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : drive_0.tif  
##               drive_1.tif  
## names       : 2023-02-23, 2023-02-24

# ploteamos
plot(forecast_rast)

Ahora debemos definir la proyección ortográfica indicando con +lat_0 y +lon_0 el centro de la proyección. Después reproyectamos y convertimos el ráster en un data.frame.

# definición de la proyección 
ortho_crs <-'+proj=ortho +lat_0=51 +lon_0=0.5 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs'

# reproyectamos el ráster
ras_ortho <- project(forecast_rast, ortho_crs)

# convertimos el ráster en un data.frame de xyz
forecast_df <- as.data.frame(ras_ortho, xy = TRUE)

# transformamos a un formato largo
forecast_df <- pivot_longer(forecast_df, 3:length(forecast_df), names_to = "date", values_to = "ta")

3.2 Límites administrativos y retícula

Importamos los límites administrativos con gisco_get_countries() los que debemos preparar para la proyección ortográfica. Del mismo modo creamos la retícula empleando st_graticule(). Con el objetivo de preservar la geometría habrá que recortar a únicamente la parte visible. El océano se crea partiendo de un punto en 0,0 con el radio de la tierra. Empleando la función st_intersection() reducimos a la parte visible y reproyectamos los límites.

# obtenemos los límites administrativos
world_poly <- gisco_get_countries(year = "2016", epsg = "4326", resolution = "10") 

# obtenemos la retícula global
grid <- st_graticule()

# definimos lo que sería oceano 
ocean <- st_point(x = c(0,0)) |>
            st_buffer(dist = 6371000) |> # radio de la tierra
              st_sfc(crs = ortho_crs)
plot(ocean)

# seleccionamos sólo visible de los límites y reproyectamos
world <- world_poly |>
            st_intersection(st_transform(ocean, 4326)) |>
            st_transform(crs = ortho_crs) # 
plot(world)

Debemos repetir la misma selección visible para la retícula, aunque aquí antes limitamos las líneas de la retícula a los océanos. El límite del océano nos sirve para crear la sombra del globo, pero para usarlo en geom_glowpath() es necesario convertirlo en un data.frame.

# eliminamos las líneas que transcurren sobre los continentes
grid_crp <- st_difference(grid, st_union(world_poly))

# seleccionamos la parte visible
grid_crp <- st_intersection(grid_crp, st_transform(ocean, 4326)) |>
                  st_transform(crs = ortho_crs)

plot(grid_crp)

# convertimos el límite del globo en un data.frame
ocean_df <- st_cast(ocean, "LINESTRING") |> st_coordinates() |> as.data.frame()

3.3 Construcción del mapa

forecast_tomorrow <- filter(forecast_df, date == today()+days(1)) |>
                        mutate(ta_limit = case_when(ta > 45 ~ 45,
                                              ta < -45 ~ -45,
                                               TRUE ~ ta))

# construimos una sombra 
ggplot() + 
   geom_glowpath(data = ocean_df, 
                aes(X, Y, group = "L1"),
                shadowcolor='grey90',
                     colour = "white",
                alpha = .01,
                shadowalpha=0.05,
                shadowsize = 1.5) +
   coord_sf() +
   theme_void()

# combinando varias capas de sobra
g <- ggplot() +
   geom_glowpath(data = ocean_df, 
                aes(X, Y, group = "L1"),
                shadowcolor='grey90',
                     colour = "white",
                alpha = .01,
                shadowalpha=0.05,
                shadowsize = 1.8) +
   geom_glowpath(data = ocean_df, 
                aes(X, Y, group = "L1"),
                shadowcolor='grey90',
                       colour = "white",
                alpha = .01,
                shadowalpha=0.02,
                shadowsize = 1) +
   geom_glowpath(data = ocean_df, 
                aes(X, Y, group = "L1"),
                shadowcolor='grey90',
                       colour = "white",
                alpha = .01,
                shadowalpha=0.01,
                shadowsize = .5) 

g2 <- g + geom_raster(data = forecast_tomorrow, aes(x, y, fill = ta_limit)) +
          geom_sf(data = grid_crp, 
                  colour = "white", 
                  linewidth = .2) +
          geom_sf(data = world, 
                   fill = NA,
                   colour = "grey10",
                   linewidth = .2) 

g2 + scale_fill_distiller(palette = "RdBu", 
                          limits = c(-45, 45),
                          breaks = c(-45, -25, 0, 25, 45)) +
     guides(fill = guide_colourbar(barwidth = 15, 
                                   barheight = .5, 
                                   title.position = "top",
                                   title.hjust = .5)) +
  coord_sf() +
  labs(fill = str_wrap("Temperatura máxima a 2 metros para el 14 de febrero", 30)) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 7),
        plot.margin = margin(10, 10, 10, 10)) 

labeling <- slice(forecast_tomorrow, which.min(ta), which.max(ta))
labeling

## # A tibble: 2 × 5
##           x         y date          ta ta_limit
##       <dbl>     <dbl> <chr>      <dbl>    <dbl>
## 1 -1397101.  3923352. 2023-02-24 -42.5    -42.5
## 2  2119554. -4121598. 2023-02-24  43.9     43.9

g2 +  geom_mark_circle(data = labeling, 
                       aes(x, y, 
                          description = str_glue('{scales::number(ta, accuracy = .1, decimal.mark = ",", style_positive = "plus", suffix = "ºC")}')
                          ), 
                   expand = unit(1, "mm"), 
                   label.buffer = unit(4, "mm"),
                   label.margin = margin(1, 1, 1, 1, "mm"),
                   con.size = 0.3,
                   label.fontsize = 8,
                   label.fontface = "bold",
                   con.type = "straight",
                  label.fill = alpha("white", .5)) +
     scale_fill_distiller(palette = "RdBu", 
                          limits = c(-45, 45),
                          breaks = c(-45, -25, 0, 25, 45)) +
     guides(fill = guide_colourbar(barwidth = 15, 
                                   barheight = .5, 
                                   title.position = "top",
                                   title.hjust = .5)) +
  coord_sf(crs = ortho_crs) +
  labs(fill = str_wrap("Temperatura mínima 2 metros para el 14 de febrero", 30)) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 7),
        plot.margin = margin(10, 10, 10, 10))