使用lidR包计算叶面积指数（LAI）遇两类错误的技术问询

我用R语言的lidR包基于离散回波航测LiDAR数据，遵循Beer-Lambert定律（假设叶片随机取向，Ross-Nilson间隙函数值为0.5）编写函数计算叶面积指数（LAI）。定义了低层级API函数Calc_LAI()，负责：

• 统计地面回波数量
• 统计植被回波数量
• 调用pixel_metrics()计算LAI
• 返回LAI栅格结果

遇到的两个问题 #

问题1：嵌套辅助函数lai_func()未被识别 #

将lai_func()嵌套在if(is(las,"LAS")){}代码块内时，运行Calc_LAI()会提示lai_func()未找到。但之前在其他调用catalog_apply()的低层级API中，同样的嵌套方式可正常运行。

问题2：pixel_metrics()报错重复元素/非数值指标 #

调用pixel_metrics()时出现错误：Error: Duplicated elements found. At least one of the metrics was not a number. Each metric should be a single number.。我已在冠层与地面回波比值的分子分母各加1，避免除零或对0取自然对数的情况，但仍报错。

可复现代码 #

## 加载必要包
library(lidR)

## 读取示例LiDAR数据
lasfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
MixedConifer <- readLAS(lasfile)

## 计算LAI的低层级API函数
Calc_LAI <- function(las)
{
  if (is(las, "LAScatalog"))  {
    options <- list(automerge = TRUE, need_buffer = TRUE)
    LAI <- catalog_apply(las, Calc_LAI, .options = options)
    return(LAI)
  }
  else if (is(las, "LAScluster")) {
    bbox <- st_bbox(las)
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    LAI <- Calc_LAI(las, param1, param2)
    LAI <- sf::st_crop(LAI, bbox)
    return(LAI)
  }
  else if (is(las, "LAS")) {
    las <- add_attribute(las, paste(las$gpstime, las$ReturnNumber, sep="_"), "UID")
    dem <- rasterize_terrain(las, res=1, algorithm = knnidw())
    norm <- las - dem
    
    lai_func <- function(ScanAngle, Z, UID){
      nground <- length(unique(UID[Z <= 6.6]))
      nveg <- length(unique(UID[Z >= 9.9]))
      out <- log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
      return(list(LAI=out))
    }
    
    LAI <- pixel_metrics(las=norm, func=lai_func(ScanAngle=ScanAngleRank, Z=Z, UID=UID), res=9.9)
    return(LAI)
  }
  else {
    stop("不支持该数据类型。")
  }
}

## 第一次尝试：lai_func()嵌套在Calc_LAI内
Calc_LAI(MixedConifer) # 报错：lai_func()未找到

## 第二次尝试：lai_func()移到Calc_LAI外
lai_func <- function(ScanAngle, Z, UID){
  nground <- length(unique(UID[Z <= 6.6]))
  nveg <- length(unique(UID[Z >= 9.9]))
  out <- log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
  return(list(LAI=out))
}

Calc_LAI(MixedConifer) # 报错：存在重复元素，至少一个指标不是数值

## 直接运行函数片段测试
las <- MixedConifer
las <- add_attribute(las, paste(las$gpstime, las$ReturnNumber, sep="_"), "UID")
dem <- rasterize_terrain(las, res=1, algorithm = knnidw())
norm <- las - dem

lai_func <- function(ScanAngle, Z, UID){
  nground <- length(unique(UID[Z <= 6.6]))
  nveg <- length(unique(UID[Z >= 9.9]))
  out <- log(nveg+1/(nground+1))*cos((ScanAngle*pi/180)/0.5)
  return(list(LAI=out))
}

LAI <- pixel_metrics(las=norm, func=lai_func(ScanAngle=ScanAngleRank, Z=Z, UID=UID), res=9.9) # 同样报重复元素错误

问题解决方法 #

问题1的原因与解决 #

lidR的pixel_metrics()会在独立环境中执行传入的函数，嵌套在Calc_LAI()分支内的lai_func()无法被这个环境识别。而catalog_apply()的执行环境逻辑不同，所以之前的嵌套方式能生效。

解决方法：

• 把lai_func()定义在Calc_LAI()外部，确保全局环境可见；
• 或在Calc_LAI()函数最顶部定义lai_func()，让整个函数作用域内都能访问。

问题2的原因与解决 #

错误根源：

1. 括号位置错误：原代码中log(nveg+1/(nground+1))实际计算的是log(nveg + (1/(nground+1)))，而非预期的log((nveg+1)/(nground+1))，会导致计算结果异常；
2. 向量输入未聚合：pixel_metrics()会把像素内所有点的ScanAngleRank作为向量传入函数，cos((ScanAngle*pi/180)/0.5)会生成向量，最终out是向量，而pixel_metrics()要求每个指标返回单个数值；
3. LAScluster分支参数错误：Calc_LAI(las, param1, param2)中的param1/param2未定义，会触发额外错误。

修正后的代码片段：

# 修正后的辅助函数
lai_func <- function(ScanAngle, Z, UID){
  nground <- length(unique(UID[Z <= 6.6]))
  nveg <- length(unique(UID[Z >= 9.9]))
  # 聚合扫描角为单个值（取均值）
  mean_scan <- mean(ScanAngle, na.rm = TRUE)
  # 修正对数计算的括号位置
  ratio <- (nveg + 1)/(nground + 1)
  # 防御性检查，确保比值大于0
  if (ratio <= 0) ratio <- 1e-6
  out <- log(ratio) * cos((mean_scan * pi/180)/0.5)
  return(list(LAI=out))
}

# 修正LAScluster分支的参数问题
Calc_LAI <- function(las)
{
  if (is(las, "LAScatalog"))  {
    options <- list(automerge = TRUE, need_buffer = TRUE)
    LAI <- catalog_apply(las, Calc_LAI, .options = options)
    return(LAI)
  }
  else if (is(las, "LAScluster")) {
    bbox <- st_bbox(las)
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    LAI <- Calc_LAI(las) # 移除未定义的param1/param2
    LAI <- sf::st_crop(LAI, bbox)
    return(LAI)
  }
  else if (is(las, "LAS")) {
    las <- add_attribute(las, paste(las$gpstime, las$ReturnNumber, sep="_"), "UID")
    dem <- rasterize_terrain(las, res=1, algorithm = knnidw())
    norm <- las - dem
    
    LAI <- pixel_metrics(las=norm, func=lai_func(ScanAngle=ScanAngleRank, Z=Z, UID=UID), res=9.9)
    return(LAI)
  }
  else {
    stop("不支持该数据类型。")
  }
}

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