title: Updating raster values in foreach loops in R? (raster time-series NA imputation)
body: <p>I am trying to impute NA values in a time series raster. Here is an reproducible example of my data:</p>
  
  <pre><code>library(raster)
  library(rgdal)
  library(doParallel)
  library(foreach)
  
  r1 &lt;- r2 &lt;- r3 &lt;- r4 &lt;- r5 &lt;- raster(nrow=100, ncol=100)
  values(r1) &lt;- runif(ncell(r1))
  values(r2) &lt;- runif(ncell(r2))
  values(r3) &lt;- runif(ncell(r3))
  values(r4) &lt;- runif(ncell(r4))
  values(r5) &lt;- runif(ncell(r5))
  
  s &lt;- stack(r1, r2, r3, r4, r5)
  time_series &lt;- brick(s)
  time_series[1, 30][2] &lt;- NA
  time_series[3, 20][3] &lt;- NA
  time_series[5, 10][5] &lt;- NA
  time_series[8, 40][4] &lt;- NA
  
  </code></pre>
  
  <p>There are packages such as gapfill but I found them too slow for my task. I found another method here as in the answer:
  <a href="https://gis.stackexchange.com/questions/279354/ndvi-time-series-with-missing-values">https://gis.stackexchange.com/questions/279354/ndvi-time-series-with-missing-values</a></p>
  
  <p>by: <a href="https://gis.stackexchange.com/users/8520/jeffrey-evans">https://gis.stackexchange.com/users/8520/jeffrey-evans</a></p>
  
  <p>I want to convert the for loop to foreach so I can calculate it for larger images. Here is the code with for loops:</p>
  
  <pre><code>impute.loess &lt;- function(y, x.length = NULL, s = 0.80, 
                           smooth.data = FALSE, ...) {
    if(is.null(x.length)) { x.length = length(y) }
    options(warn = -1)
    x &lt;- 1:x.length
    if (all(is.na(y))) {
      return(y)
    } else {
      p &lt;- loess(y ~ x, span = s, data.frame(x = x, y = y))
      if(smooth.data == TRUE) {
        y &lt;- predict(p, x)
      } else {
        na.idx &lt;- which( is.na(y) )
        if( length(na.idx) &gt; 1 ) {
          y[na.idx] &lt;- predict(p, data.frame(x = na.idx))
        }
      }   
      return(y)
    }
  }
  
  time_series_new &lt;- time_series
  
  time_series_new[] &lt;- NA
  
  for (rl in 1:nrow(time_series)) {
    v &lt;- getValues(time_series, rl, 1)
    time_series_new[rl,] &lt;- as.matrix(t(apply(v, MARGIN=1, FUN=impute.loess)))
  }
  
  </code></pre>
  
  <p>Foreach alternative I tried is this:</p>
  
  <pre><code>time_series_new2 &lt;- time_series
  
  time_series_new2[] &lt;- NA
  
  cl &lt;- parallel::makeCluster(detectCores()-1)
  doParallel::registerDoParallel(cl)
  
  time_series_new2 &lt;- foreach (rl = 1:nrow(time_series),
                               .packages = "raster",
                               .combine = 'rbind') %dopar% {
                          v &lt;- getValues(time_series, rl, 1)
                          time_series_new[rl,] &lt;- as.matrix(t(apply(v,
                                                  MARGIN=1, FUN=impute.loess)))
  }
  
  parallel::stopCluster(cl)
  
  </code></pre>
  
  <p>However, here is the difference:</p>
  
  <pre><code>&gt; class(time_series_new)
  [1] "RasterBrick"
  attr(,"package")
  [1] "raster"
  
  &gt; class(time_series_new2)
  [1] "matrix"
  
  </code></pre>
  
  <p>If I don't assign foreach loop to an object, it just export the outcome. I want an updated raster object at the end but couldn't find a solution to my problem.</p>
  
  <p>I couldn't find how to set matrix values raster object - setvalues didn't work maybe since the dimensions are different as:</p>
  
  <pre><code>&gt; dim(time_series_new)
  [1] 100 100   5
  
  &gt; dim(time_series_new2)
  [1] 10000     5
  
  </code></pre>
  
  <p>I know foreach loop works differently. Is there a way to update time_series_new2 object inside foreach loop so that I can have an updated raster object at the end?</p>
  
  <p>EDIT:</p>
  
  <p>setValues() actually works! as:</p>
  
  <pre><code>time_series_new3 &lt;- time_series
  
  time_series_new3[] &lt;- NA #empty raster object
  
  time_series_new3 &lt;- setValues(time_series_new3, time_series_new2) #filled with matrix rendered from foreach loop
  
  &gt; time_series_new3
  class      : RasterBrick 
  dimensions : 100, 100, 10000, 5  (nrow, ncol, ncell, nlayers)
  resolution : 3.6, 1.8  (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 
  min values : 1.468023e-04, 3.525158e-04, 9.689084e-05, 5.349121e-05, 4.214607e-05 
  max values :    0.9999564,    0.9999854,    0.9997795,    0.9999780,    0.9997880 
  
  &gt; time_series_new2
  class      : RasterBrick 
  dimensions : 100, 100, 10000, 5  (nrow, ncol, ncell, nlayers)
  resolution : 3.6, 1.8  (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 
  min values : 1.468023e-04, 3.525158e-04, 9.689084e-05, 5.349121e-05, 4.214607e-05 
  max values :    0.9999564,    0.9999854,    0.9997795,    0.9999780,    0.9997880
  
  &gt; all.equal(time_series_new2, time_series_new3)
  [1] TRUE
  
  </code></pre>
  
  <p>Still, I would like to know about updating in foreach.</p>