1

Suppose I have 4 raster layers each belong to every other week of the month. I want to use linear interpolation to create new layers for each day. In this case, the first 2 rasters belonging to the month of Feb with 29 days and the second 2 belong to March with 31 days. I want to know how to create daily raster objects that can fill the time period with the consideration of number of days in the month (29 rasters for Feb and 31 rasters for March). Thanks!

library(raster)
r1 <- raster(nrow=5, ncol=5)
feb15 <- setValues(r1, runif(ncell(r1)))
feb29 <- setValues(r1, runif(ncell(r1)))
mar15 <- setValues(r1, runif(ncell(r1)))
mar30 <- setValues(r1, runif(ncell(r1)))
Geo-sp
  • 1,704
  • 3
  • 19
  • 42
  • is this question a duplicate of http://stackoverflow.com/questions/18339400/r-interpolation-between-raster-layers-of-different-dates?rq=1 ? – Jealie Oct 26 '15 at 17:51

1 Answers1

1

About the same question was asked two weeks ago. Here is the adapted answer for your example. 

library(raster)
r1 <- raster(nrow=5, ncol=5)
feb15 <- setValues(r1, runif(ncell(r1)))
feb29 <- setValues(r1, runif(ncell(r1)))
mar15 <- setValues(r1, runif(ncell(r1)))
mar30 <- setValues(r1, runif(ncell(r1)))

s <- stack(feb15, feb29, mar15, mar30)
x <- calc(s, fun=function(y) approx(c(15,29,29+15,29+30), y, 1:59, rule=2)$y)
names(x) <- c(paste('feb', 1:29), paste('march', 1:30))

Unfortunately, approx fails when it only gets NA values. It "need at least two non-NA values to interpolate". Here is a work around:

s[1:2] <- NA
# fail <- calc(s, fun=function(y) approx(c(15,29,29+15,29+30), y, 1:59, rule=2)$y)

f <- function(y) {
    if (sum(! is.na(y)) < 2) {
        return(rep(NA, 59))
    } else {
        approx(c(15,29,29+15,29+30), y, 1:59, rule=2)$y
    }
}

x <- calc(s, fun=f)
Community
  • 1
  • 1
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thanks @RobertH ! Would you tell me what that 29+15, 29+30 does? I want to expand this to cover the whole year. – Geo-sp Oct 27 '15 at 14:14
  • in this case these are day numbers that start at Feb 1. To do the whole year, I would start counting at Jan 1. – Robert Hijmans Oct 27 '15 at 15:49
  • Thanks! I wrote it this way:`x <- calc(s, fun=function(y) approx(c(15,31,46,59,74,90, 105,120,135,151,166,181,196,212,227,243,258,273,288,304,319,334,349,365), y, 1:365, rule=2)$y)` but produce the following error: `Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function` – Geo-sp Oct 27 '15 at 16:08
  • This is probably because of cells with only NA values. I have added a work-around above. – Robert Hijmans Oct 27 '15 at 16:49