3

I have some rasters that I want to compute the moving average and then rewrite the rasters again.

  r1 <- r2 <- r3 <- r4 <- r5 <- r6 <- raster(nrows=10, ncols=10);
  r1 <- setValues(r1,runif(100,min=1,max=100));
  r2 <- setValues(r2,runif(100,min=1,max=100));
  r3 <- setValues(r3,runif(100,min=1,max=100));
  r4 <- setValues(r4,runif(100,min=1,max=100));
  r5 <- setValues(r5,runif(100,min=1,max=100));
  r6 <- setValues(r6,runif(100,min=1,max=100));
  # Stack them
  st1 <- stack(r1,r2,r3,r4,r5,r6)
  #compute the moving average:
  x <- calc(st1, function(x) movingFun(x, 3, mean))

Now the layers in x should be recalculated as:

  firstlayer=(first layer + fourth layer)/2
  secondlayer=(second layer + fifth layer)/2
  thirdlayer=(third layer + sixth layer)/2

These are examples but my data are more than this so I prefer a loop to do the job no matter how many layers I have in x

  • In a situation where you had 8 input layers, would the 4th x-layer be the average of 4 and 8 or 4 and 7? In case of 4 and 7, would the output be comprised of 5 layers? – Mikkel Lydholm Rasmussen Apr 15 '15 at 13:53
  • My real data are like this:three years (2007-2008-2009) so I stack all these datasets then compute the moving average.lastly I do this day1(2007)+day1(2008)+day1(2009)/3 it and the same for day 2 day2(2007)+day2(2008)+day2(2009)/3. This is what I tried to explain in my question. In my rasters I will know the number of the layer when 2008 or 2009 starts. So It depends at which layer ,the year 2008 starts.if 8 input layers, then 1+5,2+6,3+7,4+8 – usersam Apr 15 '15 at 14:12
  • By the way this is what they call it climatology in climate fields. – usersam Apr 15 '15 at 14:20
  • You are attempting, unsuccessfully, to apply a rolling average. Using this approach, follow your original attempt with movingFUN. Honestly, you would be better served using something like a polynomial regression to smooth the data. The mean is very sensitive to skewness in the data, which is quite expected. Whereas, an approach like lowess regression is robust to stochasticity in the climate series. – Jeffrey Evans Apr 15 '15 at 17:34
  • Why do you ask the same question at the same time on two different forums (here, and stackoverflow)??? I answered on stackoverflow: stackoverflow.com/questions/29651155/… – Robert Hijmans Apr 16 '15 at 5:28
2

It is fairly straight forward to set up the looping logic with an i,j index. However, I do not quite get your logic. What happens after (n - 4)? You can only calculate the adjusted mean to day 361. Why does calc or overlay, with movingFun, not work for you?

That aside, addressing your question, the missing piece is that you index rasters in a stack using double brackets; eg., st1[1] <- (st1[1] + st1[[4]])/2

# Using your example, which is somewhat questionable:
library(raster)
r1 <- r2 <- r3 <- r4 <- r5 <- r6 <- raster(nrows=10, ncols=10);
  r1[] <- runif(ncell(r1),1,100)
  r2[] <- runif(ncell(r2),1,100)
  r3[] <- runif(ncell(r3),1,100)
  r4[] <- runif(ncell(r4),1,100)
  r5[] <- runif(ncell(r5),1,100)
  r6[] <- runif(ncell(r6),1,100)
  r <- stack(r1,r2,r3,r4,r5,r6)

# use an offset of 1 (r1+r2)/2, (r2+r3)/2, ..., making the output n-1
  r.mean <- stack(r1)
  for( i in 1:(nlayers(r)-1)) {
    j=i+1
    r.mean[[i]] <- (r[[i]] + r[[j]])/2 
  }

From your description, you are trying to derive normals over a short time period. I have never seen a normalization or smoothing equation that uses the same day over multiple time periods. Where did you get this idea because it is not how climate normals or smoothing is defined? I would take a look at this description on why low-pass (mean) smoothing is not always a good idea.

Depending on the process you are modeling, smoothing will constrain the variance and could provide misleading inference. You really want to have a measure of the error so you can specify an error term associated with the smoothed data. I would encourage you to read up on climate normalization and smoothing before running forward with a bunch of analysis.

I believe that NOAA now uses a harmonic mean approach to normals following;

Arguez, A., and S. Applequist (2013) A Harmonic Approach for Calculating Daily Temperature Normals Constrained by Homogenized Monthly Temperature Normals. Journal of Atmospheric and Oceanic Technology, 30:1259–1265.

Here is a general FAQ on calculating normals with the associated missing data rules commonly applied

  • I was trying to follow this For each data set, the seasonal climatology was computed as the 31 day moving average, with the moving averages based on data from all years for the 31 day period surrounding each day of year. in this article sciencedirect.com/science/article/pii/S003442571300206X – usersam Apr 16 '15 at 7:15
  • and the same was applied here:sciencedirect.com/science/article/pii/S0034425714004490 – usersam Apr 16 '15 at 7:19
  • This is completely different from what you describe as: "I do this day1(2007)+day1(2008)+day1(2009)/3". This is not a 30 day moving average! The paper you cite is actually using the harmonic approach I cited with a +/- 31 day window for the smoothing period. You can accomplish a simplified form of this analysis using movingFun with type="around" and circular=TRUE. – Jeffrey Evans Apr 16 '15 at 17:56
1

I have made a rough code example, that uses simple structures and attempts to make the script easy to understand and follow. It is likely that it is inefficient and could be structured much better.

The fact that you have 3 years is the key that we need to consider here (I will completely disregard the potential of a leap-year). Adding additional years is reasonably easy, but not "supported" below.

#assume that x already exists
nyear <- 3
ndays <- 365
averageyear <- stack() #empty raster stack that will be filled during the loop.
for(i in 1:ndays){
  acrossyears <- mean(x[i]+x[i+ndays]+x[i+2*ndays],na.rm=T) #calculate the average between 1, 366 and 731 and so forth
  tempfilename <- paste("PlaceholderString_",i,".tif",sep="")
  writeRaster(acrossyears,filename=tempfilename,format="GTiff",overwrite=T)
}

The above will result in a raster-stack with 365 layers, one for each day. Adding more years is done by simply adding more bits to the calculation of acrossyears. Dealing with leapyears could be done manually (adding 1 to [i+ndays] here and there, or just dropping the leap-day through subsetting), or by rewriting x into separate files for each year and the doing away with the use of variations of [i+ndays].

  • Thanks but I think it would easier(in terms of memory) to just write the results of acrossyearsright a way and no need to stack them again.so the final results acrossyears will be 365 layers(is this what is expected from your code?) – usersam Apr 15 '15 at 15:38
  • Last thing what if the number of days is not the same in each year? I am using remote sensing data and there some missing days within a year. Thanks in advance. – usersam Apr 15 '15 at 16:03
  • Added "support" for writing separate rasters for each timestep. As for dealing with missing days, you have a couple of options: insert no-data rasters instead to keep the number of observations consistent, or gapfill using software like TimeSAT. If these options doesn't work for you, the next step would be to make the script "time-aware" and using a raster time-series construct, but this is up to you. – Mikkel Lydholm Rasmussen Apr 15 '15 at 16:43
  • 1
    The way you are defining an empty stack will throw and extent error. You need to define the dimensions of the stack, ideally with a raster used in the analysis, and then populate it. You also do not want to keep "re-stacking" it but rather use "addLayer" to append the existing stack. Alternately, you could also populate the stack with the rasters from the original stack (n - step) and use an index to replace each raster (see my example). – Jeffrey Evans Apr 15 '15 at 17:24
  • @JeffreyEvans No error appears to occur when I make an empty raster stack and then fill it with rasters. That part aside, the scope of the question changed into writing each raster into a new file. – Mikkel Lydholm Rasmussen Apr 15 '15 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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