如何结合海流与岛屿栅格生成gdistance过渡对象用于最短路径分析

我的最终目标是使用gdistance::shortestPath()计算岛屿两侧点位间的海上最短路径（不穿越陆地），同时考虑海流速度与流向。生成过渡层是实现该目标的中间步骤，但不清楚如何将海流速度栅格、海流方向栅格、含岛屿的海洋栅格传入gdistance::transition()，或是否需要先合并这些栅格。

准备代码 #

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)

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

卡住的步骤 #

目前卡在生成过渡对象的环节，不知道如何指定transitionFunction：

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

后续计划用生成的过渡对象计算路径：

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

解决方案 #

不需要合并栅格，而是要自定义过渡函数，同时调用三个栅格的数据来计算移动成本（gdistance的"最短路径"实际是最小成本路径，海流会影响移动的时间/成本）。

步骤说明 #

1. 将三个栅格打包为RasterStack，方便在过渡函数中同时调用；
2. 自定义过渡函数，逻辑如下：
   • 若相邻单元中存在陆地（island_mask值为999），则设置成本为Inf（不可通行）；
   • 计算当前单元到相邻单元的方位角；
   • 结合海流的速度与方向，计算实际移动速度（顺流时速度叠加，逆流时速度抵消）；
   • 用单元间的地理距离除以实际速度，得到移动成本；
3. 调用transition()时传入自定义函数，设置symm=FALSE（海流是单向的，往返成本不对称）。

完整代码 #

# 1. 将三个栅格合并为RasterStack
raster_stack <- stack(island_mask, ocean_spd, ocean_dir)

# 2. 自定义过渡函数：计算相邻单元的移动成本
current_transition <- function(x) {
  # x是2行3列的矩阵：每行对应一个单元，列分别为island_mask, ocean_spd, ocean_dir
  mask_from <- x[1,1]
  mask_to <- x[2,1]
  
  # 任意单元为陆地则设为不可通行
  if (mask_from == 999 || mask_to == 999) {
    return(Inf)
  }
  
  # 获取当前单元的海流数据
  spd <- x[1,2]
  dir <- x[1,3]
  
  # 获取栅格分辨率，用于计算单元距离和移动方向
  res <- res(raster_stack)
  dx <- res[1]
  dy <- res[2]
  
  # 获取两个单元的坐标，计算移动方位角
  from_cell <- which(raster_stack[[1]][] == mask_from)[1]
  to_cell <- which(raster_stack[[1]][] == mask_to)[1]
  coords_from <- xyFromCell(raster_stack, from_cell)
  coords_to <- xyFromCell(raster_stack, to_cell)
  
  # 计算移动方向（0-360度方位角）
  move_dir <- (atan2(coords_to[2] - coords_from[2], coords_to[1] - coords_from[1]) * 180 / pi) %% 360
  
  # 计算海流方向与移动方向的夹角（取最小角度）
  angle_diff <- abs(dir - move_dir)
  angle_diff <- ifelse(angle_diff > 180, 360 - angle_diff, angle_diff)
  angle_diff_rad <- angle_diff * pi / 180
  
  # 假设船只自身速度为1m/s，计算有效移动速度
  vessel_spd <- 1
  effective_spd <- vessel_spd + spd * cos(angle_diff_rad)
  # 避免速度为负导致成本异常
  effective_spd <- max(effective_spd, 0.01)
  
  # 计算单元中心的距离（16邻域含对角单元）
  cell_dist <- sqrt((coords_to[1]-coords_from[1])^2 + (coords_to[2]-coords_from[2])^2)
  
  # 返回移动成本（时间=距离/有效速度）
  return(cell_dist / effective_spd)
}

# 3. 生成过渡对象
trans_obj <- transition(raster_stack, transitionFunction = current_transition, directions = 16, symm = FALSE)

# 4. 转换点位格式并计算最小成本路径
pt_from <- SpatialPoints(poly_df[1,,drop=FALSE], proj4string = CRS(proj4string(raster_stack)))
pt_to <- SpatialPoints(poly_df[2,,drop=FALSE], proj4string = CRS(proj4string(raster_stack)))

pt_dist <- gdistance::shortestPath(trans_obj, pt_from, pt_to, output = "SpatialLines")

# 转换为sf对象并计算路径长度
sf::st_as_sf(pt_dist) %>%
  sf::st_length()

关键说明 #

• 自定义函数中的船只自身速度需根据实际场景调整；
• 若需计算最短距离而非最小时间，可将成本函数改为直接返回单元间地理距离，但会忽略海流影响；
• 过渡函数返回的成本值会影响路径选择，顺流区域成本更低，路径会优先选择顺流路线。

内容的提问来源于stack exchange，提问作者flee