How to interpolate noise values on a landsat image time series using cfmask flag data?

Question

I need to create a function to interpolate noise values on a landsat image time series by the mean on a time span of 2 (xt+1 and xt-1).

I´m using the fmask product to detect cloud and shadow, then interpolation is applied.

For one time series:

Since c2 is the vector of fmask time series (2 for cloud and 4 for shadow), and t2 the vector of evi time series:

    for (i in 2:(length(t2)-1)){ 
         if (c2[i]==2 | c2[i]==4) 
         t2[i]<-mean(c(t2[i-1], t2[i+1]))}

But it is not possible to do this using the calc function of raster package, because it does not works with functions with 2 parameters.

Any suggestion about how to deal with this and apply this interpolation for all the pixels of the raster time series?

I´m trying this, but it is still not working:

for (i in 2:(length(stacklist)-1)){ 
re<-raster(stacklist[i]) 
re1<-raster(stacklist[i+1]) 
re0<-raster(stacklist[i-1]) 
rc<-raster(stacklist2[i]) 
if (rc[i]==2 | rc[i]==4) re[i]<-mean(c(re0[i],re1[i])) 
writeRaster(re,filename =paste0(substr(stacklist[i], 48, 59),"_filtered.tif"))}  

Answer 1

I believe the following will meet your needs. I'm assuming:

1. Each image in time is a raster. These are placed in a RasterStack named t2.stack ordered in increasing time. The image at the i-th time period in the stack is referenced by t2.stack[[i]].
2. For each image in time there is an fmask. This allows for different fmasks for different times (although they can also all be the same). These are placed in a RasterStack named c2.stack. These are correspondingly ordered in time as t2.stack. The fmask at the i-th time period in the stack is referenced by c2.stack[[i]].
3. All rasters (both images and fmasks) are of the same size (same number of rows, columns, and therefore, pixels).

We first simulate some data to illustrate:

library(raster)
set.seed(42)     ## This is for repeatable results

## Simulate some stacks of rasters over a time period [1,nT=3], yours will be longer
## 1. Each raster is 10 x 10, yours will be larger
## 2. Each c2 raster contains integers uniformly distributed in [0, 5]
## 3. Each t2 raster contains reals unformly distributed in [0,1]
## 4. c2.stack is a stack of c2 rasters over time period, c2.raster[[i]] is c2 at time i
## 5. t2.stack is a stack of t2 rasters over time period, t2.raster[[i]] is t2 at time i
nT <- 3
c2 <- raster(ncol=10, nrow=10)                  
c2[] <- floor(runif(ncell(c2), min=1, max=6))   
c2.stack <- stack(x=c2)
for (i in 2:nT) {
  c2[] <- floor(runif(ncell(c2), min=1, max=6))
  c2.stack <- addLayer(c2.stack, c2)
}

t2 <- raster(ncol=10, nrow=10)
t2[] <- runif(ncell(t2), min=0, max=1)
t2.stack <- stack(x=t2)
for (i in 2:nT) {
  t2[] <- runif(ncell(t2), min=0, max=1)
  t2.stack <- addLayer(t2.stack, t2)
}

## Here is the t2.stack data
print(head(t2.stack[[1]]))
##            1           2          3          4          5          6          7         8          9         10
##1  0.48376814 0.444569528 0.06038559 0.32750602 0.87842905 0.93060489 0.39217846 0.1588468 0.31994760 0.30696562
##2  0.10781125 0.979334303 0.49690343 0.09307467 0.21177366 0.93050075 0.29684641 0.6532182 0.90107048 0.99079579
##3  0.43033322 0.393776922 0.14190890 0.27980670 0.56482222 0.93513951 0.35840015 0.8420072 0.72240921 0.75073599
##4  0.92398845 0.002378107 0.16042991 0.39927295 0.67531958 0.48037202 0.53382878 0.3169502 0.81475759 0.29221952
##5  0.40913209 0.090918308 0.79859664 0.35978525 0.04048758 0.04108634 0.95443424 0.3733412 0.80641967 0.91005901
##6  0.44007621 0.576336503 0.07366780 0.16462739 0.73989078 0.47571101 0.68552095 0.9515149 0.49746449 0.47050063
##7  0.56019195 0.652510121 0.27957350 0.97990759 0.64386411 0.58257844 0.61587103 0.9251403 0.39002290 0.28791969
##8  0.09073596 0.322033904 0.75827011 0.10441293 0.71027785 0.96647738 0.20149123 0.1084887 0.05540218 0.82972352
##9  0.58119776 0.470092375 0.36501412 0.28012463 0.59971585 0.81856961 0.09783228 0.9636895 0.16873644 0.08608341
##10 0.86121070 0.524790602 0.65681088 0.22951937 0.72122603 0.49075039 0.96525559 0.9069425 0.55125053 0.07559910

print(head(t2.stack[[2]]))
##        1           2         3          4         5          6         7         8          9         10
##1  0.02270001 0.513239528 0.6307262 0.41877162 0.8792659 0.10798707 0.9802787 0.2649666 0.08427752 0.38590718
##2  0.12489583 0.581554222 0.2401496 0.72188789 0.1459287 0.15283877 0.2592227 0.7778863 0.42646630 0.06004834
##3  0.11483254 0.482756897 0.9791736 0.81151679 0.5429128 0.07236709 0.4664852 0.3399056 0.68991861 0.51415737
##4  0.51492302 0.545514354 0.4474573 0.08388484 0.9301337 0.01644819 0.4140924 0.2269761 0.09964059 0.48292388
##5  0.65012867 0.921329778 0.3626018 0.85513499 0.3009062 0.46566243 0.1427307 0.8077190 0.66580763 0.06194098
##6  0.43092557 0.396855081 0.6969568 0.65931965 0.4073507 0.30692022 0.2551073 0.6725682 0.89439343 0.84573616
##7  0.39290186 0.079050540 0.8284231 0.07289182 0.1147627 0.63998427 0.3205662 0.1887495 0.39382964 0.86202602
##8  0.34791141 0.001433898 0.9112845 0.95172345 0.4909190 0.46365172 0.5964720 0.9060510 0.17300118 0.78588108
##9  0.23293438 0.577048209 0.8408770 0.13220378 0.8958912 0.45013734 0.8941425 0.2485452 0.08369529 0.04864107
##10 0.97981587 0.484167741 0.8453930 0.41629360 0.4893425 0.18328782 0.7591615 0.3051433 0.16567825 0.03280914

print(head(t2.stack[[3]]))
##           1         2          3            4         5          6           7          8         9        10
##1  0.1365052 0.1771364 0.51956045 0.8111207851 0.1153620 0.89342179 0.575352881 0.14657239 0.9028058 0.2530025
##2  0.1505976 0.7685472 0.23012333 0.3053993280 0.5185696 0.33459967 0.154434968 0.26636957 0.3507546 0.5784584
##3  0.8086018 0.9332703 0.83386334 0.1270027745 0.6494540 0.69035166 0.032044824 0.92048915 0.4784689 0.2665206
##4  0.8565107 0.2291465 0.79194687 0.6467748603 0.4243347 0.09506827 0.003467704 0.53113367 0.5243071 0.2131856
##5  0.7169321 0.9613436 0.51826660 0.1745280223 0.5625401 0.75925817 0.666971338 0.22487292 0.3458498 0.3198318
##6  0.9048984 0.1991984 0.68096302 0.1375177586 0.1069947 0.09285940 0.916448955 0.27706044 0.8857939 0.7728646
##7  0.7950512 0.2056736 0.04819332 0.0388159312 0.2845741 0.34880983 0.737453325 0.25166358 0.5174370 0.7594447
##8  0.6360845 0.2039407 0.99304528 0.0004050434 0.2065700 0.63402809 0.017291825 0.02673547 0.6078406 0.5705414
##9  0.2458533 0.9195251 0.67221620 0.6454504393 0.2082855 0.48061774 0.986518583 0.99311985 0.4507740 0.7148488
##10 0.3165616 0.8336876 0.43397558 0.9959922582 0.8058112 0.48624217 0.538772001 0.34103991 0.0520156 0.4587231

## and the fmask product at time 2 (c2.stack[[2]])
print(head(c2.stack[[2]]))
##   1 2 3 4 5 6 7 8 9 10
##1  4 2 2 2 5 5 4 4 3  1
##2  4 5 4 3 3 3 1 2 4  5
##3  2 3 3 3 4 2 5 5 2  4
##4  5 4 4 5 5 3 5 1 4  4
##5  1 1 3 4 4 5 1 5 2  1
##6  4 2 4 2 4 4 1 1 1  4
##7  5 3 4 1 3 1 3 2 1  1
##8  4 3 3 3 3 1 5 3 4  4
##9  5 5 2 2 4 4 5 4 1  2
##10 1 4 1 1 1 1 3 1 4  4

## Make a copy of the t2.stack so that we can compare results using overlay later
t3.stack <- t2.stack

The key to the processing is to use masks in place of the conditional statement. The code is as follows:

for (i in 2:(nlayers(t2.stack)-1)) { 
  cloud.shadow.mask <- (c2.stack[[i]]==2 | c2.stack[[i]]==4)
  the.mean <- mask((t2.stack[[i-1]] + t2.stack[[i+1]])/2, cloud.shadow.mask, 
           maskvalue=FALSE, updatevalue=0.)
  the.middle <- mask(t2.stack[[i]], cloud.shadow.mask, 
         maskvalue=TRUE, updatevalue=0.)
  t2.stack[[i]] <- the.mean + the.middle
}

Notes:

1. cloud.shadow.mask is a raster that is TRUE when there is cloud or shadow at the pixel and false otherwise.
2. Use mask on the raster that is the mean between t2.stack[[i-1]] and t2.stack[[i+1]] to set those pixels for which cloud.shadow.mask is FALSE to zero.
3. Conversely, mask the raster t2.stack[[i]] to set those pixels for which cloud.shadow.mask is TRUE to zero.
4. Then add these two rasters to update t2.stack[[i]].

Here is the result for nT=3 where only t2.stack[[2]] is computed:

print(head(t2.stack[[2]]))
##           1           2         3          4         5          6         7         8          9         10
##1  0.3101367 0.310852970 0.2899730 0.56931340 0.8792659 0.10798707 0.4837657 0.1527096 0.08427752 0.38590718
##2  0.1292044 0.581554222 0.3635134 0.72188789 0.1459287 0.15283877 0.2592227 0.4597939 0.62591255 0.06004834
##3  0.6194675 0.482756897 0.9791736 0.81151679 0.6071381 0.81274558 0.4664852 0.3399056 0.60043905 0.50862828
##4  0.5149230 0.115762292 0.4761884 0.08388484 0.9301337 0.01644819 0.4140924 0.2269761 0.66953235 0.25270254
##5  0.6501287 0.921329778 0.3626018 0.26715664 0.3015139 0.46566243 0.1427307 0.8077190 0.57613471 0.06194098
##6  0.6724873 0.387767442 0.3773154 0.15107258 0.4234427 0.28428520 0.2551073 0.6725682 0.89439343 0.62168264
##7  0.3929019 0.079050540 0.1638834 0.07289182 0.1147627 0.63998427 0.3205662 0.5884019 0.39382964 0.86202602
##8  0.3634102 0.001433898 0.9112845 0.95172345 0.4909190 0.46365172 0.5964720 0.9060510 0.33162139 0.70013243
##9  0.2329344 0.577048209 0.5186152 0.46278753 0.4040007 0.64959367 0.8941425 0.9784047 0.08369529 0.40046610
##10 0.9798159 0.679239081 0.8453930 0.41629360 0.4893425 0.18328782 0.7591615 0.3051433 0.30163307 0.26716111

For large images that may not fit into memory, use overlay instead of calc for efficiency. Here we repeat the computation using the original t2.stack that was copied to t3.stack

for (i in 2:(nlayers(t3.stack)-1)) { 
  cloud.shadow.mask <- overlay(c2.stack[[i]], fun = function(x) x == 2 | x == 4)
  the.mean <- overlay(t3.stack[[i-1]], t3.stack[[i+1]], fun = function(x,y) (x+y)/2)
  the.mean <- mask(the.mean, cloud.shadow.mask, maskvalue=FALSE, updatevalue=0.)
  the.middle <- mask(t3.stack[[i]], cloud.shadow.mask, maskvalue=TRUE, updatevalue=0.)
  t3.stack[[i]] <- overlay(the.mean, the.middle, fun=sum)
}

The results are the same.

print(all.equal(t2.stack[[2]],t3.stack[[2]]))
##[1] TRUE