在R中为GPS点位提取对应年份及前后年的USDA CropScape作物值

1. 确认数据结构与预处理 #

假设你的GPS数据框为gps_data，包含：

• 日期列date（需为Date类型）
• 经纬度列lon、lat
  预处理后的CropScape栅格已按年份命名（如cropscape_2014至cropscape_2021），先将所有栅格存入列表统一管理：

# 把各年份栅格存入列表，命名与年份一致
cropscape_list <- list(
  "2014" = cropscape_2014,
  "2015" = cropscape_2015,
  "2016" = cropscape_2016,
  "2017" = cropscape_2017,
  "2018" = cropscape_2018,
  "2019" = cropscape_2019,
  "2020" = cropscape_2020,
  "2021" = cropscape_2021
)

2. 提取GPS点位的日历年度 #

用lubridate包从日期中提取年份（未安装则先执行install.packages("lubridate")）：

library(lubridate)
gps_data$cal_year <- year(gps_data$date)

3. 批量提取对应年份作物值 #

方式一：用dplyr逐行处理（简洁直观） #

library(dplyr)
library(raster)

gps_data <- gps_data %>%
  rowwise() %>%
  mutate(
    # 前一年作物值：仅当年份在2014-2021范围内时提取，否则设为NA
    crop_calYr_minus1 = ifelse(cal_year - 1 >= 2014, 
                              extract(cropscape_list[[as.character(cal_year - 1)]], 
                                      c(lon, lat)), 
                              NA),
    # 当前日历年度作物值
    crop_calYr = extract(cropscape_list[[as.character(cal_year)]], c(lon, lat)),
    # 后一年作物值：2021年点位直接设为NA，其余年份正常提取
    crop_calYr_plus1 = ifelse(cal_year + 1 <= 2021, 
                             extract(cropscape_list[[as.character(cal_year + 1)]], 
                                     c(lon, lat)), 
                             NA)
  ) %>%
  ungroup()

方式二：用for循环处理（适合习惯循环逻辑的场景） #

# 先初始化三列为NA
gps_data$crop_calYr_minus1 <- NA
gps_data$crop_calYr <- NA
gps_data$crop_calYr_plus1 <- NA

# 逐行遍历GPS数据
for(i in 1:nrow(gps_data)){
  curr_yr <- gps_data$cal_year[i]
  curr_lon <- gps_data$lon[i]
  curr_lat <- gps_data$lat[i]
  
  # 提取前一年值
  if(curr_yr - 1 >= 2014){
    gps_data$crop_calYr_minus1[i] <- extract(cropscape_list[[as.character(curr_yr - 1)]], c(curr_lon, curr_lat))
  }
  
  # 提取当前年份值
  gps_data$crop_calYr[i] <- extract(cropscape_list[[as.character(curr_yr)]], c(curr_lon, curr_lat))
  
  # 提取后一年值，2021年点位跳过设为NA
  if(curr_yr + 1 <= 2021){
    gps_data$crop_calYr_plus1[i] <- extract(cropscape_list[[as.character(curr_yr + 1)]], c(curr_lon, curr_lat))
  }
}

大数据量优化方案 #

如果GPS点位数量多，逐行提取效率低，可转为空间点对象批量提取：

# 转为空间点对象（确保与栅格投影一致）
gps_sp <- SpatialPointsDataFrame(
  coords = gps_data[, c("lon", "lat")], 
  data = gps_data,
  proj4string = crs(cropscape_list[["2014"]])
)

# 批量提取所有年份的作物值
for(yr in 2014:2021){
  gps_data[[paste0("crop_", yr)]] <- extract(cropscape_list[[as.character(yr)]], gps_sp)
}

# 匹配到目标三列并清理中间列
gps_data <- gps_data %>%
  mutate(
    crop_calYr_minus1 = .data[[paste0("crop_", cal_year - 1)]],
    crop_calYr = .data[[paste0("crop_", cal_year)]],
    crop_calYr_plus1 = ifelse(cal_year + 1 <=2021, .data[[paste0("crop_", cal_year + 1)]], NA)
  ) %>%
  select(-starts_with("crop_20"))

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