Create transition object based on ocean current and direction for use in shortest path analysis

Question

My ultimate goal is to use gdistance::shortestPath() to calculate the shortest path distance from a point on one side of an island to a point on the other side of the same island, while not travelling over land (i.e. travel by sea), while accounting for ocean current speed and current direction.

An intermediate of the above is to generate a transition layer, which is subsequently passed to shortestPath(). I am having trouble understanding how to pass multiple rasters to gdistance::transition().

I have the following raster layers:

• ocean current speed
• ocean current bearing
• map of ocean with an island

A small example:

library(dplyr)
library(sf)
library(terra)
library(gdistance)
library(raster)

# Ocean current speed raster
ocean_spd <- terra::rast(nrow = 100,
                         ncol = 100, 
                         xmin = 0,
                         xmax = 1,
                         ymin = 0,
                         ymax = 1,
                         crs = NA) 

# set all ocean values to 0.5 (m/s)
ocean_spd <- terra::setValues(ocean_spd, 0.5) 

# Ocean current bearing raster
# set all current bearing to 90 (left to right)
ocean_dir <- terra::setValues(ocean_spd, 90) 


# Create island raster
# island corners
poly_df <- data.frame(x = c(0.5, 0.5, 0.8, 0.8),
                      y = c(0.5,0.8, 0.8, 0.5))

# island as a polygon
poly_sf <- poly_df %>% 
  sf::st_as_sf(coords = c("x","y")) %>%
  summarise(geometry = sf::st_combine(geometry)) %>%
  sf::st_cast("POLYGON") %>%
  sf::st_make_valid()

# convert polygon to raster
island_mask <- terra::rasterize(poly_sf, ocean_spd)
plot(island_mask)

# island cells = 999 and ocean cell = 1
island_mask[island_mask == 1] <- 999
island_mask[is.nan(island_mask)] <- 1

# remove island from current and bearing rasters
ocean_spd[island_mask == 999] <- NA
ocean_dir[island_mask == 999] <- NA

# convert rast objects to raster - needed by transition()
island_mask <- raster(island_mask)
ocean_spd <- raster(ocean_spd)
ocean_dir <- raster(ocean_dir)

# have a look at the layers
plot(island_mask)
text(island_mask, digits=1)

plot(ocean_spd)
text(ocean_spd, digits=1)

plot(ocean_dir)
text(ocean_dir, digits=1)

Calculate transition layer - this is where I am stuck, it is not clear to me how I integrate the island raster, the ocean speed raster and the ocean current bearing raster to get the transition object.

trans_obj <- transition(island_mask, transitionFunction = ?, 16, symm = FALSE)

From here I would then use the transition object (trans_obj) to calculate shortest path between two points on the edge of the island (e.g. ports) - via the ocean

pt_dist <- gdistance::shortestPath(trans_obj,
                                   poly_df[1,],
                                   poly_df[2,],
                                   output = "SpatialLines")

# get the distance
sf::st_as_sf(pt_dist) %>%
  sf::st_length()

My question:

how can I pass the island raster, the ocean speed raster and the ocean current bearing raster to transition() to generate the transition object? Or is there some other middle step that I need to do to combine my three rasters before passing them to transition()?

Any advice would be much appreciated.

Answer 1

One approach is to

• build a transition matrix for the speed vector's horizontal and vertical components each
• build transition matrices for horizontal and vertical position
• geoCorrect above
• reconstruct resultant conductance from horizontal and vertical components

Example:

I used some mapping (with {purrr}), but that could be achieved with .apply variants or other if you don't want to keep the number of packages low

    library(dplyr)
    library(sf)
    library(terra)
    library(gdistance)
    library(raster)
    library(purrr) ## for convenient list mapping

set some global variables (for the sake of demonstration only)

const_ocean_dir = 60 ## remove when passing this as a function argument
const_ocean_spd = .5 ## remove when passing this as a function argument
raster_template <- rast(nrow = 100, ncol = 100, xmin = 0, xmax = 1, ymin = 0, ymax = 1, crs = NA)
nr = nrow(raster_template); nc = ncol(raster_template)
dir_count = 16 ## number of transition directions

create island mask (code from your question)

    poly_df <- data.frame(x = c(0.5, 0.5, 0.8, 0.8),
                      y = c(0.5,0.8, 0.8, 0.5))
    poly_sf <- poly_df %>% 
      sf::st_as_sf(coords = c("x","y")) %>%
      summarise(geometry = sf::st_combine(geometry)) %>%
      sf::st_cast("POLYGON") %>%
      sf::st_make_valid()
    island_mask <- terra::rasterize(poly_sf, raster_template)
    island_mask[island_mask == 1] <- 999
    island_mask[is.nan(island_mask)] <- 1

create a list "rasters" and start populating it:

    rasters <- list()
    rasters$ocean_spd <- raster_template |> setValues(const_ocean_spd)
    rasters$ocean_dir <- raster_template |> setValues(const_ocean_dir)

create rasters with horizontal and vertical speed components, depending on current direction:

    ## a helper to convert bearing (degrees clockwise from North)
    ## to radiant (clockwise from East):
    bearing_to_rad <- function(bearing) (bearing - 90) * pi/180
    
    ## raster ocean speed into longitudinal and latitudinal component:
    rasters$ocean_spd_x <- with(rasters, ocean_spd * cos(bearing_to_rad(ocean_dir)))
    rasters$ocean_spd_y <- with(rasters, ocean_spd * sin(bearing_to_rad(ocean_dir)))

Create rasters of longitudinal and latitudinal position (values increasing from West to East and North to South). We need these to multiply position shift with speed components further down:

    rasters$pos_x <- raster_template |> setValues(matrix(1:nr, nr, nc, byrow = TRUE))
    rasters$pos_y <- raster_template |> setValues(matrix(1:nr, nr, nc))

Create rasters with speed at position (these are the basis for the transition matrices):

    rasters$spd_at_pos_x <- rasters$ocean_spd_x * rasters$pos_x
    rasters$spd_at_pos_y <- rasters$ocean_spd_y * rasters$pos_y

mask all rasters with island mask:

    ## mask rasters with island raster:
    rasters <- rasters |> map(~ {.x[island_mask == 999] = NA; .x})  

Create transition matrix. Transition is asymmetric (depends on direction of movement relative to speed components), thus symm = FALSE. Add a slight offset (.001) to avoid zero matrix for x- or y-speed component at angles of 90°, 180° etc. (Note: when transitionFunction is called with the argument x, x[1] is the "from" cell's value and x[2]that of the "to" cell.)

    tr_layers <- 
        rasters[c('spd_at_pos_x', 'spd_at_pos_y')] |>
        map(~ transition(.x |> raster(),
                         transitionFunction = \(from_to) .001 + from_to[2] - from_to[1],
                         directions = dir_count,
                         symm = FALSE
                         )
            )

geoCorrect transition layers (as diagonal neighbour cells are farther off than orthogonal ones):

    tr_layers <- tr_layers |> map(~ .x |> geoCorrection(type = 'c'))

calculate a conductance proxy (to be maximised) from longi- and latitudinal components:

    tr_layers$ocean_spd_resultant <-  with(tr_layers, spd_at_pos_x + spd_at_pos_y)

create a helper function to shift values of a transition matrix to entirely positive (negative values not allowed for some gdistance conditions):

     shift_to_positive <- function(tr_matrix){
         adj <- adjacencyFromTransition(tr_matrix)
         tmp <- tr_matrix[adj]
         tmp[is.na(tmp)] = 0
         tmp <- tmp - min(tmp, na.rm = TRUE)
         tr_matrix[adj] <- tmp
         tr_matrix[!adj] <- 0
         tr_matrix
    }

shift transition matrix and return shortest path and cost distance in a list:

    tr_layers$ocean_spd_resultant <- tr_layers$ocean_spd_resultant |> shift_to_positive()

    list(shortest_path = shortestPath(tr_layers$ocean_spd_resultant, A, B, output="SpatialLines"),
         cost_distance = costDistance(tr_layers$ocean_spd_resultant, A, B)
         )

output

after wrapping the code into a function get_shortest_path:

    get_shortest_path <- function(const_ocean_spd, const_ocean_dir, A = c(0, 0), B = c(1, 1){
      ## above code here
    }

... and calculating shortest paths and cost distances for bearings 30, 60 ... 360°:

• blue arrows indicate direction of the current
• cost: relative cost of moving from A to B