0

我们有一个栅格，它表示对应于生长季节开始的序数日期。也就是说，栅格中的每个像素值都位于 1:365 之间，表示序号日期。

我还计算了相应年份所有 365 天的累积生长期天数。这些数据作为 365 层的栅格堆栈加载到 R 中。

我的目标是在生长季节层开始时随机抽样地理位置。然后我希望从那些相同的坐标中提取累积生长期的值，但只能从对应于季节开始像素值的生长期堆栈层中提取。

例如，如果给定像素的季节开始是一年中的第 100 天，我想在一年中的第 100 天（nlayers = 100）从该位置提取生长期天数。

我一直在尝试编写一个函数来实现这一点，但我似乎无法让它正常工作。我想最终得到一个数据框或矩阵，显示我的 x 位置、y 位置、赛季开始日和那天的 GDD。我不是在一列中有许多 GDD 值，而是得到一个 GDD 值的许多列。

看来问题出在我对提取物的使用上。我已经尝试了参数 nl、layer，并用 [[]] 索引 x 参数。它们似乎产生了相同的结果。这是一个仅需要 5 天的简化代码，以及我正在尝试构建的功能。

任何帮助/建议表示赞赏！

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

1 回答 1

2

这似乎没有引起很多关注，但我能够解决它。使用原始问题中发布的示例，最简单的解决方案是stackSelect函数。Robert Hijmans 在 R-sig-geo 上向我指出了这一点。

x <- stackSelect(GDD, SOST)

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

如果您的数据具有相同的范围和分辨率，这将非常有用。但是，我没有提到并包括我的数据没有完全对齐。因此，据我所知，我仍然需要创建一个函数。稍加思考，我就能够想出以下示例和功能并解决问题。

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
于 2013-07-26T21:47:17.900 回答