ApproxNA producing different and incorrect results in stocks of large rasters

I'm trying to confirm if there is a bug or at least an unexplained behavior change in the approxNA function in the latest version (2.5-8) of the R Raster package. Later, please note, the correct answer is seen as a comment in the answer line below, informing about the update of the package that was released

I have a small reproducible example below. The initial symptoms were that most of the approxNA command results on the annual daily eMODIS brick turned out to be the same as the input - it had many, many NA values ​​that were not interpolated ... but when I switched to an older version of the package ( 2.3-12) and older R version (3.3.1, Bug-in-your-hair), it returns correct interpolated results. I tried this using the raster datatype set to INT2S and no raster dataset defined.

I have checked the documentation in case there is some reason to expect the results to be different, but there is no change in the documentation between the old version of the package that worked correctly, and the new version of the package does not work correctly.

In my reproducible example, there are fewer errors, but they are obvious to view only from the graph (in the example). In these plots, you can see that the northern portions of several layers have NA values ​​that should have been interpolated.

I'm running Windows 7 64-bit and R 3.3.2 64-bit with 2.5-8 bitmap with R-Studio interface, but I tested also on R console with the same results and on another computer that has the same configuration and also tested in 32 bit R with the same results.

The trigger seems to be the presence of NA values ​​in the first or last layer, but errors happen in cells that do not have NA values ​​in the first or last layer: for example, the top-left cells are only missing in layer 2 in the original stack, so it must be very by simple interpolation for these cells from level 1 to level 3.

I followed the R instructions for submitting a bug report a few months ago, but I'm not sure if this contact information is correct, or if it might have ended up in the spam folder.

I would also be happy with answers suggesting direction towards a workaround, as long as it's not too slow and uses regular packages or packages that are regularly maintained.

I cannot tag them correctly, but below:

input brick schedule

good weekend brick schedule

bad exit brick schedule

Thanks for any suggestions.

##### start  Reproducible example

library(raster)
# make a large raster, it must be quite big to have a problem
r <- raster(ncols=8142, nrows=3285)
# make six rasters to stack, use different ranges to better see the interpolation results
r1 <- setValues(r, round(10 * runif(ncell(r)),0))
r2 <- setValues(r, NA)
r3 <- setValues(r, round(100 * runif(ncell(r)),0))
r4 <- setValues(r, round(1000 * runif(ncell(r)),0))
r5 <- setValues(r, round(50 * runif(ncell(r)),0))
r6 <- setValues(r, round(100 * runif(ncell(r)),0))

# insert some NA values
r1[100:600,] <- NA
r3[,1500:1950] <- NA
r5[,400:800] <- NA
r6[2500:2900,] <- NA

# insert some constant values to make it easier to see if interpolation is working
r4[300:500,] <- 750
r6[300:400,] <- 20
r6[400:500,] <- 100
# make a stack
s <- stack(r1,r2,r3,r4,r5,r6)

# see what the pre-interpolation stack looks like
#plot(s)
plot(s, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))

# interpolate, this takes a while
x2 <- approxNA(s, method = "linear", rule=2)

# see what the post interpolation looks like, there should be filled in values for all
# parts of all layers but instead it’s retaining many of the NA values
#plot(x2)
plot(x2, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))


## End reproducible example