---
title: "Anisotropic shortest paths across the Strait of Gibraltar"
author: "Javier Fernández-López and Klaus Schliep"
date: "`r format(Sys.Date(), '%B %d, %Y')`"
bibliography: rWind.bib
nocite: |
  @terra
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Anisotropic shortest paths across the Strait of Gibraltar}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 6
)
suppressPackageStartupMessages(library(rWind))
if (utils::packageVersion("rWind") < "1.2.0") {
  stop(
    "This vignette requires rWind >= 1.2.0. Install the current package ",
    "source before knitting it."
  )
}
suppressPackageStartupMessages(library(terra))
suppressPackageStartupMessages(library(gdistance))
```

## Overview

This vignette shows the core rWind workflow: obtain wind data, represent them
as modern spatial rasters, translate wind speed and direction into anisotropic
movement conductance, and calculate directional least-cost paths.

The example crosses the Strait of Gibraltar. Because wind-assisted movement is
direction dependent, the least-cost path from Spain to Morocco does not
necessarily match the path in the opposite direction.

## Obtain wind data

`wind.dl()` automatically selects the appropriate data source. Current dates
are requested from PacIOOS, whereas older dates are retrieved from the
NOAA/NCEI GFS archive. GFS data are available at three-hour UTC intervals.

The download call used for this example is:

```{r download, eval=FALSE}
w <- wind.dl(
  2015, 2, 12, 12,
  lon1 = -7, lon2 = -4,
  lat1 = 34.5, lat2 = 37.5
)
```

To make the vignette reproducible without a network connection, it uses the
result saved with the package.

```{r cached-data}
w <- readRDS("w.rds")
head(w)
```

## Create a `SpatRaster`

`wind2raster()` converts the tabular result into a two-layer
`terra::SpatRaster` [@terra]. The layers are named `direction` and `speed`, and
the coordinate reference system is WGS 84 (EPSG:4326).

```{r create-raster}
wind_layer <- wind2raster(w)
wind_layer
```

Individual layers use standard `terra` syntax:

```{r inspect-speed}
wind_layer[["speed"]]
```

## Build directional conductance

For passive movement, `flow.dispersion()` applies a restricted version of the
formulation of Felicísimo et al. [-@felicisimo2008]. Deviations of 90 degrees
or more from the flow are treated as barriers. The resulting `TransitionLayer`
stores directional conductance and can be used by `gdistance` [@gdistance].
Internally, rWind converts only the direction layer to the legacy raster
representation required by `gdistance`; the public spatial representation
remains a `SpatRaster`.

```{r conductance}
conductance <- flow.dispersion(wind_layer, type = "passive")
conductance
```

Every transition cost is the mean of the local costs at the origin and
destination cells. Diagonal transitions are multiplied by `sqrt(2)`. This
corrects diagonal length on an abstract square grid, but not the changing
physical distances of longitude/latitude cells.

### Active movement

The optional active model is deliberately general rather than mechanistic. It
uses the signed component of flow parallel to movement. Let
`z = speed * cos(angle) / speed.scale`: supportive flow (`z >= 0`) has cost
`1 / (1 + z)`, while opposing flow has cost `1 - z`. Perpendicular flow is
neutral and movement against the flow remains possible.

```{r active-conductance}
active_conductance <- flow.dispersion(
  wind_layer,
  type = "active",
  speed.scale = median(w$speed[w$speed > 0], na.rm = TRUE)
)
active_conductance
```

If `speed.scale` is omitted, rWind uses the median positive speed of the raster.
Supply one fixed value when comparing dates or regions. Active costs are
relative connectivity indices: they are not estimates of travel time, energy,
or an organism's intrinsic movement speed.

## Calculate paths in both directions

```{r paths}
spain <- c(-5.5, 37)
morocco <- c(-5.5, 35)

spain_to_morocco <- shortestPath(
  conductance, spain, morocco,
  output = "SpatialLines"
)
morocco_to_spain <- shortestPath(
  conductance, morocco, spain,
  output = "SpatialLines"
)
```

The following map uses `terra` for the raster and path representation. Wind
arrows are drawn from the downloaded U and V components; their lengths are
scaled for display.

```{r map}
plot(
  wind_layer[["speed"]],
  col = hcl.colors(20, "YlOrRd", rev = TRUE),
  main = "Wind-assisted paths across the Strait of Gibraltar",
  xlab = "Longitude", ylab = "Latitude"
)

lines(vect(spain_to_morocco), col = "#C62828", lwd = 3, lty = 2)
lines(vect(morocco_to_spain), col = "#1565C0", lwd = 3, lty = 2)
points(spain[1], spain[2], pch = 19, cex = 1.3, col = "#C62828")
points(morocco[1], morocco[2], pch = 19, cex = 1.3, col = "#1565C0")

arrow_scale <- 0.12 / max(w$speed, na.rm = TRUE)
arrows(
  w$lon, w$lat,
  w$lon + w$ugrd10m * arrow_scale,
  w$lat + w$vgrd10m * arrow_scale,
  length = 0.035, col = "grey20"
)

text(spain[1], spain[2] + 0.18, "Spain", col = "#C62828", font = 2)
text(morocco[1], morocco[2] - 0.18, "Morocco", col = "#1565C0", font = 2)
legend(
  "topleft",
  legend = c("Spain to Morocco", "Morocco to Spain"),
  col = c("#C62828", "#1565C0"),
  lwd = 3, lty = 2, cex = 0.8, bg = "white"
)
```

The two paths differ because the conductance matrix is asymmetric: travelling
with the wind has a different cost from travelling against it. For larger
areas or long time series, both raster size and the number of possible cell
transitions increase rapidly, so start with the smallest spatial extent and
resolution needed for the biological question.

## References