[R-sig-Geo] Extract Value from Raster Stack Layer Determined by Different Raster's Pixel Value

Robert J. Hijmans r.hijmans at gmail.com
Sat Jul 27 02:38:15 CEST 2013 

Andrew,

To avoid looping you could do:

names(SOST) <- 'SOST'
set.seed(232)
n = 5
samp <- sampleRandom(SOST, size = n, xy = TRUE)
e <- extract(GDD, samp[,1:2])[cbind(1:nrow(samp),samp[,3])]
out <- data.frame(samp, 'GDD' = e)
out

Robert


On Fri, Jul 26, 2013 at 2:40 PM, Andrew Vitale <vitale232 at gmail.com> wrote:
> 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]]
>>> >
>>> > _______________________________________________
>>> > R-sig-Geo mailing list
>>> > R-sig-Geo at r-project.org
>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>> --
>> 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

More information about the R-sig-Geo mailing list