Funciones externas

Las funciones para estimar la línea separador entre día y noche se basan en un javascript L.Terminator.js del paquete {Leaflet} que encontré en stackoverflow. El script con las funciones lo podéis descargar aquí o acceder en github.

source("terminator.R") # importamos las funciones

Funciones personalizadas

La función principal terminator() basada en el javascript de {Leaflet} necesita como argumentos: la fecha con la hora, la extensión mínima y máxima de longitud así como la resolución o el intervalo de longitud.

t0 <- Sys.time() # fecha y hora de nuestro sistema operativo
t0

## [1] "2022-03-27 13:10:39 CEST"

coord_nightday <- terminator(t0, -180, 180, 0.2) # estimamos la línea día-noche

# lo convertimos en un objecto espacial de clase sf
line_nightday <- st_linestring(as.matrix(coord_nightday)) %>% st_sfc(crs = 4326) 

# ploteamos
plot(line_nightday)

En el sigiente paso obtenemos los polígonos que corresponden al día y la noche que separa la línea estimada anteriormente. Para ello creamos un rectángulo cubriendo todo el planeta y empleamos la función st_split() del paquete {lwgeom} que divide el rectángulo.

# rectángulo 
wld_bbx <- st_bbox(c(xmin = -180, xmax = 180,
                       ymin = -90, ymax = 90), 
                     crs = 4326) %>% 
             st_as_sfc()

# divisón con la línea día-noche
poly_nightday <-  st_split(wld_bbx, line_nightday) %>% 
                      st_collection_extract(c("POLYGON")) %>% 
                       st_sf() 

# ploteamos 
plot(poly_nightday)

La pregunta que ahora surge es cuál de los dos poligonos corresponde a la noche y cúal al día. Eso dependerá en qué día del año estamos, dado los cambios de posición de la Tierra con respecto al Sol. Entre el equinoccio de primeravera y el equinoccio de otoño corresponde con el primer polígono, cuándo también podemos observar el día polar en el polo norte, y en el caso contrario sería el segundo. El paquete {terra} sólo acepta la clase vectorial propia llamada SpatVector, por eso convertimos el objeto vectorial sf con la función vect().

# selecionamos el segundo polígono
poly_nightday <- slice(poly_nightday, 2) %>% 
                    mutate(daynight = 1)

# creamos el raster con una resolución de 0,5º y la extensión del mundo
r <- rast(vect(wld_bbx), resolution = .5)

# rasterizamos el polígono de noche 
night_rast <- rasterize(vect(poly_nightday), r) 

# resultado en formato raster
plot(night_rast)

En el último paso reproyectamos el raster a Mollweide.

# definimos la proyección del raster (WGS84)
crs(night_rast) <- "EPSG:4326"

# reproyectamos
night_rast_prj <- project(night_rast, "ESRI:54009", 
                          mask = TRUE, 
                          method = "near")
# mapa
plot(night_rast_prj)

Finalmente incluimos los pasos individuales que hemos hecho en una función personalizada.

rast_determiner <- function(x_min, date, res) {
  
  # crea fecha con hora añadiendo el número de minutos
  t0 <- as_date(date) + minutes(x_min) 
  # estimamos las coordenadas de la línea que separa día y noche 
  night_step <- terminator(t0, -180, 180, 0.2) %>% as.matrix()
  # pasamos los puntos a línea
  night_line <- st_linestring(night_step) %>% st_sfc(crs = 4326)
  
  # definimos el rectángulo del planeta
  wld_bbx <- st_bbox(c(xmin = -180, xmax = 180,
                       ymin = -90, ymax = 90), 
                     crs = 4326) %>% 
             st_as_sfc()
  
  # dividimos el polígono con la línea de día-noche
  poly_nightday <-  st_split(wld_bbx, night_line) %>% 
                      st_collection_extract(c("POLYGON")) %>% 
                       st_sf()  
  
  # seleccionamos el polígono según la fecha
  if(date <= make_date(year(date), 3, 20) | date >= make_date(year(date), 9, 23)) {
    
    poly_nightday <- slice(poly_nightday, 2) %>% 
      mutate(daynight = 1)
    
  } else {
    
    poly_nightday <- slice(poly_nightday, 1) %>% 
      mutate(daynight = 1)
  }
  
  # creamos el raster con la resolución del argumento res
  r <- rast(vect(wld_bbx), resolution = res)
  
  # rasterizamos el polígono de noche
  night_rast <- rasterize(vect(poly_nightday), r) 
  
  return(night_rast)
  
}

Dado que queremos obtener el área de noche para diferentes horas del día construimos una segunda función para aplicar la primera sobre diferentes intervalos del día (en minutos).

night_determinator <- function(time_seq, # minutos 
                               date = Sys.Date(), # fecha (por defecto del sistema)
                               res = .5) { # resolución del raster 0.5º

# aplicamos la primera función sobre un vector de minutos
night_raster <-  map(time_seq, 
                     rast_determiner,
                     date = date, 
                     res = res)

# convertimos el raster en un objeto de tantas capas como unidades de minutos
night_raster <- rast(night_raster)

# definimos la proyección WGS84
crs(night_raster) <- "EPSG:4326"

return(night_raster)
  
}

Preparación

Para crear una animación de 24 horas mostrando el movimiento de la noche sobre la Tierra debemos hacer unos pasos previos. Primero obtenemos los límites del mundo con la función ne_countries() y los reproyectamos a la nueva proyección Winkel II. Después convertimos los datos raster en un data.frame indicando que mantenga valores ausentes. Podemos observar que cada capa del raster (de cada intervalo de 30 minutos) es una columna en el data.frame. Renombramos las columnas y convertimos la tabla en un formato largo empleando la función pivot_longer(). Lo que hacemos es fusionar todas las columnas de las capas en una única. Como último paso excluimos los valores ausentes con la función filter().

# límites de países
wld <- ne_countries(scale = 10, returnclass = "sf") %>% 
         st_transform("ESRI:54019")

# convertimos el raster a un data.frame con xyz
df_winkel <- as.data.frame(night_raster_winkel, xy = TRUE, na.rm = FALSE)

# renombramos todas las columnas correspondientes a los intervalos del día
names(df_winkel)[3:length(df_winkel)] <- str_c("H", as_hms(seq(0, 1410, 30)*60))

# cambiamos a un formato largo de tabla
df_winkel <- pivot_longer(df_winkel, 3:length(df_winkel), names_to = "hour", values_to = "night") 

# excluimos los valores ausentes para reducir el tamaño de la tabla
df_winkel <- filter(df_winkel, !is.na(night))

Sólo resta crear una retícula y obtener la extensión del mapamundi.

# retícula
grid <- st_graticule() %>%   st_transform("ESRI:54019")

# obtenemos la extensión del mundo
bbx <- st_bbox(wld)

El mapa de cualquier hora lo construimos con ggplot2 añadiendo la geometría vectorial con geom_sf() (los límites y la retícula) y con geom_raster() los datos raster. En el título estamos usando un símbolo unicode como reloj. Además definimos la extensión del mapa en coord_sf() para mantenerlo constante sobre todos los mapas en la animación. Por último, hacemos uso de {{ }} de {rlang} dentro de la función filter() para poder filtrar nuestros datos raster en forma de tabla. Con el objetivo de que nuestra función pueda evaluar correctamente los valores que pasamos en x (los intervalos del día) es necesario usar esta gramatica de tidy evaluation por data masking de tidyverse. Es un tema para otro post.

# ejemplo 5 UTC
x <- "H05:00:00"

# mapa
ggplot() +
  # límites
  geom_sf(data = wld,
        fill = "#74a9cf", 
        colour = "white",
          size = .1) +
  # retícula
  geom_sf(data = grid, size = .1) +
  # datos raster filtrados 
  geom_raster(data = filter(df_winkel, hour == {{x}}), 
              aes(x, y), 
            fill = "grey90",
            alpha = .6) +
  # título
  labs(title = str_c("\U1F551", str_remove(x, "H"), " UTC")) + 
  # límites de extensión
  coord_sf(xlim = bbx[c(1, 3)], 
           ylim = bbx[c(2, 4)])  +
  # estilo del mapa
  theme_void() +
  theme(plot.title = element_text(hjust = .1, vjust = .9))

## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.

Animación

Creamos la animación aplicando la función walk(), que a su vez pasará por el vector de intervalos para filtrar nuestros datos y mapear cada paso haciendo uso de ggplot.

walk(str_c("H", as_hms(seq(0, 1410, 30)*60)), function(step){
  
g <- ggplot() +
  geom_sf(data = wld,
        fill = "#74a9cf", 
        colour = "white",
          size = .1) +
  geom_sf(data = grid,
          size = .1) +
  geom_raster(data = filter(df_winkel, hour == {{step}}), 
              aes(x, y), 
            fill = "grey90",
            alpha = .6) +
  labs(title = str_c("\U1F551", str_remove(x, "H"), " UTC")) + 
  coord_sf(xlim = bbx[c(1, 3)], ylim = bbx[c(2, 4)])  +
  theme_void() +
  theme(plot.title = element_text(hjust = .1, vjust = .9))


ggsave(str_c("wld_night_", str_remove_all(step, ":"), ".png"), g,
       height = 4.3, width = 8.4, bg = "white", dpi = 300, units = "in")

})

La creación del gif final lo hacemos con gifski() pasándole los nombres de las imagenes en el orden como deben aparecer en la animación.

files <- str_c("wld_night_H", str_remove_all(as_hms(seq(0, 710, 30)*60), ":"), ".png")

gifski(files, "night_day.gif", width = 807, height = 409, loop = TRUE, delay = 0.1)