rvad o cómo calcular el viento
By Pao Corrales in meteorology Español
September 10, 2019
Cómo parte de mi tesis de licenciatura implementé un algorítmo para calcular el viento horizontal a partir de observaciones de radares. Lo escribí en python y era super lento, seguramente por mi culpa! La buena noticia es que este año pasé el código al lado del bien y ahora hay un paquete de R.
El objetivo del paquete rvad es aproximar las componentes del viento horizontal a partir del viento radial medido por un radar Doppler usando el método Velocity Azimith Displya o VAD de Browning and Wexler (1968).
Installation
Por ahora está diponible la versión en desarrollo en GitHub, podés instalarla así:
# install.packages("remotes")
remotes::install_github("paocorrales/rvad")
## Skipping install of 'rvad' from a github remote, the SHA1 (31d393d3) has not changed since last install.
## Use `force = TRUE` to force installationEjemplo
Usualmente los datos de radar vienen en formatos de los más diversos y con suerte son transformados a netCDF. Dejo esa tarea para la casa. En este caso es necesario tener los datos en formato tidy, o sea una obsevación de viento radial para cada azimut, ángulo de elevación y distancia al radar o rango. A modo de nuestra hay un set de datos disponibles con el paquete:
library(rvad)
str(radial_wind) # datos de ejemplo, vienen con el paquete!
## Classes 'data.table' and 'data.frame': 2076960 obs. of 4 variables:
## $ range : int 250 750 1250 1750 2250 2750 3250 3750 4250 4750 ...
## $ radial_wind: num NA NA -4.189 0.419 -4.556 ...
## $ azimuth : num 214 214 214 214 214 ...
## $ elevation : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## - attr(*, ".internal.selfref")=<externalptr>El campo de viento radial para un ángulo de elevación tiene esta pinta:
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.3
one_elevation <- subset(radial_wind, elevation == unique(elevation)[3] &
!is.na(radial_wind))
ggplot(one_elevation, aes(azimuth, range)) +
geom_point(aes(color = radial_wind, size = range^2)) +
scale_radius(range = c(0, 0.05), guide = "none") +
scale_color_gradient2(low = "blue", high = "red") +
scale_y_continuous(limits = c(0, 100000)) +
coord_polar()
## Warning: Removed 259 rows containing missing values (geom_point).
Colores rojos indican que el viento se aleja del radar y los azules indican que el viento se mueve hacia el radar. Entonces, a partir del gráfico se puede ver que el viento es aproximadamente del noreste. Esto es bastante común durante la madrugada de un día de verano en las llanuras de Argentina (los datos corresponden al radar INTA Paraná el 14 de enero de 2016 a las 06 UTC).
vad_fit() hace algunos controles de calidad y luego ajusta un modelo sinusoidal a los datos. Como resultado devuelve el viento zonal (u) y meridional (v) para cada ángulo de elevación y rango. El paquete también incluye una función que aproxima la propagación del haz del radar por lo que también se obtiene la altura de cada observación.
VAD <- with(radial_wind, vad_fit(radial_wind, azimuth, range, elevation))
str(VAD)
## Classes 'rvad_vad' and 'data.frame': 5760 obs. of 7 variables:
## $ height : num 2.19 6.58 11 15.45 19.93 ...
## $ u : num NA NA NA NA NA NA NA NA NA NA ...
## $ v : num NA NA NA NA NA NA NA NA NA NA ...
## $ range : int 250 750 1250 1750 2250 2750 3250 3750 4250 4750 ...
## $ elevation: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ r2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ rmse : num NA NA NA NA NA NA NA NA NA NA ...
## - attr(*, "rvad_raw")= logi TRUELa función tiene un método de graficado rápido para revisar el resultado:
plot(VAD)
Los parámetros optimos para el control de calidad (R^2 mínimo, máximo de valores faltantes, etc.) puede depender de cada caso. Por eso es importante revisar distintos valores y chequear el resultado.
La segunda función importante vad_regrid() resume el resultado anterior en un único perfil de viento definiendo una grilla regular a partir de una regresión local.
wind_profile <- vad_regrid(VAD, layer_width = 100)La función también tiene un método de gráficado, así que rápidamente este es el resultado:
plot(wind_profile)
Si bien alguien con buen ojo puede interpretar como será el perfil vertical de viento a partir del campo de viento radial, nada mejor como un paquete de R para hacer cálculos y darnos la respuesta.
Por ahora eso es todo, con suerte el doctorado me va a dejar tiempo para seguir incluyendo nuevas herramientas. Por ejemplo, probar el algoritmo para otras estrategias de escaneo!