Как интерполировать значения шума во временные ряды изображений местности, используя данные флага cfmask?

Question

Мне нужно создать функцию для интерполяции значений шума во временных рядах изображений LandSat по среднему значению за промежуток времени 2 (xt+1 и xt-1).

Я использую продукт fmask для обнаружения облаков и теней, затем применяется интерполяция.

Для одного временного ряда:

Поскольку c2 - вектор временного ряда fmask (2 для облака и 4 для тени), а t2 - вектор временного ряда e vi:

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

Но это невозможно сделать с помощью функции 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"))}  

Answer 1

Я верю, что следующее удовлетворит ваши потребности. Я предполагаю, что:

1. Каждое изображение во времени является растром. Они помещены в RasterStack названный t2.stack заказал во все большее время. Изображение на iна период времени в стеке ссылается t2.stack[[i]],
2. Для каждого изображения во времени есть маска. Это позволяет использовать разные маски для разных времен (хотя они также могут быть одинаковыми). Они помещены в RasterStack названный c2.stack, Они соответственно упорядочены во времени как t2.stack, Fmask на iна период времени в стеке ссылается c2.stack[[i]],
3. Все растры (и изображения, и маски) имеют одинаковый размер (одинаковое количество строк, столбцов и, следовательно, пикселей).

Сначала мы смоделируем некоторые данные для иллюстрации:

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.mask является FALSE в ноль.
3. Наоборот, mask растр t2.stack[[i]] установить те пиксели, для которых cloud.shadow.mask является TRUE в ноль.
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 для эффективности. Здесь мы повторяем вычисления, используя оригинал t2.stack это было скопировано в 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)
}

Результаты одинаковы.

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