Extract Value from Raster Stack Layer Determined by Different Raster's Pixel Value

Hello,

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!
Andrew

#============================================================
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)



--
*Andrew P. Vitale*
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632
vitale232 at gmail.com

        [[alternative HTML version deleted]]

After a little more effort, I was able to get my function to run.  I'd still
be interested in a simpler solution, but here's what I came up with.

Thanks again, Robert.

What works:


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



On Thu, Jul 25, 2013 at 2:00 PM, Andrew Vitale <vitale232 at gmail.com> wrote:

Robert and list,

Thanks for such a quick, simple solution.

I should have included one more detail in the original post.  The GDD and
SOST objects with my actual data do not have the same extents or number of
cells.  My GDD data are derived from reprojecting and diaggregating a 1000m
resolution raster to as close to 30 as possible.  I then crop the GDD stack
to the SOST extent, and end up with pretty different extents and cell
counts.


Hence, I'm unable to use the stackSelect() solution without more
(forceful) data manipulation.  Do you have another suggestion?

Here's my example again, but with extents and cell numbers representing my
actual data.  Sorry these important details did not come out in my original
post.


#============================================================================
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)


#As I'm sure you've guessed, your solution doesn't work with my original
data due to differing extents.
x <- stackSelect(GDD, SOST)

set.seed(232)
samp <- sampleRandom(SOST, 5, xy = TRUE, na.rm = TRUE)[, -3]
extract(x, samp)


On Wed, Jul 24, 2013 at 6:00 PM, Robert J. Hijmans <r.hijmans at gmail.com>
wrote:

Andrew,
I think you can use the stackSelect function for this to shortcut the
problem:

# your data
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)
#

x <- stackSelect(GDD, SOST)

set.seed(232)
samp <- sampleRandom(SOST, 5, xy = TRUE, na.rm = TRUE)[, -3]
extract(x, samp)

Best, Robert

On Wed, Jul 24, 2013 at 4:02 PM, Andrew Vitale <vitale232 at gmail.com>
wrote:

Hello,

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!
Andrew

#============================================================
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)



--
*Andrew P. Vitale*
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632
vitale232 at gmail.com

        [[alternative HTML version deleted]]

--
Andrew P. Vitale
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632
vitale232 at gmail.com

--
Andrew P. Vitale
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632
vitale232 at gmail.com