使用terra::buffer后调用terra::intersect出现几何拓扑错误求助

问题与错误详情 #

使用terra::intersect()时触发错误：

Error: IllegalArgumentException: point array must contain 0 or >1 elements

确认该错误由terra::buffer()操作导致——经buffer处理的对象，调用terra::intersect()、terra::as.polygons()或terra::is.valid()均会报错，但可正常绘图。另有测试者运行代码出现拓扑错误：

Error: TopologyException: Input geom 0 is invalid: Nested shells at -120.40542544090771 34.654180422003208

可复现代码 #

数据获取（Linux环境） #

library(CropScapeR)
httr::set_config(httr::config(ssl_verifypeer = 0L))
tif_file <- tempfile(fileext = '.tif')
# 下载加州圣巴巴拉县2021年CDL数据
ST_CDL <- GetCDLData(aoi = '06083', year = 2021, type = 'f', save_path = tif_file) 
terra::writeRaster(ST_CDL, "ST_CDL.tif", overwrite=TRUE)
ST_CDL = terra::rast("ST_CDL.tif")

问题复现代码 #

# 分类自然土地
nat_lands_codes <- rbind(c(63, 63), c(64, 64), c(141, 141), c(142, 142), c(143, 143), c(152, 152))
ST_CDL_nat_lands <- terra::classify(ST_CDL, nat_lands_codes, others=NA)

# 分类作物
aux <- c(1:62, 65:140, 144:151, 153:255)
aux <- t(t(aux))
crops_codes <- cbind(aux, aux)
ST_CDL_crop <- terra::classify(ST_CDL, crops_codes, others=NA)

# 转为多边形
ST_CDL_nat_lands <- terra::as.polygons(ST_CDL_nat_lands)
ST_CDL_crop <- terra::as.polygons(ST_CDL_crop)

# 投影到WGS84
ST_CDL_nat_lands <- terra::project(ST_CDL_nat_lands, "EPSG:4326")
ST_CDL_crop <- terra::project(ST_CDL_crop, "EPSG:4326")

# 生成缓冲（耗时5-10分钟）
ST_CDL_crop_buffer <- terra::buffer(ST_CDL_crop, 500)

# 执行相交操作触发错误
intersection <- terra::intersect(ST_CDL_crop_buffer, ST_CDL_nat_lands)

核心原因与解决方法 #

1. 投影顺序错误（根本问题） #

buffer操作需要在**投影坐标系（平面坐标系，单位为米）**下执行，原数据NAD83的单位是米，你却先将数据转到WGS84（地理坐标系，单位为度），再设置buffer宽度500——这实际是缓冲500度，会生成极端异常的几何图形，直接导致拓扑错误。

修正操作顺序：先在原投影下做buffer（500米），再投影到WGS84：

# 先对原投影的多边形做buffer（500米）
ST_CDL_crop_buffer <- terra::buffer(ST_CDL_crop, 500)
# 再投影到WGS84
ST_CDL_crop_buffer <- terra::project(ST_CDL_crop_buffer, "EPSG:4326")

2. 修复拓扑异常 #

即使调整顺序，buffer操作仍可能生成无效几何（如零面积多边形、嵌套壳），可通过terra::makeValid()修复：

# 修复两个对象的拓扑
ST_CDL_crop_buffer <- terra::makeValid(ST_CDL_crop_buffer)
ST_CDL_nat_lands <- terra::makeValid(ST_CDL_nat_lands)

# 再执行相交操作
intersection <- terra::intersect(ST_CDL_crop_buffer, ST_CDL_nat_lands)

3. 简化分类代码（可选优化） #

原分类代码可简化，避免繁琐的矩阵拼接：

# 自然土地分类
nat_lands_codes <- c(63,64,141,142,143,152)
ST_CDL_nat_lands <- terra::classify(ST_CDL, function(x) ifelse(x %in% nat_lands_codes, x, NA))

# 作物分类（取自然土地以外的所有类别）
crop_codes <- setdiff(1:255, nat_lands_codes)
ST_CDL_crop <- terra::classify(ST_CDL, function(x) ifelse(x %in% crop_codes, x, NA))

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