如何用R语言筛选处于.SHP多边形范围内的.LAZ文件？

问题描述 #

我有一个包含数千个.LAZ文件的文件夹，希望筛选出处于指定.SHP多边形范围内的文件，给符合条件的.LAZ文件添加特定命名标记。因为无法将大量文件导入GIS软件，想知道能不能用R语言实现。我试了下面的脚本但没效果，疑惑是不是必须转为.LAS格式（这会大幅增加数据量）：

# Install and load necessary packages
install.packages("lidR")
install.packages("rgdal")
library(lidR)
library(rgdal)

# Read the .shp file with the extent
extent <- readOGR(dsn = "/path/to/folder", layer = "filename")

# Read the .laz files in the folder
folder <- "/path/to/folder"
laz_files <- list.files(folder, pattern = "\\.laz$", full.names = TRUE)

# Iterate through the .laz files and select those that overlap with the extent
for (laz_file in laz_files) {
  points <- readLAS(laz_file)
  if (extent(points) %over% extent) {
    cat("File", laz_file, "overlaps with the extent.\n")
    # Create a new file name with the selected file name appended
    new_name <- sub(".laz", paste0("_selected_", basename(laz_file)), laz_file)
    # Save the file under the new name
    writeLAS(points, new_name)
  } else {
    cat("File", laz_file, "does not overlap with the extent.\n")
  }
}

解决方案 #

首先明确：完全不需要转为LAS格式，lidR原生支持读写LAZ文件（依赖系统安装的LASzip，通常lidR安装时会自动配置），不会额外增加数据量。原脚本失效主要是逻辑错误和效率问题，以下是修正后的实现方案：

关键问题分析 #

1. 原脚本读取整个LAZ文件来获取范围，对于数千个文件来说效率极低，只需读取文件头部就能获取范围信息。
2. 空间交集判断逻辑错误，extent(points)的bbox对象与多边形用%over%对比不适用，需用sf包的空间函数实现准确判断。
3. 重命名逻辑错误，会导致文件名格式混乱。
4. rgdal已被弃用，建议用更稳定的sf包处理空间数据。

修正后的代码 #

# 安装并加载依赖包
install.packages(c("lidR", "sf"))
library(lidR)
library(sf)

# 读取目标SHP多边形
target_polygon <- st_read(dsn = "/path/to/your/shp/folder", layer = "your_shp_layer_name")

# 获取文件夹下所有LAZ文件的完整路径
laz_folder <- "/path/to/your/laz/folder"
laz_files <- list.files(laz_folder, pattern = "\\.laz$", full.names = TRUE)

# 批量处理每个LAZ文件
for (laz_path in laz_files) {
  # 仅读取LAZ文件头部，快速获取点云范围（无需加载全部点云）
  las_header <- readLASheader(laz_path)
  las_bbox <- st_bbox(las_header)
  las_extent_poly <- st_as_sfc(las_bbox)
  
  # 判断LAZ范围与目标多边形是否存在交集
  has_overlap <- st_intersects(las_extent_poly, target_polygon, sparse = FALSE)[1, 1]
  
  if (has_overlap) {
    cat("处理文件：", laz_path, "→ 符合条件\n")
    # 生成带标记的新文件名：原文件名_selected.laz
    file_dir <- dirname(laz_path)
    base_name <- tools::file_path_sans_ext(basename(laz_path))
    new_laz_path <- file.path(file_dir, paste0(base_name, "_selected.laz"))
    
    # 直接复制原文件（仅标记无需修改点云时，此方式效率最高）
    file.copy(laz_path, new_laz_path, overwrite = TRUE)
    
    # 如果需要裁剪点云到多边形范围内再保存，取消以下注释：
    # las_data <- readLAS(laz_path)
    # las_cropped <- clip_roi(las_data, target_polygon)
    # writeLAS(las_cropped, new_laz_path)
  } else {
    cat("处理文件：", laz_path, "→ 不符合条件\n")
  }
}

核心优化点 #

• 效率提升：用readLASheader替代readLAS，仅读取文件元数据获取范围，速度提升显著，适合批量处理大量文件。
• 准确判断：通过sf包的st_intersects实现空间交集判断，逻辑更严谨。
• 无损处理：直接复制原LAZ文件实现标记，或裁剪后保存为LAZ格式，全程无需转为LAS，避免数据量膨胀。
• 规范命名：新文件名格式为原文件名_selected.laz，清晰易识别。

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