4

I am using a moving window average to estimate change between historical climate (here, a 36-year average for a period) and the moving window averages of climate into the future. Each moving window period should start +1 year at each iteration and should be a 36-year average. My raster stack for the "future" period for which moving windows will be calculated looks something like this:

library(raster)
library(rgdal)
library(gdm)

b <- brick(ncol=5, nrow=5, nl=131)
values(b) <- 1:(25*131)
yr<-seq(1970, 2100, 1)
setZ(b, yr)

   class      : RasterBrick
dimensions : 5, 5, 25, 131  (nrow, ncol, ncell, nlayers)
resolution : 72, 36  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ...
min values :       1,      26,      51,      76,     101,     126,     151,     176,     201,      226,      251,      276,      301,      326,      351, ...
max values :      25,      50,      75,     100,     125,     150,     175,     200,     225,      250,      275,      300,      325,      350,      375, ...
time       : 1970, 2100 (min, max)

b2 <- brick(ncol=5, nrow=5, nl=131) values(b2) <- 1

I need to define the starting and ending periods and window lengths for the moving windows:

start.year<-1970 
end.year.hist<-2005 #(the last year in the "historical" 36-year period)
end.year.future<-2100
window.length<-end.year.hist-start.year+1 #(+1 to be inclusive of the starting year)
> window.length
[1] 36 #confirm the number of windows to be correct

I then define the first period (in other words the historical period that I base the size of my moving window on, which I've called "focal years"):

focalYears<-start.year:end.year.hist
> focalYears
 [1] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
[16] 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
[31] 2000 2001 2002 2003 2004 2005

Based on the the above, this is the number of windows that need to be calculated:

num.windows<-end.year.future-end.year.hist
[1] 95

And here is my loop:

    for(i in 1:length(num.windows)){
      
      futureYears<-focalYears+1 #Define the moving window period as "futureYears", which 
#is the historical window period +1 (i.e. the first iteration should be for 1971 -- 2006)
          
#Make the moving window average for the future (i.e. the mean of the 36-year future window)

      future.Rast.pred<-mean(b[[futureYears]])
      
 

   #Predict function in `gdm`, which uses the historical raster (`current.Rast.mask`) and calculates the change given the future raster, 
#being the moving window average defined above. The historical period raster is simply the average of 1970--2005. 
     

     timePred <- predict(G15.gdm.bio1, data=setNames(b2[[i]],"bio1"), time=TRUE,
                              predRasts=setNames(b[[i]],"bio1"))
      
    rm(timePred)
    }

However, I get the error:

Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'mean': not a valid subset

What is the issue with how I've tried to subset to calculate the mean period here?

P.S. However, If I change focalYears to focalYears<-seq(1,36,1)

> focalYears
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34 35 36

Then I don't get the error about the mean anymore and the loop does run, but it only iterates once for the first year and then stops.

(P.P.S. b2 is an historical raster brick which is the average of the historical period for which the gdm results were based on i.e. the 36-year period. The number of layers needs to match the future period, which is why there are 131 despite the historical period only being 36 years.)

 > b2
class      : RasterBrick
dimensions : 5, 5, 25, 131  (nrow, ncol, ncell, nlayers)
resolution : 72, 36  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ...
min values :       1,       1,       1,       1,       1,       1,       1,       1,       1,        1,        1,        1,        1,        1,        1, ...
max values :       1,       1,       1,       1,       1,       1,       1,       1,       1,        1,        1,        1,        1,        1,        1, ...

1 Answer 1

3

Your code (the loop) does not run and does not show the error. Without that is very difficult to help, but I will try.

The error occurs here

future.Rast.mask[[futureYears]]

Because you are using years (1970, 1972, 1973, ...) as an index. You would need to have that many (>1970) layers for that to work. It should be an index like (1, 2, 3, 4, ...). That is why it works with focalYears <-seq(1,36,1). If your index goes outside the allowed range you get an error like this

b <- brick(system.file("external/rlogo.grd", package="raster"))
b[[3:4]]
#Error in .local(x, ...) : not a valid subset

The loop running only once is because you do

for(i in 1:length(num.windows))

The length of num.windows is 1. You are looking for

for(i in 1:num.windows)
4
  • Ok sure, will do. I think my loop looks incomplete because I added comments in-between the different parts of the loop. I'll try clarify this question with something reproducible. Thanks. Dec 7, 2020 at 6:59
  • Thank you for trying to help this asker with writing their code snippet. To avoid flags on your answer (for being Not An Answer; I've already seen one) please try to turn your "long comment in an answer" into an "answer to the question" as soon as you can.
    – PolyGeo
    Dec 7, 2020 at 11:22
  • 1
    Once you know how to ask an R question, and take the time to do it, you can answer most of them yourself. Dec 7, 2020 at 16:47
  • @RobertHijmans Sorry, of course its not reproducible because I haven't provided the output of G15.gdm.bio1. But, as you point out, the issues were all in the first part of the loop. Thank you, your answer is super helpful. Dec 7, 2020 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.