匹配NLCD土地覆被等级与坐标时的不一致问题求助

在为物种出现点分配NLCD土地覆被类型时，发现相同坐标点的提取结果会随数据集长度变化而改变，疑似raster::extract函数在处理不同大小的数据集时重排了数据顺序。旧版本代码可正常运行，但当前使用raster 3.6-3版本出现该问题。

成因分析 #

1. 分块处理的顺序偏差：raster::extract处理较大数据集时会自动启用分块机制，分块过程中若未严格保留原始点的顺序，会导致提取结果与原始数据行错位。旧版本raster的分块逻辑不同，因此未出现该问题。
2. 空间对象转换的隐性变化：从tibble转换为SpatialPoints或进行投影转换时，若行索引未被正确保留，可能间接导致提取结果顺序混乱。

解决方法 #

方法1：禁用分块处理

在extract函数中设置chunk=FALSE，强制一次性处理所有点，避免分块导致的顺序错位：

# 修改data1的提取代码
data1$Cover.names1 <- raster::extract(NLCD, Datos_transformed1, chunk = FALSE)
# 修改data2的提取代码
data2$Cover.names2 <- raster::extract(NLCD, Datos_transformed2, chunk = FALSE)

方法2：通过唯一ID关联结果

为每个点添加唯一标识符，确保提取结果与原始数据行严格对应，即使顺序出现变化也能正确匹配：

# 处理data1
data1 <- data1 %>% mutate(id = row_number())
coords1 <- data1[, c("long", "lat", "id")]  
coordinates(coords1) <- ~long + lat
proj4string(coords1) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
Datos_transformed1 <- spTransform(coords1, CRS(crs_args))

# 提取时返回包含ID的数据框
extract_result1 <- raster::extract(NLCD, Datos_transformed1, df = TRUE)
# 关联回原始数据
data1 <- data1 %>% left_join(extract_result1, by = c("id" = "ID")) %>% rename(Cover.names1 = layer)

# 处理data2同理
data2 <- bind_rows(data1, data2) %>% mutate(id = row_number())
coords2 <- data2[, c("long", "lat", "id")]  
coordinates(coords2) <- ~long + lat
proj4string(coords2) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
Datos_transformed2 <- spTransform(coords2, CRS(crs_args))

extract_result2 <- raster::extract(NLCD, Datos_transformed2, df = TRUE)
data2 <- data2 %>% left_join(extract_result2, by = c("id" = "ID")) %>% rename(Cover.names2 = layer)

方法3：迁移到terra包（推荐）

raster包已被terra包替代，terra::extract函数的顺序处理更稳定，默认保留输入点的顺序：

library(tidyverse)
library(terra)

# 读取NLCD数据（替换为你的文件路径）
NLCD_terra <- rast("path/to/NLCD_2006_Land_Cover_CONUS.tif")

# 处理data1
data1 <- data1 %>% mutate(id = row_number())
# 创建空间点对象
points1 <- vect(data1, geom = c("long", "lat"), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84")
# 投影转换
points1_transformed <- project(points1, crs(NLCD_terra))
# 提取土地覆被
data1$Cover.names1 <- extract(NLCD_terra, points1_transformed)[,1]

# 处理data2
data2 <- bind_rows(data1, data2) %>% mutate(id = row_number())
points2 <- vect(data2, geom = c("long", "lat"), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84")
points2_transformed <- project(points2, crs(NLCD_terra))
data2$Cover.names2 <- extract(NLCD_terra, points2_transformed)[,1]

验证 #

修改后运行以下代码，确认相同坐标点的提取结果一致：

all.equal(data1$Cover.names1, data2$Cover.names2[1:5])
all.equal(data1$long, data2$long[1:5])
all.equal(data1$lat, data2$lat[1:5])

内容的提问来源于stack exchange，提问作者J. Lan