对于处理栅格数据，我建议使用terra包。它是raster包的更新，速度更快，在大多数情况下更易于使用（你可以用install.packages('terra')安装它）。terra有一组apply函数，允许你在每一堆像素上应用函数。在下面的例子中，我使用一些虚拟数据创建了一个2x2像素的光栅，并使用terra::app函数执行您在问题中概述的过程：

# Create random data for 181 rasters.
x = seq(-0, 10, length.out=181)

# `dnorm` can be used to simulate a bell curve, like the NDVI data. 
vals = rep(dnorm(x, 5, 1), 4)

# I couldn't get the raster to reshape correctly so I am just taking the first
# 50 values so the data are a nice bell curve. 
arr <- array(vals, dim = c(2, 2, 50)) 

# Replace this with your data. From your code, you can just use:
# bn <- terra::rast(fn)
r <- terra::rast(arr)

# Define a custom function that operates on one stack of pixels at a time.
fill_max <- function(timeseries) {
  
  # Find the index with the max value
  idx = which.max(timeseries)
  
  # Set all values past the max to the max value
  timeseries[idx:length(timeseries)] <- timeseries[idx]
  
  return(timeseries)
}

# Apply the function to the raster.
out <- terra::app(r, fun=fill_max)

# Check that process worked. Take the timeseries of the first pixel. 
vals <- terra::values(out)[1,]
# Plot the timeseries. 
plot(1:length(vals), vals)

创建于2023年1月18日，使用reprex v2.0.2

如果希望曲线永远不会下降，可以使用terra::cummax
示例数据

library(terra)
r <- rast(ncol=10, nrow=10, nlyr=30)
set.seed(1)
values(r) <- runif(size(r))

溶液

x <- cummax(r)

图示为四个单元格

par(mfrow=c(2,2), mar=rep(2,4))
for (i in c(2,3,8,9)) {
  plot(unlist(r[i]))
  lines(unlist(x[i]), col="blue", lwd=2)
}