使用raster::extract计算点周边比例土地覆盖时Buffer功能失效

问题原因 #

核心问题是raster包的extract()函数对terra的SpatRaster/SpatVector对象兼容性有限，传入Spat类型对象时，buffer参数未被正确解析，导致仅提取了点所在的单个像元值，而非缓冲区内的所有像元。

解决方案1：转换为raster包原生对象 #

将Spat类型对象转为raster包的原生类型后再提取：

library(dplyr)
library(rnaturalearth)
library(sf)
library(FedData)
library(raster)

# 基础地图
state = ne_states(iso_a2 = "US", returnclass = "sf") %>%    
  filter(iso_3166_2 == "US-MA") %>%                         
  st_transform(crs = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

# 点数据
Lon<-c(-69.97623, -69.97674)
Lat<-c(41.90238, 41.90267)
x<-cbind(Lon,Lat)
v<-vect(x, crs="+proj=longlat")

# 获取NLCD数据
nlcd = get_nlcd(
  template = state,
  label = "4pland",
  year = 2011,
  force.redo = F
) 

# 转换对象类型
r_nlcd <- raster(nlcd)  # SpatRaster转RasterLayer
r_p <- as(v, "Spatial") # SpatVector转SpatialPoints
r_p <- spTransform(r_p, crs(r_nlcd)) # 对齐坐标系

# 提取缓冲区内所有像元
ex_ras <- extract(r_nlcd, r_p, buffer = 1000, df = TRUE)

# 计算每个点的土地覆盖比例
cover_proportions <- ex_ras %>%
  group_by(ID) %>%
  count(layer) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup()

解决方案2：使用terra包原生函数（推荐） #

terra是raster的升级替代，原生支持Spat对象，extract()对buffer参数的支持更稳定：

library(dplyr)
library(rnaturalearth)
library(sf)
library(FedData)
library(terra)

# 基础地图
state = ne_states(iso_a2 = "US", returnclass = "sf") %>%    
  filter(iso_3166_2 == "US-MA") %>%                         
  st_transform(crs = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

# 点数据
Lon<-c(-69.97623, -69.97674)
Lat<-c(41.90238, 41.90267)
x<-cbind(Lon,Lat)
v<-vect(x, crs="+proj=longlat")

# 获取NLCD数据
nlcd = get_nlcd(
  template = state,
  label = "4pland",
  year = 2011,
  force.redo = F
) 

# 对齐点与栅格坐标系
p <- project(v, crs(nlcd))

# 提取缓冲区内所有像元
ex_terra <- extract(nlcd, p, buffer = 1000)

# 转换为数据框并计算比例
cover_proportions <- ex_terra %>%
  as.data.frame() %>%
  rename(ID = ID, land_cover = nlcd_2011_landcover) %>%
  group_by(ID) %>%
  count(land_cover) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup()

关键说明 #

• terra的extract()默认返回列表，每个元素对应一个点的缓冲区内所有像元值，转成数据框后即可用dplyr统计。
• 两种方案都必须确保点与栅格坐标系一致，这是准确提取的前提。
• 若要排除NoData像元，可在统计时添加filter(!is.na(land_cover))。

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