计算瑞士分形维度的R代码报错：Ruler is longer than maximum distance found

问题说明 #

运行R代码计算瑞士海岸线分形维度时，出现报错：

Error in measure_with_ruler(b, rulers[i] * 1000) : Ruler is longer than maximum distance found

但相同逻辑的哈萨克斯坦计算代码可正常运行，期望实现多尺子测量海岸线的可视化效果。

报错原因 #

瑞士国土面积远小于哈萨克斯坦，当前设置的最小尺子长度（25km）换算为米后（25000），在瑞士的矢量边界数据中，存在从当前点出发到后续所有点的距离都小于该尺子长度的情况，导致measure_with_ruler函数中which(pd > stick_length)[1]返回NA，触发报错。

解决办法 #

1. 缩小尺子长度范围 #

将尺子长度调整为更适合瑞士尺度的数值，比如把rulers从c(25,50,100,150,200,250)改为c(5,10,15,20,25,30)，确保最小尺子长度小于边界上相邻点的最大距离。

2. 修改measure_with_ruler函数，兼容长尺子场景 #

修改函数逻辑，当剩余路径的最大距离小于尺子长度时，直接将终点加入测量点，避免报错。修改后的函数如下：

measure_with_ruler <- function(pols, stick_length, lonlat=FALSE) {
    stopifnot(inherits(pols, "SpatVector"))
    stopifnot(length(pols) == 1)
    g <- geom(pols)[, c('x', 'y')]
    nr <- nrow(g)
    pts <- 1
    newpt <- 1
    while(TRUE) {
        p <- newpt
        j <- p:(p+nr-1)
        j[j > nr] <- j[j > nr] - nr
        gg <- g[j,]
        pd <- distance(gg[1,,drop=FALSE], gg, lonlat)
        pd <- as.vector(pd)
        i <- which(pd > stick_length)[1]
        if(is.na(i)) {
            # 当尺子过长时，直接跳到终点
            newpt <- nr
            pts <- c(pts, newpt)
            break
        }
        newpt <- i + p
        if(newpt >= nr) {
            break
        }
        pts <- c(pts, newpt)
    }
    pts <- c(pts, 1)
    g[pts,]
}

3. 修正代码笔误 #

瑞士代码的绘图部分存在笔误：rules[i]应改为rulers[i]，否则会出现变量未定义的错误。

修改后的完整瑞士代码 #

library(terra)
library(geodata)
w <- world(path=".", resolution = 3)
switz <- w[w$GID_0=="CHE",]
plot(switz)
as.data.frame(switz)
prj <- "epsg:4149"
gswitz <- project(switz, prj)
dswitz <- disagg(gswitz)
head(dswitz)
a <- expanse(dswitz)
i <- which.max(a)
a[i] / 1000000
b <- dswitz[i,]
par(mai=rep(0,4))
plot(b)

measure_with_ruler <- function(pols, stick_length, lonlat=FALSE) {
    stopifnot(inherits(pols, "SpatVector"))
    stopifnot(length(pols) == 1)
    g <- geom(pols)[, c('x', 'y')]
    nr <- nrow(g)
    pts <- 1
    newpt <- 1
    while(TRUE) {
        p <- newpt
        j <- p:(p+nr-1)
        j[j > nr] <- j[j > nr] - nr
        gg <- g[j,]
        pd <- distance(gg[1,,drop=FALSE], gg, lonlat)
        pd <- as.vector(pd)
        i <- which(pd > stick_length)[1]
        if(is.na(i)) {
            # 处理尺子过长的情况，直接跳到终点
            newpt <- nr
            pts <- c(pts, newpt)
            break
        }
        newpt <- i + p
        if(newpt >= nr) {
        break
        }
        pts <- c (pts, newpt)
    }
    pts <- c(pts, 1)
    g[pts,]
}

y <- list()
# 调整为适合瑞士的尺子长度
rulers <- c(5,10,15,20,25,30)
for (i in 1:length(rulers)) {
    y[[i]] <- measure_with_ruler(b, rulers[i]*1000)
}

par(mfrow=c(2,3), mai=rep(0,4))
for(i in 1:length(y)) {
    plot(b, col='lightgray', lwd=2)
    p <- y[[i]]
    lines(p, col='red', lwd=3)
    points(p, pch=20, col='blue', cex=2)
    # 修正笔误：rules[i]改为rulers[i]
    bar <- rbind(cbind(525000, 900000), cbind(525000, 900000-rulers[i]*1000))
    lines(bar, lwd=2)
    points(bar, pch=20, cex=1.5)
    text(525000, mean(bar[,2]), paste(rulers[i], ' km'), cex=1.5)
    text(525000, bar[2,2]-50000, paste0('(', nrow(p),')'), cex=1.25)
}

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