一种方法是

• 为速度向量的水平和垂直分量分别构建转换矩阵
• 构建水平和垂直位置的转换矩阵
• geoCorrect以上
• 从水平和垂直分量重建合成电导

示例：

• 我使用了一些map ping命令（带有{purrr}），但如果您不想保持较低的软件包数量，可以使用.apply变体或其他方法来实现 *

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

设置一些全局变量（仅用于演示）

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

创建岛屿掩码（问题中的代码）

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

创建“栅格”列表并开始填充：

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

根据当前方向创建具有水平和垂直速度分量的栅格：

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

创建经度和纬度位置的栅格（值从西向东和从北向南递增）。我们需要这些栅格将位置偏移与速度分量进一步相乘：

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

在位置处创建具有速度的光栅（这些是转换矩阵的基础）：

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

创建转换矩阵。转换是不对称的（取决于相对于速度分量的运动方向），因此symm = FALSE。添加一个轻微的偏移量（.001），以避免在90°、180°等Angular 处x或y速度分量的零矩阵（* 注：当使用参数x调用transitionFunction时，x[1]是“from”单元格的值，x[2]是“to”单元格的值。）*

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过渡层（因为对角相邻单元比正交单元更远离）：

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

根据经度和纬度分量计算电导率代理（待最大化）：

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

创建辅助函数以将过渡矩阵的值完全移动为正值（某些几何距离条件不允许负值）：

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
    }

移位转移矩阵，返回列表中最短路径和代价距离：

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

输出

在将代码 Package 到函数get_shortest_path中之后：

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

...并计算方位角30、60... 360°的最短路径和成本距离：

• 蓝色箭头表示电流的方向
• cost：从A移动到B的相对成本