Mapping wind data with R

August 28, 2022 | Milos Popovic

From mills to modern turbines, humanity has used wind energy for thousands of years. Wind is the fastest-growing renewable and is inexhaustible and reduces the use of fossil fuels. Visualizing wind direction and speed helps us understand the potential of this fascinating energy source.

Recently, I bumped into a lovely animation of the wind direction on windy.com 😍. Immediately, it made me think of the way to remake this map in R.

While I encountered many interesting tutorials here, here, and here, they all used satellite imagery to show wind movement in the form of heatmap rather than vectors. So, I decided to give it a try and make a different wind map.

On my voyage, I was lucky to find package rWind that brings wind data from Global Forecast System to your drive. This fascinating library allows you to donwload wind data for any place on the planet in just a few minutes. metR is another great package that helped me turn wind data into streamlines. Finally, I learned a great deal about mapping streamlines using metR from this excellent post on ocean current.

Standing on the shoulders of these important sources, let’s jump straight into the creation of a wind map of Europe!

Import libraries

We will make our wind map of Europe using several critical libraries: tidyverse and sf for spatial analysis and data wrangling; giscoR for importing the shapefile of Europe; rWind for downloading wind data in tidy format; classInt for creating bins; lubridate for defining timestamps for downloading wind data; metR for the creation of streamlines; and oce to interpolate the U/V components.

# libraries we need
libs <- c(
  "tidyverse", "sf", "giscoR",
  "lubridate", "classInt",
  "rWind", "metR", "oce"
)

# install missing libraries
installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == F)) {
  install.packages(libs[!installed_libs])
}

# load libraries
invisible(lapply(libs, library, character.only = T))

Get wind data

We first fetch hourly wind data using rWind. In this case, we download the yesterday’s the data but you can retrieve data for multiple days.

Since we’ll make a static wind map of Europe, we define the rough coordinates of Europe and calculate the average wind metrics.

# 1. GET WIND DATA
#----------------------

get_wind_data <- function(time_range, mean_wind_data, eur_wind_df) {

  time_range <- seq(ymd_hms(paste(2022, 8, 27, 00, 00, 00, sep = "-")),
    ymd_hms(paste(2022, 8, 28, 00, 00, 00, sep = "-")),
    by = "1 hours"
  )

  mean_wind_data <- rWind::wind.dl_2(time_range, -28.5, 58.5, 34.0, 73.5) %>%
    rWind::wind.mean()

  eur_wind_df <- as.data.frame(mean_wind_data)
  return(eur_wind_df)
}

Package rWind returns several important features for every available coordinate point. One of them is wind direction (dir) while another is wind speed in meters per second. Wind speed is computed from the U (ugrd10m) and the V (vgrd10m) vector components using the Pythagorean Theorem. The U and V vector components show where the wind comes from. The V component is parallel to the latitude with negative values for the wind that originates from the north and positive values for the wind that comes from the south. The U component is parallel to the longitude with negative values for the wind that originates from the east and positive values for the wind that comes from the west.

     time                lat   lon     ugrd10m   vgrd10m      dir    speed
4504 2022-06-24 23:00:00 73.5 -28.5  0.61497181 -1.664331 159.7207 1.774313
4505 2022-06-24 23:00:00 73.5 -28.0  0.50460491 -1.593729 162.4312 1.671705
4506 2022-06-24 23:00:00 73.5 -27.5  1.30171985 -2.266284 150.1276 2.613526
4507 2022-06-24 23:00:00 73.5 -27.0  1.13085511 -3.003650 159.3690 3.209478
4508 2022-06-24 23:00:00 73.5 -26.5  0.82255141 -1.608670 152.9183 1.806768
4509 2022-06-24 23:00:00 73.5 -26.0 -0.06343056 -1.179921 183.0772 1.181624

Transform wind data

We can calculate wind trajectory only after we have transformed wind data to coordinates.

get_wind_points <- function() {
  # convert wind data into points
  eur_wind_pts <- eur_wind_df %>%
    sf::st_as_sf(coords = c("lon", "lat")) %>%
    sf::st_set_crs(4326)
  
  return(eur_wind_pts)
}

eur_wind_pts <- get_wind_points()

Once we got the coordinates we will assign every wind point to a grid. The points are currently scattered across the map so we need an equal area to compute average values. In the chunk below, I choose a random grid size of 80x100 but you can toy with different sizes. Also note that I assigned a unique id to every grid.

get_wind_grid <- function() {
  eur_wind_grid <- eur_wind_pts %>%
    sf::st_make_grid(n = c(80, 100)) %>%
    sf::st_sf() %>%
    dplyr::mutate(id = row_number())
  
  return(eur_wind_grid)
}

eur_wind_grid <- get_wind_grid()

In the following step, we determine where each wind point falls on the grid and convert the output into a data frame. Then we compute the average wind speed as well as the U and V components for every grid unit. Finally, we merge the output with the grid object.

get_wind_grid_aggregated <- function() {

  eur_wind_grid_agg <- 
    sf::st_join(eur_wind_pts, eur_wind_grid, 
      join = sf::st_within) %>%
    sf::st_drop_geometry() %>%
    dplyr::group_by(id) %>%
    dplyr::summarise(
      n = n(), u = mean(ugrd10m), 
      v = mean(vgrd10m), speed = mean(speed)
    ) %>%
    dplyr::inner_join(eur_wind_grid, by="id") %>%
    dplyr::select(n, u, v, speed, geometry) %>%
    sf::st_as_sf() %>%
    na.omit()
    
  return(eur_wind_grid_agg)
}

eur_wind_grid_agg <- get_wind_grid_aggregated()

Interpolation

So far so good! Now that we have wind data points clustered by an equal size grid we are ready to interpolate the U and V components. Interpolation is the method for assigning the same wind speed and direction to the nearest coordinate point. Since we ponly have information about the wind speed and direction from several measuring stations, we use this approach to infer values for other coordinates. Once we infer these values, we will be able to dsiplay wind as a particle flow.

But, to do so, we first need to retrieve the centroid points for every grid unit.

get_wind_coords <- function() {
  coords <- eur_wind_grid_agg %>%
    st_centroid() %>%
    st_coordinates() %>%
    as_tibble() %>%
    rename(lon = X, lat = Y)
  return(coords)
}

coords <- get_wind_coords()

Then we merge the centroid points and the gridded wind data and convert the object to a data.frame to facilitate interpolation.

get_wind_df <- function() {
  eur_df <- coords %>%
    bind_cols(sf::st_drop_geometry(eur_wind_grid_agg))
  return(eur_df)
}

eur_df <- get_wind_df()

First we interpolate the U component using the Barnes algorithm implemented in the oce package. This algorithm, in short, calculates the wavelength values at every grid point given a distance weighting between the grid and wind points. These weights are then used to ascertain the importance of every measurement on the calculation of the function values. Finally, the algorithm uses the spectral response of the interpolated points to make corrections passes and optimize the function values.

Essentially, the Barnes algorithm takes three components: a vector of x and y grid locations and a vector of z values, one at each (x,y) location. So, below we pass grid centroids as x and y values and U component at every x,y location as our z values.

## interpolate the U component
get_u_interpolation <- function() {
    wu <- oce::interpBarnes(
        x = eur_df$lon,
        y = eur_df$lat,
        z = eur_df$u
    )
    return(wu)
}

wu <- get_u_interpolation()

Then we define the width and length for tranforming the wide table of values into a long table

dimension <- data.frame(lon = wu$xg, wu$zg) %>% dim()

This will help us make a U component data table from the interpolated matrix. The chunk below creates a data.frame interpolated U component values for every grid coordinate point.

## make a U component data table from interpolated matrix
get_u_table <- function() {
    udf <- data.frame(
        lon = wu$xg,
        wu$zg
    ) %>%
        gather(key = "lata", value = "u", 2:dimension[2]) %>%
        mutate(lat = rep(wu$yg, each = dimension[1])) %>%
        select(lon, lat, u) %>%
        as_tibble()

    return(udf)
}

udf <- get_u_table()

In the code below, we repeat the previous steps for the V component.

## interpolate the V component
get_u_interpolation <- function() {
    wv <- oce::interpBarnes(
        x = eur_df$lon,
        y = eur_df$lat,
        z = eur_df$v
    )
    return(wv)
}

wv <- get_u_interpolation()

## make the V component data table from interpolated matrix
get_v_table <- function() {
    vdf <- data.frame(lon = wv$xg, wv$zg) %>%
        gather(key = "lata", value = "v", 2:dimension[2]) %>%
        mutate(lat = rep(wv$yg, each = dimension[1])) %>%
        select(lon, lat, v) %>%
        as_tibble()
    return(vdf)
}

vdf <- get_v_table()

In the final step, we merge the V and U data table and compute the velocity.

## merge the V and U component tables and compute velocity
get_final_table <- function() {
    df <- udf %>%
        bind_cols(vdf %>% select(v)) %>%
        mutate(vel = sqrt(u^2 + v^2))
    return(df)
}

df <- get_final_table()

Here is the output that you should get if you followed the steps

head(df)
# A tibble: 6 x 5
    lon   lat       u     v   vel
  <dbl> <int>   <dbl> <dbl> <dbl>
1   -28    34  0.853  -5.00  5.08
2   -26    34  0.0866 -3.73  3.73
3   -24    34 -0.507  -5.84  5.86
4   -22    34 -1.92   -6.04  6.34
5   -20    34 -2.91   -5.41  6.14
6   -18    34 -1.68   -6.14  6.37

Mapping wind data

Hurray, we have the wind data in proper shape! The time is ripe to prepare the map of Europe and display the values. We start off by calling the Europe vector map and defining the bounding box for Europe.

# 5. MAP WIND DATA
#-------------------------

get_europe_sf <- function() {
    eur_sf <- giscoR::gisco_get_countries(
        year = "2016", epsg = "4326",
        resolution = "10", region = c("Europe", "Asia")
    )

    return(eur_sf)
}

# bounding box
crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"

get_bounding_box <- function(bbox, bb) {
    bbox <- st_sfc(
        st_polygon(list(cbind(
            c(-25, 48.5, 48.5, -25, -25),
            c(32.000, 32.000, 69.5, 69.5, 32.000)
        ))),
        crs = crsLONGLAT
    )

    bb <- sf::st_bbox(bbox)

    return(bb)
}

In the final step before mapping, we customize the color palette and define breaks for wind velocity.

# colors
cols <- c(
  '#feebe2', '#d84594', '#bc2b8a', '#7a0177'
)

newcol <- colorRampPalette(cols)
ncols <- 6
cols2 <- newcol(ncols)

# breaks
vmin <- min(df$vel, na.rm = T)
vmax <- max(df$vel, na.rm = T)

brk <- classInt::classIntervals(df$vel,
  n = 6,
  style = "fisher"
)$brks %>%
  head(-1) %>%
  tail(-1) %>%
  append(vmax)

breaks <- c(vmin, brk)

It’s time to make the map, people!

Below is quite a long chunk where we essentially use metR::geom_streamline to create wavelengths with the input values for every coordinate point. This function takes several values. Apart from x and y, there are dx and dy variables that allow us to pass the U and V components.

The remaining code plugs the national map of Europe and integrates the customized palette choice.

make_wind_map <- function(eur_sf, bb) {
    eur_sf <- get_europe_sf()
    bb <- get_bounding_box()

    p <- df %>%
        ggplot() +
        metR::geom_streamline(
            data = df,
            aes(
                x = lon, y = lat, dx = u, dy = v,
                color = sqrt(..dx..^2 + ..dy..^2)
            ),
            L = 2, res = 2, n = 60,
            arrow = NULL, lineend = "round",
            alpha = .85
        ) +
        geom_sf(
            data = eur_sf,
            fill = NA,
            color = "#07CFF7",
            size = .25,
            alpha = .99
        ) +
        coord_sf(
            crs = crsLONGLAT,
            xlim = c(bb["xmin"], bb["xmax"]),
            ylim = c(bb["ymin"], bb["ymax"])
        ) +
        scale_color_gradientn(
            name = "Average speed (m/s)",
            colours = cols2,
            breaks = breaks,
            labels = round(breaks, 1),
            limits = c(vmin, vmax)
        ) +
        guides(
            fill = "none",
            color = guide_legend(
                override.aes = list(size = 3, alpha = 1, shape = 15),
                direction = "horizontal",
                keyheight = unit(2.5, units = "mm"),
                keywidth = unit(15, units = "mm"),
                title.position = "top",
                title.hjust = .5,
                label.hjust = .5,
                nrow = 1,
                byrow = T,
                reverse = F,
                label.position = "bottom"
            )
        ) +
        theme_bw() +
        theme(
            text = element_text(family = "Montserrat"),
            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(),
            legend.position = c(.5, .05),
            legend.text = element_text(size = 60, color = "white"),
            legend.title = element_text(size = 80, color = "white"),
            legend.key = element_blank(),
            plot.margin = unit(c(t = 1, r = -2, b = -1, l = -2), "lines"),
            panel.grid.major = element_line(color = "#070C33", size = 0.2),
            panel.grid.minor = element_blank(),
            plot.title = element_text(
                face = "bold", size = 120,
                color = "white", hjust = .5
            ),
            plot.caption = element_text(
                size = 50, color = "grey80",
                hjust = .25, vjust = -10
            ),
            plot.subtitle = element_blank(),
            plot.background = element_rect(fill = "#070C33", color = NA),
            panel.background = element_rect(fill = "#070C33", color = NA),
            legend.background = element_rect(fill = "#070C33", color = NA),
            legend.spacing.y = unit(.5, "pt"),
            panel.border = element_blank()
        ) +
        labs(
            x = "",
            y = "",
            title = "Average wind speed in Europe (27 August 2022)",
            subtitle = "",
            caption = ""
        )
    return(p)
}

p <- make_wind_map()

ggsave(
    filename = "eur_wind_27august2022.png",
    width = 8.5, height = 7, dpi = 600, device = "png", p
)

Aaaaand, here is the final map! 😍

Alright, that’s all, my friend! In this tutorial, you learned how to download, reshape and visualize wind data in the form of streamlines and all that in R! Feel free to check the full code here, clone the repo and reproduce, reuse and modify the code as you see fit.

I’d be happy to hear your view on how this map could be improved or extended to other geographic realms. To do so, please follow me on Twitter, Instagram or Facebook! Also, feel free to support my work by buying me a coffee here!

wind

maps

ggplot2

dataviz