如何用R批量读取GeoTIFF并提取像素经纬度及波段值

前置准备 #

先安装并加载所需R包：

# 安装包（首次运行需要执行）
install.packages(c("terra", "dplyr", "purrr", "tidyr", "stringr"))

# 加载包
library(terra)
library(dplyr)
library(purrr)
library(tidyr)
library(stringr)

1. 批量获取并解析文件信息 #

定位目标文件夹，提取所有GeoTIFF文件的路径、拍摄日期和波段信息：

# 设置文件夹路径
tif_dir <- "D:/mytiffs"

# 获取所有tif文件的完整路径
tif_files <- list.files(tif_dir, pattern = "\\.tif$", full.names = TRUE)

# 从文件名中解析日期和波段
file_metadata <- tibble(file_path = tif_files) %>%
  mutate(
    file_name = basename(file_path),
    capture_date = str_extract(file_name, "\\d{4}-\\d{2}-\\d{2}"),
    band_code = str_extract(file_name, "B\\d{2}")
  )

2. 批量读取影像并提取像素经纬度与波段值 #

编写单个文件处理函数，批量读取影像、转换为经纬度投影、提取像素级数据：

# 定义单文件处理逻辑
process_single_tif <- function(file_path, capture_date, band_code) {
  # 读取GeoTIFF影像
  img <- rast(file_path)
  # 转换为WGS84经纬度投影（EPSG:4326），若原影像已是该投影可删除此步
  img_wgs84 <- project(img, "EPSG:4326")
  # 提取像素坐标（x=经度, y=纬度）和波段值
  pixel_data <- as.data.frame(img_wgs84, xy = TRUE)
  # 重命名波段列并添加日期信息
  pixel_data <- pixel_data %>%
    rename(!!band_code := 3) %>%  # 将默认波段列名替换为BXX
    mutate(date = capture_date)
  return(pixel_data)
}

# 批量处理所有文件，合并为单个数据框
raw_pixel_data <- pmap_dfr(file_metadata, process_single_tif)

3. 按拍摄日期整合多波段数据 #

将长格式数据转换为宽格式，每个日期对应所有波段的数值，最终输出符合要求的data.frame：

final_dataframe <- raw_pixel_data %>%
  pivot_wider(
    id_cols = c(x, y, date),
    names_from = band_code,
    values_from = all_of(band_code)
  ) %>%
  rename(longitude = x, latitude = y)  # 重命名坐标列，更直观

注意事项 #

• 确保所有影像的空间范围、分辨率一致，否则合并时会出现数据不匹配；
• 处理大尺寸影像时，可考虑分块读取（参考terra包的readStart/readStop方法），避免内存溢出；
• 若原影像投影已为WGS84，删除project步骤可大幅提升处理速度。

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