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"))}
I believe the following will meet your needs. I'm assuming:
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]]. 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]].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:
cloud.shadow.mask is a raster that is TRUE when there is cloud or shadow at the pixel and false otherwise.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.mask the raster t2.stack[[i]] to set those pixels for which cloud.shadow.mask is TRUE to zero.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