Movimiento de nubes (precipitaciones) en la Península Ibérica

Cargamos y trabajamos los datos

Los datos con los que vamos a trabajar han sido obtenidos del servidor THREDSS de Meteogalicia el cuál proporciona acceso a través de diferentes protocolos a variables de superficie del modelo Weather Research Forecast (WRF), un sistema numérico de predicción meteoróligica mesoescalar. Servidor THREDDS de MeteGalicia: https://www.meteogalicia.gal/web/modelos/threddsIndex.action?request_locale=es&msclkid=998d2299cfad11ecbdb8f8e59fac9015

“El análisis mesoescala utiliza las herramientas de simulación atmosférica más avanzadas para obtener una serie temporal de datos meteorológicos.

Las series temporales mesoescala de datos meteorológicos se generan a partir de un sistema numérico de predicción meteorológica, WRF (Weather Research and Forecasting), ejecutado en un clúster de alto rendimiento.”

• http://www.normawind.com/servicios-detalle.php?idservicio=7

El modelo WRF se ejecuta operativamente dos veces al día iniciando a las 00h y a las 12h. El primero se ejecuta para un total de 96 horas y el segundo para un total de 84 horas. Asimismo, hay configurados tres dominios para resoluciones de 36km, 12km y 4km.

Del conjunto de variables disponibles, vamos a utilizar la predicción horaria de nubosidad a niveles medios y bajos (precipitaciones). Para ello, vamos a seleccionar los datos con un horizonte temporal de 96 horas y resolución espacial de 12km que se corresponden con el día 28/05/2022. Link para descargar los datos: http://mandeo.meteogalicia.es/thredds/catalog/wrf_2d_12km/fmrc/files/catalog.html

Los datos vienen en formato netCDF y son de clase “RasterBrick”. Usamos la función brick del paquete {raster} para cargarlos:

cft <- brick('wrf_arw_det_history_d02_20220527_0000.nc4')
cft

## class      : RasterBrick 
## dimensions : 138, 171, 23598, 96  (nrow, ncol, ncell, nlayers)
## resolution : 12, 12  (x, y)
## extent     : -6, 2046, -6, 1650  (xmin, xmax, ymin, ymax)
## crs        : +proj=lcc +lat_0=34.8230018615723 +lon_0=-14.1000003814697 +lat_1=43 +lat_2=43 +x_0=536.40234 +y_0=-18.55861 +datum=WGS84 +units=m +no_defs 
## source     : wrf_arw_det_history_d02_20220527_0000.nc4 
## names      : X3600, X7200, X10800, X14400, X18000, X21600, X25200, X28800, X32400, X36000, X39600, X43200, X46800, X50400, X54000, ... 
## z-value    : 3600, 345600 (min, max)
## varname    : prec

Con la función getValues obtenemos los valores. Al ser un objeto RasterBrick, los valores que devuelve esta función son una matriz en la que las filas representan las celdas y las columnas las capas.

cft[] <- getValues(cft)

A continuación fijamos la proyección, y la asignamos a nuestro objeto cft.

projLCC2d <- "+proj=lcc +lon_0=-14.1 +lat_0=34.823 +lat_1=43 +lat_2=43 +x_0=536402.3 +y_0=-18558.61 +units=km +ellps=WGS84"
projection(cft) <- projLCC2d

Ahora fijamos el índice de tiempo. Haciendo uso de la función as.POSIXct creamos una variable que se llama timeIndex en la que almacenamos una secuencia con la fecha y la hora utilizando la zona horaria UTC. Dicha secuencia tiene una longitud de 96 (que es el número de capas que contiene nuestro objeto cft) e incrementa por horas.

Añadimos esta información temporal asociada con las capas de nuestro objeto a este con la función setZ. Finalmente, usamos la función format para que el nombre de las capas de nuestro objeto se correspondan con su fecha y hora pero en un formato de manera que se muestre únicamente el día detrás de la letra D, y la hora detrás de la letra H.

##set time index
timeIndex <- seq(as.POSIXct('2022-05-27 01:00:00', tz = 'UTC'), length = 96, by = 'hour')
cft <- setZ(cft, timeIndex)
names(cft) <- format(timeIndex, 'D%d_H%H')

Damos contexto espacial

Para dar contexto espacial haremos uso del paquete map. Este paquete trabaja con una projección longitud-latitud y nuestros datos vienen en projección raster. Por lo que necesitamos hacer una trasnformación.

Para ello, creamos un objeto llamado projLL al cual le asignamos un sistema de coordenadas de referencia (CRS) que se corresponde con el WGS 84, un sistema de coordenadas geográficas usado munidalmente que permite localizar cualquier punto de la Tierra. Utilizamos la proyección longitud-latitud, el “datum” WGS84 (el datum describe el punto vertical y horizontal de referencia del sistema de coordenadas), el elipsoide (la manera en la que la redondez de la Tierra está calculada) para los datos es el WGS84 y finalmente, el parámetro “towgs84” describe los parámetros de conversión lo cuales, en este caso son 3: las coordenadas x, y y z en coordenadas geocéntricas.

A continuación, usamos la función porjectExtent para proyectar los valores de un objeto raster a otro objeto raster nuevo que contiene otra proyección (CRS). Concretamente, proyectamos los valores del objeto cft al objeto projLL.

Finalmente, extraemos el cuadro delimitador espacial (spatial bounding box) de nuestro objeto espacial cftLL, lo convertimos en vector y lo asignamos a una nueva variable que toma el nombre de cftExt. Es decir, creamos una variable que recoge el alcance espacial del objeto cftLL en forma de vector.

projLL <- CRS('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0')
cftLL <- projectExtent(cft, projLL)
cftExt <- as.vector(bbox(cftLL))

A continuación, utilzamos este alcance espacial almacenado en la variable cftExt para extraer las líneas del mapa. Para ello, dentro de la libreria map utilizamos la base de datos ‘worldHires’ la cual contiene aproximadamente dos millones de puntos representando las costas del mundo así como los límites nacionales. En definitiva, en este paso se extraen las fronteras estatales de la penísnula ibérica y Baleares así como Portugal y parte de los países colindantes: Francia y el norte de África.

boundaries <- map('worldHires',
                  xlim = cftExt[c(1, 3)], ylim = cftExt[c(2, 4)],
                  plot = FALSE)

El siguiente paso consiste en convertir el resultado del paso anterior en un objeto SpatialLines para poder representarlo a la hora de hacer el GIF.

boundaries <- map2SpatialLines(boundaries, proj4string = projLL)

Posteriormente, proyectamos la misma proyección que contiene el objeto cft (projLCC2d) para que ambos estén en el mismo sistema de coordenadas de referencia.

Si mostramoslor resultados, se puede observar el mapa creado y en el que se reflejará la información.

## Project to the projection of the cft object
boundaries <- spTransform(boundaries, CRS(projLCC2d))
plot(boundaries)

Creamos los frames y el GIF

En primer lugar indicamos la paleta de colores que queremos usar y lo asignamos al objeto cloudTheme.

A continuación, creamos los frames. Para ello, primero utilizamos la función trellis.device que sirve como incialización de un dispositivo de visualización con los parámetros gráficos adecuados. Indicamos el tipo de objeto, en este caso imágenes en formato png, el nombre que recibirán las imágenes (R28Plot seguido del número de la imagen), el directorio donde se van a almacenar y las dimensiones.

Posteriormente, usamos la función leverplot para la creación del mapa y añadimos los bordes fronterizos creados anteriormente.

cloudTheme <- rasterTheme(region = brewer.pal(n = 9, 'Blues'))

tmp <- tempdir()
old <- setwd(tmp)
trellis.device(png, file = paste0(old, '/R28Plot%02d.png'),
               res = 300, width = 1500, height = 1500)
levelplot(cft, layout = c(1, 1), par.settings = cloudTheme) +
  layer(sp.lines(boundaries, lwd = 0.6))
dev.off()

## png 
##   2

El últmo paso consiste en recoger las 96 imágenes creadas y almacenarlas en un objeto (png_files) y con la función av_encode_video del paquete av creamos el vídeo a partir del objeto que contiene las imágenes, indicamos el nombre que tomará el resultado y la velocidad a la que pasarán las imágenes creando la secuencia.

Con la función browserURL del paquete utils cargamos el vido obtenido.

setwd("C:/Users/sanch/OneDrive/Escritorio/BIA/2021-2022/Datos_Espaciales_y_Espacio_Temporales/practica10/gif")

png_files <- sprintf("R28Plot%02d.png", 1:96)
av::av_encode_video(png_files, 'output28_05.mp4', framerate = 6)

## [1] "C:\\Users\\sanch\\OneDrive\\Escritorio\\BIA\\2021-2022\\Datos_Espaciales_y_Espacio_Temporales\\practica10\\gif\\output28_05.mp4"

utils::browseURL('output28_05.mp4')

Si proyectamos el vídeo obtenido, se aprecia lo siguiente:

El dia 27 se caracteriza por la ausencia de nubes en la región de interés. No es hasta las útlimas horas del día y hasta las 12 horas del día 28 cuando se observan ciertos niveles de nubosidad en la región del Océano Atlántico asi como algunas nubes en la zona del País Vasco y el Mar Cantábrico. A partir del medio día del día 28 aparecen algunas nubes en la zonas de Teruel, Cuenca, Salamanca y Segovia las cuáles, para el final del día ya han desaparecido. Al mismo tiempo se van formando nuevos intervalos nubosos en el Océano Atlántico. A partir del día 29 se observa la presencia de nubes en más puntos de la península procedentes de un frente noboso que evoluciona en el Océano Cantábrico.

Comparando el mapa de interpolación en el que se muestran los niveles de concentración de NO2 en el aire en la penisnula iberica para el día 28 de mayo de este año con el vídeo creado en esta tarea para ese mismo día, se podría decir que el aumento de la presencia del gas contaminante se debe a la ausencia de precipitaciones. No obstante, la afirmarmación de esta declaración requiere de un trabajo de investigación en el que se estudie con mayor profundidad la posible relación contaminación-metereología.