We have a raster which represents the ordinal date corresponding to the start of growing season. That is, each pixel value in the raster lies between 1:365, representing the ordinal date.
I have also calculated cumulative growing degree days for all 365 days in the corresponding year. These data are loaded into R as a raster stack with 365 layers.
My goal is to randomly sample geographic locations on the start of growing season layer. I then hope to extract the value of cumulative growing degree days from those same coordinates, but only from the growing degree days stack's layer which corresponds to the start of season pixel value.
For example, if the start of season at a given pixel was the 100th day of the year, I would like to extract the growing degree days from that location on the 100th day of the year (nlayers = 100).
I have been attempting to write a function to accomplish this, but I can't seem to get it to work right. I would like to end up with a data frame or matrix showing my x location, y location, start of season day, and GDD for that day. Instead of many GDD values in one column, I get many columns of one GDD value.
It seems the problem is in my use of extract. I've experimented with the arguments nl, layer, and indexing the x argument with [[]]. They seem to produce the same result. Here's a simplified code with only 5 days to consider, and the function I am trying to construct.
Any help/suggestions is appreciated!
#============================================================
library(raster)
SOST <- raster()
SOST[] <- 1:5
r1 <- r2 <- r3 <- r4 <- r5 <- raster()
r1[] <- 10
r2[] <- 20
r3[] <- 30
r4[] <- 40
r5[] <- 50
GDD <- stack(r1,r2,r3,r4,r5)
getGDD <- function(sost, gdd, n){set.seed(232)
samp <- sampleRandom(sost, n, xy = TRUE,
na.rm = TRUE)
df <- data.frame('x'=as.numeric(), 'y'=
as.numeric(), 'SOST'=as.numeric(),
'GDD'=as.numeric())
df.temp <- data.frame('x' = samp[1:n,1], 'y' =
samp[1:n,2], 'SOST' = samp[,3],'GDD' =
extract(gdd, samp[1:n,1:2], nl = samp[1:n,3]))
df <- rbind(df, df.temp)
return(df)
}
getGDD(sost = SOST, gdd = GDD, n = 5)
It doesn't seem like this gathered a lot of attention, but I was able to solve it. Using the sample posted in the original question, the easiest solution is the stackSelect function. This was pointed out to me by Robert Hijmans on R-sig-geo.
x <- stackSelect(GDD, SOST)
set.seed(232)
samp <- sampleRandom(SOST, 5, xy = TRUE, na.rm = TRUE)[, -3]
extract(x, samp)
This works great if your data have the same extent and resolution. However, I failed to mention and include that my data do not align perfectly. Thus, as far as I know, I still need to create a function. With a little more thought, I was able to come up with the following example and function and solve the problem.
library(raster)
#SOST and GDD simulations with same ncell and extents as actual data:
SOST <- raster(nrow = 3991, ncol = 3025, xmn = 688635, xmx = 779385,
ymn = 4276125, ymx = 4395855, crs = "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
SOST[] <- 1:5
r1 <- r2 <- r3 <- r4 <- r5 <- raster(nrow = 3951, ncol = 2995, xmn = 688620.2, xmx = 779377.8, ymn = 4276121, ymx = 4395848, crs = "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
r1[] <- 10
r2[] <- 20
r3[] <- 30
r4[] <- 40
r5[] <- 50
GDD <- stack(r1,r2,r3,r4,r5)
getGDD <- function(sost, gdd, n){
set.seed(232)
samp <- sampleRandom(sost, size = n, xy = TRUE)
extr <- NULL
for(i in 1:n){
extr[i] <- extract(gdd[[samp[i,3]]], cbind(as.matrix(samp[i,1]),
as.matrix(samp[i,2])))
}
out <- data.frame(x = samp[,1], y = samp[,2], 'SOST' = samp[,3], 'GDD' = extr)
return(out)
}
test <- getGDD(sost = SOST, gdd = GDD, n = 5)
test