Getting started with windscape

Wind fields

The basic unit of wind data is a wind field: the wind across a grid at a single moment, stored as two raster layers giving the eastward (u) and northward (v) components of the wind vector in each cell. wind_field() creates one from a two-layer SpatRaster, or from one time step of a wind time series. Here’s the first time step of the example data, for midnight UTC on January 1, 2000:

series <- windscape_example("wind_series")
field <- wind_field(series, step = 1)

ggplot(field, aes(x, y)) +
      geom_raster(aes(fill = speed)) +
      geom_path(data = states, aes(long, lat, group = group), color = "white", linewidth = 0.2) +
      geom_wind_arrow(color = "white") +
      scale_fill_viridis_c(name = "wind speed\n(m/s)") +
      coord_quickmap(xlim = c(-120, -90), ylim = c(30, 50), expand = FALSE) +
      theme_void()

geom_wind_arrow() draws wind vectors, and geom_wind_trail() draws trails that follow the airflow. Wind fields are useful for visualizing the wind at a given moment, but windscape’s connectivity models are built from many of them.

Wind time series

A wind_series is a sequence of wind fields stored as a single SpatRaster: all the u layers, followed by all the v layers, in the same time order. The example series covers the western and central US at about 0.3 degree resolution, with 96 time steps from the year 2000 (every 6 hours on the 1st and 15th of each month):

series
#> class       : SpatRaster
#> size        : 64, 96, 192  (nrow, ncol, nlyr)
#> resolution  : 0.3157895, 0.3174603  (x, y)
#> extent      : -120.1579, -89.84211, 29.84127, 50.15873  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> sources     : wind_usa.tif
#> names       : u 2000-01-01, u 200~00:00, u 200~00:00, u 200~00:00, u 2000-01-15, u 200~00:00, ...
#> min values  :         -0.4,        -0.4,        -0.6,        -0.6,         -0.4,        -0.5, ...
#> max values  :          0.7,           1,         0.9,           1,          0.7,         1.1, ...

Layer names record each time step. A real analysis would usually use far more data: windscape’s models describe the average behavior of the wind over the time series, and hourly data from many months or years captures its variability much better. Downsampled or averaged winds (e.g. monthly means) give poor results, because averaging cancels out winds blowing in opposite directions. How long a record to use depends on the question: decades for processes like gene flow that integrate over many generations, or a single season for dispersal during a flowering or sporulation period.

download_wind_data() downloads hourly wind data from NCAR’s Geoscience Data Exchange, with no account needed: ERA5 (1940 to present), CFSR (1979-2010), and CFSv2 (2011 to present). Data are clipped to a bounding box on the server and saved as one cached file per month, which wind_series() combines into one series:

files <- download_wind_data("era5", xlim = c(-120, -90), ylim = c(30, 50),
                       years = 2011:2020, hours = seq(0, 21, 3), dir = "~/wind_data")
series <- wind_series(files)

Wind data from other sources can be used too. Load it as a SpatRaster on a longitude/latitude grid, with wind components in m/s, and convert it with wind_series(). Its order argument describes how the layers are arranged: "uuvv" (all u layers, then all v layers) or "uvuv" (alternating).

series <- wind_series(terra::rast("my_wind_data.tif"), order = "uvuv")

Least-cost paths

A least-cost model treats the wind rose as a network of grid cells, where the cost of moving from a cell to a neighbor is the inverse of the conductance between them. The least-cost route between two places is the one with the lowest total cost, and its cost is a travel time, in hours if the rose was built with trans = 1 from wind speeds in m/s.

These travel times are best understood as a measure of accessibility by wind. Each step’s cost is the time to cross it at the long-run average wind speed in that direction, which includes zeros for all the time the wind spends blowing in other directions. A travel time is therefore the time to follow the best route, including time waiting at a standstill for the wind direction to align roughly with the direction of travel. Travel times are comparable across places and directions, which is what most analyses need, but they don’t predict when real particles arrive.

wind_graph() builds the network, in one of two directions. A "downwind" graph measures travel from the site to other places, and an "upwind" graph measures travel from other places to the site. least_cost() maps travel times between a site and every grid cell across the landscape:

down <- least_cost(rose, site, "downwind")
up <- least_cost(rose, site, "upwind")

d <- rbind(data.frame(as.data.frame(down, xy = TRUE), direction = "downwind: travel time from the site"),
           data.frame(as.data.frame(up, xy = TRUE), direction = "upwind: travel time to the site"))

ggplot(d, aes(x, y)) +
      geom_raster(aes(fill = pmax(hours, 10))) + # floor at 10 hours for the log scale
      geom_path(data = states, aes(long, lat, group = group), color = "white", linewidth = 0.15) +
      geom_contour(aes(z = hours), color = "black", breaks = c(100, 200, 500, 1000), linewidth = 0.15) +
      annotate("point", site[1], site[2], color = "black", size = 1.5) +
      facet_wrap(~direction) +
      scale_fill_gradientn(name = "hours", trans = "log10", values = c(0, .5, .7, .85, 1),
                           colors = c("cyan", "dodgerblue", "purple", "red", "orange")) +
      coord_quickmap(xlim = c(-120, -90), ylim = c(30, 50), expand = FALSE) +
      theme_void() +
      theme(strip.text = element_text(margin = margin(4, 0, 4, 0)))

Places to the east are reached quickly from the site, while places to the west are reached slowly, and the pattern reverses for travel to the site. least_cost_paths() returns the routes themselves, which can be drawn with geom_wind_trail().

Random walks

A random walk model instead simulates particles moving through the wind rose. At each time step, particles move from each cell to its neighbors at rates set by the conductances, so they spread along all routes in proportion to how much wind flows along them, rather than only along the fastest one. Particles are removed from the air at a constant rate, set by a half-life in hours, and the particles removed from each cell make up its deposition.

random_walk() has two modes. In "pulse" mode, particles are released once and tracked over time, which shows how a release spreads and drifts. Here are the airborne particles from a single release after one, three, and seven days, with a half-life of three days, each shown relative to its peak density (the total airborne mass declines as particles are deposited):

pulse <- random_walk(rose, site, mode = "pulse", half_life = 72,
                     iter = 168, record = c(24, 72, 168))

d <- as.data.frame(pulse$airborne, xy = TRUE)
d <- data.frame(d[c("x", "y")], hours = rep(c(24, 72, 168), each = nrow(d)),
                density = unlist(d[-(1:2)]))
d$hours <- factor(paste(d$hours, "hours"), levels = paste(c(24, 72, 168), "hours"))
d$density <- d$density / ave(d$density, d$hours, FUN = max) # relative to each panel's peak
d$density <- pmax(d$density, 1e-3) # cells the particles haven't reached are exactly zero

ggplot(d, aes(x, y)) +
      geom_raster(aes(fill = density)) +
      geom_path(data = states, aes(long, lat, group = group), color = "white", linewidth = 0.15) +
      annotate("point", site[1], site[2], color = "red", size = 1) +
      facet_wrap(~hours) +
      scale_fill_viridis_c(name = "relative\ndensity", trans = "log10",
                           limits = c(1e-3, NA), oob = scales::squish) +
      coord_quickmap(xlim = c(-120, -90), ylim = c(30, 50), expand = FALSE) +
      theme_void() +
      theme(strip.text = element_text(margin = margin(4, 0, 4, 0)))

In "stream" mode, particles are released continuously, and the model returns the long-run steady state: the share of particles airborne over each cell (residence), and the share deposited in each cell (deposition). This suits processes that integrate over many releases, such as seed rain or gene flow. Like least-cost models, random walks can run in either direction. A downwind walk shows where particles released at the site end up, while an upwind walk shows where particles deposited at the site came from (its origin layer):

down <- random_walk(rose, site, mode = "stream", direction = "downwind", half_life = 48)
up <- random_walk(rose, site, mode = "stream", direction = "upwind", half_life = 48)

The half-life controls how far particles travel: short half-lives keep them close to the source, while long half-lives let them spread across the landscape, as they would for long-lived propagules or for gene flow over many generations. With trans = 1 and winds in m/s, it is in hours.

Testing hypotheses

Pairwise wind connectivity can be compared with pairwise ecological data, such as genetic differentiation or gene flow, to test whether wind shapes ecological patterns. windscape provides tools for testing several kinds of hypotheses (see Kling and Ackerly 2021 for examples):

• Flow: is directional wind connectivity related to directional ecological flow, such as gene flow? Compare the matrices directly.
• Isolation: are sites with weaker wind connectivity in both directions more differentiated, for example genetically? pairwise_means() converts an asymmetric matrix into a symmetric one by averaging the two directions.
• Asymmetry: are directional imbalances in wind connectivity related to imbalances in ecological flow? pairwise_ratios() converts a matrix into log ratios of the two directions.

mantel_test() tests these relationships with Mantel tests, which assess significance by permutation since the values in a pairwise matrix aren’t independent. Unlike most Mantel implementations, it handles asymmetric matrices and multiple control variables. Here we run each test against simulated random data, using geographic distance as a control variable for the flow and isolation tests:

set.seed(1)
distance <- point_distance(sites)
gene_flow <- matrix(runif(100), 10) # simulated gene flow
gene_diff <- pairwise_means(matrix(runif(100), 10)) # simulated genetic differentiation

flow <- mantel_test(deposition, gene_flow, z = list(distance))
isolation <- mantel_test(pairwise_means(hours), gene_diff, z = list(distance))
asymmetry <- mantel_test(pairwise_ratios(hours), pairwise_ratios(gene_flow))

sapply(list(flow = flow, isolation = isolation, asymmetry = asymmetry),
       function(x) round(c(stat = x$stat, p.value = x$p.value), 3))
#>           flow isolation asymmetry
#> stat    -0.103     0.045     0.154
#> p.value  0.358     0.757     0.280

As expected for random data, none of the relationships here are significant. With only ten sites, these tests also have little power; real analyses often need more sites.

The pairwise connectivity article compares the two models’ matrices, covers site placement in more detail, and works through these tests with the birch landscape genetic data set.