R语言gdistance计算避岸两点距离返回Inf的transition函数问题排查

问题背景 #

计算西法边境Banyuls岸段两点的避岸航行距离时，已完成依赖包加载、岸线数据下载、矢量预处理，生成了投影为EPSG:25831、分辨率100m的栅格对象r，预设规则为栅格值1代表不可穿越陆地，小于1代表可通行水域。调用gdistance包计算成本距离时返回结果为Inf，已确认水域视觉上连通。

预处理代码 #

library(sf)
library(fasterize)
library(raster)
library(dplyr)
library(tidyverse)
library(stars)
library(gdistance)

spainurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E3417%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(spainurl, "spain.zip", mode = "wb")
unzip("spain.zip", exdir = "spain", junkpaths = TRUE)

franceurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E19888%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(franceurl, "france.zip", mode = "wb")
unzip("france.zip", exdir = "france", junkpaths = TRUE)

spainCoast_CoteBanyuls <- list.files("spain",
                                     pattern = "shp$",
                                     full.names = TRUE) %>% st_read()

frenchCoast_CoteBanyuls <- list.files("france",
                                      pattern = "shp$",
                                      full.names = TRUE) %>% st_read()

lines_spain <- st_geometry(spainCoast_CoteBanyuls) %>% st_cast("LINESTRING")
spainCoast_l <- st_sf(n = as.character(seq_len(length(lines_spain))), lines_spain)

lines_france <- st_geometry(frenchCoast_CoteBanyuls) %>% st_cast("LINESTRING")
franceCoast_l <- st_sf(n = as.character(seq_len(length(lines_france))), lines_france)

spainmax <- spainCoast_l[which.max(st_length(spainCoast_l)), ]
spainrest <- spainCoast_l[-which.max(st_length(spainCoast_l)), ]

joined <- c(st_geometry(spainmax), st_geometry(franceCoast_l)) %>%
  st_union()

join_end <- st_union(joined, st_geometry(spainrest))

bbox_all <- st_bbox(joined) %>%
  st_as_sfc()

polygon_joined <- bbox_all %>%
  lwgeom::st_split(join_end) %>%
  st_collection_extract("POLYGON")

#Polygons on position 2 and 3 need to be removed (visual inspection)
polygon_end <- polygon_joined[-c(2:3)]
polyCombin_df <- st_sf(var = 1, polygon_end)
polyCombin_df_t <- polyCombin_df %>% st_transform(25831)

r <- raster(polyCombin_df_t, res = 100)
r <- fasterize(polyCombin_df_t, r, fun = "max")

par(mar=c(2,2,1,1))
plot(r)

报错的计算代码 #

# 待计算两点：A = c(505470.9,4709364)，B = c(513902.9,4697052)
tr <- transition(r,transitionFunction= function(x) 1/mean(x), directions= 8)
costDistance(tr, c(505470.9,4709364), c(513902.9,4697052))

根因定位 #

• 核心问题是栅格值不符合预期：fasterize仅会对陆地面要素覆盖的栅格赋值为1，未被陆地覆盖的水域栅格默认填充为NA，并非预设的“小于1的可通行值”。gdistance计算连通性时会将NA自动识别为绝对障碍，直接阻断所有通行路径，因此返回无穷大结果。
• 缺少转移矩阵地理校正步骤：8邻域生成的原始转移矩阵没有校正斜向、横向/纵向移动的实际几何距离差，就算排除NA问题，计算出的距离也存在系统偏差。
• 坐标传入不规范：直接给costDistance传入数值向量时不会校验坐标参考系，若坐标偏移落到陆地栅格上，也会触发Inf结果。

修复方案 #

1. 修正栅格通行成本赋值：将水域的NA替换为通行成本1，将陆地赋值为极高成本（等效不可穿越，避免使用NA）
2. 生成转移矩阵后做地理校正，修正不同方向移动的距离偏差
3. 传入坐标时绑定对应坐标系，提前校验两点落在可通行水域

修复后可运行代码 #

# 1. 修正栅格成本值
r[is.na(r)] <- 1       # 水域赋值基础通行成本1
r[r == 1] <- 10^6      # 陆地赋值极高通行成本，等效不可穿越

# 2. 构建并校正转移矩阵
tr <- transition(r, transitionFunction = function(x) 1/mean(x), directions = 8)
tr <- geoCorrection(tr, type = "c", multpl = FALSE)

# 3. 校验两点是否落在可通行区域
pts <- matrix(c(505470.9,4709364,513902.9,4697052), ncol = 2, byrow = TRUE)
point_cost <- extract(r, pts)
if(any(point_cost > 1000)) stop("计算点落在陆地上，请检查坐标")

# 4. 计算避岸距离
dist <- costDistance(tr, pts[1,], pts[2,])
print(dist)

若需要严格的避岸要求，可以在岸线外围设置梯度递增的成本带，距离岸线越近通行成本越高，即可得到离岸有安全距离的最优路径。

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