适配沿海两点海上路径代码计算海运距离及相关问题求助

问题1：距离输出不符合预期的原因 #

• 核心错误是经纬度参数顺序倒置：你定义起点、终点坐标时把纬度放在了前、经度放在了后，但project()函数要求输入顺序为c(经度, 纬度)，你传入的两个点实际并不在西班牙境内，因此计算结果和预期的1600km偏差极大。
• 次要影响因素是栅格分辨率过低：300×300的全球栅格单个单元格对应约110km的实际宽度，本身就会带来百km级的计算误差。
• 修复方案：
  1. 调整坐标顺序为c(经度, 纬度)，比如西班牙北部桑坦德港坐标为c(-3.81, 43.46)，南部阿尔赫西拉斯港坐标为c(-5.45, 36.14)
  2. 可根据研究区域提升栅格分辨率，或改用UTM等平面投影坐标系，进一步降低距离计算误差
  3. costDistance()返回值单位为米，除以1000即可得到公里单位的海运距离。

问题2：支持绕行小岛的实现方案 #

现有框架无需修改核心逻辑，仅替换陆地基础数据即可实现需求：

• 原有代码使用的wrld_simpl是低精度全球矢量数据集，未包含小型岛屿的轮廓，你可以改用GSHHG全球高精度海岸线数据集，它提供5个精度等级的矢量数据，最高可识别百米级小岛，直接替换原有wrld_simpl做栅格化即可
• 如果仅需计算特定区域的路径，可先裁剪对应区域的高精度陆地矢量，还能大幅提升计算效率。

修正后参考代码 #

library(raster)
library(gdistance)
library(maptools)
library(rgdal)
library(maps)

# 加载全球陆地矢量数据，可替换为更高精度的GSHHG数据
data(wrld_simpl)
world_crs <- crs(wrld_simpl)
world <- wrld_simpl
worldshp <- spTransform(world, world_crs)

# 可根据研究区域调整栅格范围和分辨率，此处以伊比利亚半岛区域为例提升精度
ras <- raster(nrow=1000, ncol=1000, xmin=-10, xmax=3, ymin=35, ymax=44)

# 栅格化处理，海洋设为1，陆地设为999（高通行成本）
worldmask <- rasterize(worldshp, ras)
worldras <- is.na(worldmask)
worldras[worldras==0] <- 999

# 构建过渡层并做地理校正
tr <- transition(worldras, function(x) 1/mean(x), 16)
tr = geoCorrection(tr, scl=FALSE)

# 西班牙北部桑坦德港（经度、纬度顺序）
port_origin <- structure(c(-3.81, 43.46), .Dim = 1:2)
port_origin <- project(port_origin, crs(world_crs, asText = TRUE))

# 西班牙南部阿尔赫西拉斯港（经度、纬度顺序）
port_destination <- structure(c(-5.45, 36.14), .Dim = 1:2)
port_destination <- project(port_destination, crs(world_crs, asText = TRUE))

# 计算路径和距离
distance <- costDistance(tr, port_origin, port_destination)
# 输出公里级距离
print(paste0("海运距离：", round(distance[1,1]/1000, 2), " km"))

path <- shortestPath(tr, port_origin, port_destination, output = "SpatialLines")

# 可视化
A <- accCost(tr, port_origin)
A <- mask(A, worldmask, inverse=TRUE)
plot(A)
points(rbind(port_origin, port_destination))
lines(path)

# 封装函数
RouteShip <- function(from_port, to_port, cost_raster, land_mask) {
  tr <- transition(cost_raster, function(x) 1/mean(x), 16)
  tr = geoCorrection(tr, scl=FALSE)
  A <- accCost(tr, from_port)
  A <- mask(A, land_mask, inverse=TRUE)
  path <- shortestPath(tr, from_port, to_port, output = "SpatialLines")
  plot(A)
  points(rbind(from_port, to_port))
  lines(path)
}

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