## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE,
                      fig.width = 7, fig.height = 4.5, out.width = "100%", dpi = 96)
set.seed(1)

## ----workflow, echo = FALSE, fig.alt = "Diagram of the windscape workflow: hourly wind data downloaded from NCAR (or your own files) become a wind_series, which yields wind fields and a wind rose; pre-built wind roses can instead be downloaded directly. The wind rose feeds least-cost and random walk models, which produce path data, windshed maps, and pairwise matrices for hypothesis tests."----
knitr::include_graphics("workflow.svg")

## ----setup--------------------------------------------------------------------
library(windscape)
library(ggplot2)

states <- map_data("state")

## ----field--------------------------------------------------------------------
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()

## ----series-------------------------------------------------------------------
series

## ----download, eval = FALSE---------------------------------------------------
# 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)

## ----own-data, eval = FALSE---------------------------------------------------
# series <- wind_series(terra::rast("my_wind_data.tif"), order = "uvuv")

## ----rose---------------------------------------------------------------------
rose <- wind_rose(series, trans = 1)
rose

## ----rose-map, fig.height = 5-------------------------------------------------
rose <- windscape_example("wind_rose")

ggplot(rose, aes(x, y)) +
      geom_path(data = states, aes(long, lat, group = group), color = "gray70", linewidth = 0.2) +
      geom_wind_rose(res = 12) +
      coord_quickmap(xlim = c(-120.5, -89.5), ylim = c(29.5, 51), expand = FALSE) +
      theme_void()

## ----net-flow-----------------------------------------------------------------
ggplot(net_flow(rose), aes(x, y)) +
      geom_raster(aes(fill = speed)) +
      geom_path(data = states, aes(long, lat, group = group), color = "white", linewidth = 0.2) +
      geom_wind_trail(color = "white", hours = 24) +
      scale_fill_viridis_c(name = "net flow\n(km/h)") +
      coord_quickmap(xlim = c(-120, -90), ylim = c(30, 50), expand = FALSE) +
      theme_void()

## ----site---------------------------------------------------------------------
site <- cbind(-105, 40)

## ----least-cost, fig.height = 3.4---------------------------------------------
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)))

## ----pulse, fig.height = 2.8--------------------------------------------------
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)))

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

## ----windsheds, fig.height = 3.4----------------------------------------------
d <- rbind(data.frame(fortify(down)[c("x", "y")], value = fortify(down)$deposition,
                      windshed = "downwind: where particles released here land"),
           data.frame(fortify(up)[c("x", "y")], value = fortify(up)$origin,
                      windshed = "upwind: where particles landing here come from"))
d$value <- d$value / ave(d$value, d$windshed, FUN = max) # relative to each windshed's peak

ggplot(d, aes(x, y)) +
      geom_raster(aes(fill = value)) +
      geom_path(data = states, aes(long, lat, group = group), color = "white", linewidth = 0.15) +
      annotate("point", site[1], site[2], color = "red", size = 1.5) +
      facet_wrap(~windshed) +
      scale_fill_viridis_c(name = "relative\ndensity", trans = "log10",
                           limits = c(1e-4, 1), 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)))

## ----sites--------------------------------------------------------------------
sites <- cbind(lon = runif(10, -115, -95), lat = runif(10, 33, 47))

## ----pairwise-lc--------------------------------------------------------------
hours <- pairwise_least_cost(wind_graph(rose), sites)
round(hours[1:5, 1:5])

## ----check, message = TRUE----------------------------------------------------
check_cell_distance(rose, sites)

## ----pairwise-rw--------------------------------------------------------------
deposition <- pairwise_random_walk(rose, sites, half_life = 48)
signif(deposition[1:5, 1:5], 2)

## ----mantel-------------------------------------------------------------------
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))