我无法重新创建栅格数据。我只是下载了几年的Fremont County，ID数据，因为它看起来像是您的点所在的位置。此处我使用了raster包，但我更愿意使用terra。

library(raster)

l1 <- list.files("Downloads/polygonclip", pattern="tif$", full.names
=T)
cd1 <- raster::stack(l1)

# Should be able to change the categories in the stack
# so that the extract returns character value, e.g. 
# "barley" instead of 152 with something like this. 
# Get levels
ca1 <- levels(cd1) 
# MOdify to get desired factors
ca2 <- ca1[[1]][,c(1,5)]
# Assign
levels(cd2, layer=c(1:4)) <- ca2

# ... But that way doesn't work for me for some reason, 
# so I did it this way. You'll need to do this for each
# year/layer in the stack
levels(cd1[[1]]) <- ca2
levels(cd1[[2]]) <- ca2
levels(cd1[[3]]) <- ca2
levels(cd1[[4]]) <- ca2

# Recreate the points data shared
pts1 <- new("SpatialPointsDataFrame", data = structure(list(id = c("GT74", 
"GT48", "GT52", "GT82", "GT74", "GT52", "GT82", "GT9", "GT45", 
"GT43"), date = structure(c(1607320849, 1544097667, 1588500050, 
1554854505, 1618632063, 1549476050, 1577080837, 1449784850, 1542376850, 
1548309669), class = c("POSIXct", "POSIXt"), tzone = "MST7MDT"), 
    id_bioYr = c("GT74_2020", "GT48_2018", "GT52_2019", "GT82_2018", 
    "GT74_2020", "GT52_2018", "GT82_2019", "GT9_2015", "GT45_2018", 
    "GT43_2018"), state = c(1L, 1L, 3L, 2L, 1L, 1L, 2L, 2L, 2L, 
    1L), used = c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0), landcover = c("agricultural", 
    "agricultural", "agricultural", "agricultural", "agricultural", 
    "shadow", "water", "xericGrass", "agricultural", "agricultural"
    ), calYr = c(2020, 2018, 2020, 2019, 2021, 2019, 2019, 2015, 
    2018, 2019)), row.names = c(471519L, 253950L, 310790L, 49469L, 
479572L, 293620L, 573180L, 627720L, 238551L, 151426L), class = "data.frame"), 
    coords.nrs = numeric(0), coords = structure(c(483309.821132995, 
    480304.323982496, 479889.104783906, 469493, 458709.035703965, 
    475920.645473197, 470714.848866241, 474984.326637862, 473390.631355455, 
    479941.915819768, 4866298.15517148, 4867412.65105531, 4865629.72751859, 
    4865043, 4864209.38410713, 4863155.30805073, 4865581.17021548, 
    4865450.06432659, 4865829.2649267, 4861889.90054007), .Dim = c(10L, 
    2L), .Dimnames = list(NULL, c("coords.x1", "coords.x2"))), 
    bbox = structure(c(458709.035703965, 4861889.90054007, 483309.821132995, 
    4867412.65105531), .Dim = c(2L, 2L), .Dimnames = list(c("coords.x1", 
    "coords.x2"), c("min", "max"))), proj4string = new("CRS", 
        projargs = "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))

# Removed one of the points because it was for the year 2015 
# which threw everything off. 
pts1 <- pts1[-8,]

# Need same projection for extract        
pts1 <- spTransform(pts1, crs(cd1))

# Set Z dimension (years) for matching in next function
# Note the years used
cd1 <- setZ(cd1, 2018:2021, "years")

# This takes a year and returns the indices of that year
# from the list of years. The indices are for the year
# prior to the selected year, the year selected, and the
# year after selected year.
# Note the years used. 2018 was my oldest, but I think
# your data goes to 2015, so you'll have to change each
# 2018
rysel <- function(yind) {
    if(yind==2018) {
        yex = c(1, 2)
    } else if(yind==2021) {
        yex = c(3, 4)
    } else if(2018 < yind & yind < 2021) {
        yex = c(yind-2018, yind-2018+1, yind-2018+2)
    }
    return(yex)
}

# This takes the index of the row in the points data
# and uses the rysel function to subset the stack and 
# then extract.
extryr <- function(x) {
pty = pts1$calYr[x]

ytmp = rysel(pty)

yextr = extract(cd1[[ytmp]], pts1[x,])

# Pad the extract for either end
# Note the years used
if(pty==2018){
    yextr = c(yextr, "NA")
} else if(pty==2021) {
    yextr = c("NA",yextr)
}

return(yextr)
}

# Data frame with column labels for holding extract values
df2 <- data.frame(calYrprior = NA, calYrcur = NA, calYrnext = NA)

# 
for(i in 1:nrow(pts1)) df2[i,] = extryr(i)

# Convert the values in df2 to numeric
# You can ignore the warning
df2 <- sapply(df2, as.numeric)
Warning messages:
1: In extryr(i) : NAs introduced by coercion
2: In extryr(i) : NAs introduced by coercion
3: In extryr(i) : NAs introduced by coercion

> df2
      calYrprior calYrcur calYrnext
 [1,]        152      152       152
 [2,]         21       21        NA
 [3,]        152      152       152
 [4,]        152      152        21
 [5,]         NA       23        43
 [6,]          0        0         0
 [7,]          0        0         0
 [8,]        152      152        NA
 [9,]         NA       NA        NA

# Combine the pts1 data and the extract data
pts1 <- cbind(pts1, df2)

你需要修改我标注的年份。如果栅格中的年份和点中的年份有差异，这就不起作用了。例如，如果堆栈中有2014、2016-2019，而点覆盖了2014-2019年，你就必须修改它。
提取仍然返回数值而不是字符/因子值。但是您可以根据需要执行替换来更改这些值。

正如约翰·波罗所建议的，使用"terra"（代替"raster"）要容易得多。
示例数据

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
x <- c(r*1, r*2, r*3, r*4, r*5)
time(x, "years") <- 2001:2005
names(x) <- paste0("Y", 2001:2005)

d <- data.frame(lon=c(5.867235,  6.319561), 
       lat=c(49.85114, 49.55199), year=c(2003, 2004))

溶液

# match the year in your data.frame with that of the SpatRaster layers
i <- match(d$year, time(x))
## or 
# i <- match(d$year, 2001:2005)

如果只需要匹配年份的值，则可以执行以下操作

ey <- extract(x, d[, c("lon","lat")], layer=i)
ey
#  ID layer value
#1  1 Y2003  1473
#2  2 Y2004   692

但也要得到前后的年份：

## extract values for all years
e <- extract(x, d[, c("lon","lat")], ID=FALSE)

## find the layers (years) for each point
j <- cbind(rep(1:length(d$year), each=3), 
           rep(i, each=3) + c(-1:1))

## subset the extracted data
e[j]
#[1]  982 1473 1964  519  692  865

## to keep track of what is what
b <- data.frame(point=j[,1], year=time(x)[j[,2]], value=e[j])
b
#  point year value
#1     1 2002   982
#2     1 2003  1473
#3     1 2004  1964
#4     2 2003   519
#5     2 2004   692
#6     2 2005   865

reshape(b, direction="wide", idvar="point", timevar="year")
#  point value.2002 value.2003 value.2004 value.2005
#1     1        982       1473       1964         NA
#4     2         NA        519        692        865

这是几天前的一个类似（但更复杂）的问题。