如何使用 cfmask 标志数据在陆地卫星图像时间序列上插入噪声值?

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

我需要创建一个函数,根据时间跨度 2(xt+1 和 xt-1)的平均值对 landsat 图像时间序列的噪声值进行插值。

我正在使用 fmask 产品检测云和阴影,然后应用插值。

对于一个时间序列:

因为c2是fmask时间序列的向量(2代表云,4代表阴影),t2是evi时间序列的向量:

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

但是无法使用 raster 包的 calc 函数来执行此操作,因为它不适用于具有 2 个参数的函数。

关于如何处理此问题并将此插值应用于栅格时间序列的所有像素的任何建议?

我正在尝试这个,但它仍然不起作用:

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"))}  

相信以下会满足您的需求。我假设:

  1. 时间上的每个图像都是一个光栅。这些被放置在一个名为 t2.stackRasterStack 中,按增加的时间排序。堆栈中第 i 个时间段的图像由 t2.stack[[i]] 引用。
  2. 对于每个时间图像,都有一个 fmask。这允许不同时间使用不同的 fmask(尽管它们也可以全部相同)。这些被放置在名为 c2.stackRasterStack 中。这些在时间上相应地排序为 t2.stack。堆栈中第 i 个时间段的 fmask 被 c2.stack[[i]].
  3. 引用
  4. 所有光栅(图像和 fmasks)都具有相同的大小(相同的行数、列数,因此像素数)。

我们先模拟一些数据来说明:

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

处理的关键是使用掩码代替条件语句。代码如下:

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
}

备注:

  1. cloud.shadow.mask 是一个光栅,当像素处有云或阴影时为 TRUE,否则为 false。
  2. 在栅格上使用 mask,即 t2.stack[[i-1]]t2.stack[[i+1]] 之间的平均值,将 cloud.shadow.maskFALSE 的那些像素设置为零。
  3. 相反,mask 栅格 t2.stack[[i]]cloud.shadow.maskTRUE 的那些像素设置为零。
  4. 然后将这两个栅格相加更新t2.stack[[i]].

这是 nT=3 的结果,其中只计算了 t2.stack[[2]]

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

对于可能无法放入内存的大图像,使用 overlay 而不是 calc 以提高效率。在这里,我们使用复制到 t3.stack

的原始 t2.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)
}

结果是一样的

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