¿Cómo puedo obtener dos rásteres que den el valor máximo para cada celda de la cuadrícula por año y también den las fechas en que se produjo ese valor máximo? A continuación se muestra el ejemplo reproducible con algunos pasos que he implementado.

library(raster)

# Create a raster

r1 <- raster(nrow=10, ncol=7)
r <- stack(setValues(r1, runif(ncell(r1))),
           setValues(r1, runif(70 ,0.6,0.9)),
           setValues(r1, runif(70 ,0.2,0.4)),
           setValues(r1, runif(70 ,1,2)),
           setValues(r1, runif(70 ,0.5,1.0)),
           setValues(r1, runif(70 ,0.3,0.9)),
           setValues(r1, runif(70 ,1,2)))
r

# Make Dates. This is random, i have about 24000 values.

Dates<-data.frame(Date=c("2000-01-02","2000-01-03","2000-02-03",
           "2001-09-02","2001-09-03","2001-10-01",
           "2001-10-02"))

Date_val<-as.Date(Dates$Date,format="%Y-%m-%d")
Date_val

r.dt<-setZ(r,Date_val)

# Get indices to make annual maxima value for each grid cells

indices <- format(as.Date(getZ(r.dt), format = "%Y-%b-%d"), format = "%Y")

# Implement stackApply to get maximum value each year.

rmax<-stackApply(r.dt,indices = indices,fun=max,na.rm=T)
plot(rmax)