R语言：内存受限下基于raster/terra的农田周边林地面积计算

需求：计算县域内蓝莓田周边500米/5公里范围内的保护林地面积，数据为栅格格式。当前使用SF的st_buffer()+st_intersection()方法内存占用极高，需替换为仅基于terra/raster的低内存方案，同时避免重复计算（即重叠缓冲区覆盖的同一片林地只算一次）。

核心思路 #

用栅格距离分析替代矢量缓冲区，通过连通斑块识别解决重复计算问题：

1. 生成蓝莓田的距离栅格（记录每个栅格到最近蓝莓田的距离）
2. 提取距离≤阈值的保护林地栅格
3. 识别保护林地的连通斑块，判断每个斑块是否有部分在距离范围内，仅统计符合条件的斑块总面积（避免重复计算重叠区域）

代码实现（基于terra） #

library(terra)
library(CropScapeR)
library(httr)

# 1. 下载并加载CDL数据（复用原有逻辑）
httr::set_config(httr::config(ssl_verifypeer = 0L))
tif_file <- tempfile(fileext = '.tif')
ST_CDL <- GetCDLData(aoi = '34001', year = 2021, type = 'f', save_path = tif_file)
terra::writeRaster(ST_CDL, "ST_CDL.tif", overwrite=TRUE)
ST_CDL <- terra::rast("ST_CDL.tif")

# 2. 分类生成目标二值栅格
# 保护林地编码集合
conserved_codes <- c(63,64,141,142,143,152)
# 蓝莓田编码
blueberry_code <- 242

# 保护林地：1=保护林地，0=其他
ST_CDL_conserved <- terra::classify(ST_CDL, 
                                    rbind(cbind(conserved_codes, 1), cbind(NA, 0)), 
                                    others=0)
# 蓝莓田：1=蓝莓田，0=其他
ST_CDL_blueberries <- terra::classify(ST_CDL, 
                                      cbind(blueberry_code, 1), 
                                      others=0)

# 3. 计算蓝莓田的距离栅格（terra自动分块处理，低内存友好）
# 注：CDL默认Albers等面积投影，单位为米，无需额外转换
distance_to_blueberry <- terra::distance(ST_CDL_blueberries, target=1)

# 4. 设置距离阈值（500米，可替换为5000即5公里）
threshold <- 500
# 提取距离范围内的保护林地
conserved_in_buffer <- terra::mask(ST_CDL_conserved, distance_to_blueberry <= threshold, maskvalue=0)

# 5. 识别连通斑块，解决重复计算问题
# 给每个保护林地连通斑块分配唯一ID
conserved_patches <- terra::patches(ST_CDL_conserved, directions=8)
# 提取所有被缓冲覆盖的斑块ID
patch_ids_in_buffer <- unique(terra::extract(conserved_patches, conserved_in_buffer[conserved_in_buffer == 1], cells=TRUE)$ID)
# 生成仅包含目标斑块的栅格
target_patches <- terra::classify(conserved_patches, cbind(patch_ids_in_buffer, 1), others=0)

# 6. 计算总面积
# 自动计算每个栅格单元的实际面积（适配投影）
cell_area <- terra::cellSize(target_patches, unit="m2")
total_area <- sum(cell_area[target_patches == 1], na.rm=TRUE)

# 转换为公顷（可选）
total_area_ha <- total_area / 10000
cat("蓝莓田周边", threshold, "米范围内的保护林地总面积：", total_area_ha, "公顷\n")

关键优化说明 #

• 低内存特性：全程基于栅格操作，避免矢量多边形化、缓冲区计算的内存爆炸，terra会自动分块处理大栅格，无需一次性加载全部数据。
• 避免重复计算：通过patches()识别连通斑块，只要斑块有部分在蓝莓田缓冲范围内，就完整统计整个斑块面积；若需求改为仅统计被缓冲覆盖的斑块部分，可跳过步骤5，直接计算sum(cell_area[conserved_in_buffer == 1], na.rm=TRUE)。
• 阈值灵活调整：只需修改threshold参数即可切换500米/5公里的计算范围。

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