How to plot wind vectors on wind map

Question

I have NAM wind models that I create map from but I am struggling to add average wind direction vectors. How do I convert the U and V vectors so they give average direction for a sub sample of wind angles? Currently when I plot wind direction the arrows are all swirling. Do I need to some how average wind direcions over chunks of area.

[map][1]

library(here)
library(tidyverse)
library(raster)
library(sf)
library(lubridate)
library(ggquiver)
library(magick)
library(metR)
library(rCAT)
library(circular)
library(rNOMADS)

crop_ras <- raster(xmn = -121.404,
             xmx = -117.064,
             ymn = 32.36,
             ymx = 34.75)

date <- as_date(Sys.Date())

mod.list <- c("00", "03", "06", "09", "12", "15", "18", "21", "27", 
              "30", "33", "36", "39", "42", "45", "51", "54", "57", "60", 
              "63", "66", "69", "75", "78", "81", "84")

gribstring <- "wind_u_%s.grb"
wind.u <- vector("list", length(mod.list))

#U Vector
unlink(here("grib_u/*"))
for(i in seq_along(mod.list)) {
  url_ucomp <- paste0("https://nomads.ncep.noaa.gov/cgi-bin/filter_nam.pl?file=nam.t00z.awphys" , mod.list[[i]] , ".tm00.grib2&lev_10_m_above_ground=on&var_UGRD=on&subregion=&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%2Fnam.20220512")
  gribfile <- sprintf(gribstring, mod.list[i])
  download.file(url_ucomp, here("grib_u", gribfile))
  
}

wind_u_brick <- stack(here("grib_u", list.files(path = here("grib_u"), pattern = "*.grb"))) %>% 
  projectRaster(crs = 4326) %>% 
  crop(crop_ras)

wind_u_csv <- rasterToPoints(wind_u_brick) %>% 
  as.data.frame() %>% 
  pivot_longer(!c(x,y), names_to = "date", values_to = "u")


# V vector
unlink(here("grib_v/*"))
for(i in seq_along(mod.list)) {
  url_ucomp <- paste0("https://nomads.ncep.noaa.gov/cgi-bin/filter_nam.pl?file=nam.t00z.awphys" , mod.list[[i]] , ".tm00.grib2&lev_10_m_above_ground=on&var_VGRD=on&subregion=&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%2Fnam.20220512")
  gribfile <- sprintf(gribstring, mod.list[i])
  download.file(url_ucomp, here("grib_v", gribfile))
  
}

wind_v_brick <- stack(here("grib_v", list.files(path = here("grib_v"), pattern = "*.grb"))) %>% 
  projectRaster(crs = 4326) %>% 
  crop(crop_ras)
  

wind_v_csv <- rasterToPoints(wind_v_brick) %>% 
  as.data.frame() %>% 
  pivot_longer(!c(x,y), names_to = "date", values_to = "v")

wind_data <- left_join(wind_u_csv, wind_v_csv, by = c("x", "y")) %>% 
  dplyr::select(!c("date.y")) %>% 
  rename(date = date.x)
  
wind_data$date <- str_replace(wind_data$date, "wind_u_", "")
  

wind_ras <- wind_data %>%
  mutate(wind = ((sqrt((u^2) + (v^2)))),
         angle = atan2(u/wind, v/wind),
         angle = (angle * 180/pi),
         date = as.factor(date)) %>%
  dplyr::select(-c(u, v)) %>% 
  dplyr::select(x, everything()) %>% 
  relocate(wind, .before = date) %>% 
  mutate(wind = wind*1.94)

wind_angle <- wind_ras %>%
  dplyr::select(c("x", "y", "angle", "date", "wind")) %>%
  group_by(x, y, date) %>%
  summarize(angle = mean(angle),
         wind = mean(wind),
         x =x,
         y=y,
         date = date) %>%
  mutate(angle = angle + 180)

wind_ras <- wind_ras %>% 
  dplyr::select(!c(angle))

date <-  unique(wind_ras$date)

ras <- brick(xmn = -121.4329,
             xmx = -117.1239,
             ymn = 32.3719,
             ymx = 34.7869,
             nrows = 23,
             ncol = 31)

for(i in 1:length(levels(wind_ras$date))) {

  ras[[i]] <- wind_ras %>% 
    filter(date == levels(date)[i]) %>% 
    rasterFromXYZ() %>% 
    dropLayer(2)
}

names(ras) <-  date

wind_pal <- colorRampPalette(c( "#7635b8", "#3a57ba", "#4999b3","#5bc752", "#c9b84d", "#cc8247", "#ab3c4d", "#bf4182" ))

wind_limits <- c(0,50)

wind_values <- c(0, 10, 20, 30, 40, 50)

wind.list = lapply(sort(unique(wind_ras$date)), function(i) {
  ggplot(data = wind_ras[wind_ras$date==i,], aes(y=y, x=x)) +
  geom_raster(aes(fill = wind), interpolate = T) +
  facet_wrap(~date) +
  scale_fill_gradientn(colours = wind_pal(100), limits = wind_limits, breaks = wind_values) +
  theme_classic() +
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank(),
        plot.margin=unit(c(0,0,0,0), "null"),
        strip.background = element_blank(),
        strip.text.x = element_blank()) +
  coord_sf(xlim = c(-121.45, -117.06), ylim = c(32.36, 34.75), expand = F)
})


wind.arrows <- wind_angle %>%
  ungroup() %>%
    slice(which(row_number() %% 100 ==1))

test_wind <- wind.arrows %>%
  ungroup() %>%
  filter(date == date[[2]]) %>%
  dplyr::select(!c("date"))

a <- wind.list[[2]]


a +
  geom_spoke(data = test_wind, aes( x = x, y = y, angle = angle,
           radius = scales::rescale(wind, c(.02, .15))), arrow = arrow(length = unit(.2, 'cm')))


  [1]: https://i.stack.imgur.com/o60am.png
  [2]: https://i.stack.imgur.com/aD6d1.png

Answer 1

Here is an example of such a map, a think. This may not directly answer your question, but it shows how to provide some minimal example code. Your code is unnecessarily complex for the question at hand. Perhaps you can expand the code below a bit and use it in your question.

library(rasterVis)
library(raster)
proj <- CRS('+proj=longlat +datum=WGS84')
df <- expand.grid(x = seq(-2, 2, .01), y = seq(-2, 2, .01))
  
df$z <- with(df, (3*x^2 + y)*exp(-x^2-y^2))
r <- rasterFromXYZ(df, crs=proj)
vectorplot(r, par.settings=RdBuTheme())

But perhaps your question is really about combining U and V values? In that case, the code you would show would just create two vectors, one representing plausible U values, and one with plausible V value; and what you have tried to combine the values.