Wie interpoliert man Rauschwerte in einer Landsat-Bild-Zeitreihe mit cfmask-Flag-Daten?

Question

Ich muss eine Funktion erstellen, um Rauschwerte in einer Landsat-Bild-Zeitreihe zu interpolieren, und zwar in einer Zeitspanne von 2 (xt + 1 und xt-1).

Ich verwende das fmask-Produkt, um Wolken und Schatten zu erkennen, dann wird die Interpolation angewendet.

Für eine Zeitreihe:

Da c2 der Vektor der fmask Zeitreihe (2 für Cloud und 4 für Schatten) und t2 der Vektor der evi Zeitreihe:

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

Aber es Dies ist mit der Berechnungsfunktion des Raster-Pakets nicht möglich, da es nicht mit Funktionen mit 2 Parametern funktioniert.

Irgendwelche Vorschläge, wie Sie damit umgehen und diese Interpolation für alle Pixel der Raster-Zeitreihe anwenden?

Ich versuche dies, aber es funktioniert immer noch nicht:

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

Ich glaube, die nach Ihren Bedürfnissen gerecht wird. Ich nehme an:

1. Jedes Bild in der Zeit ist ein Raster. Diese werden in einer RasterStack mit dem Namen t2.stack geordnet in zunehmender Zeit platziert. Das Bild mit dem i-ten Zeitraum im Stapel wird durch t2.stack[[i]] referenziert.
2. Für jedes Bild in der Zeit gibt es eine fmask. Dies ermöglicht verschiedene fmasks für verschiedene Zeiten (obwohl sie auch alle gleich sein können). Diese werden in eine RasterStack mit dem Namen c2.stack platziert. Diese sind zeitlich entsprechend als t2.stack geordnet. Die fmask bei der i-ten Zeitperiode in dem Stapel wird durch c2.stack[[i]] referenziert.
3. Alle Raster (sowohl Bilder als auch Fmasks) haben die gleiche Größe (gleiche Anzahl von Zeilen, Spalten und daher Pixel).

Wir simulieren zunächst einige Daten veranschaulichen:

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.23.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 

Der Schlüssel für die Verarbeitung ist Masken anstelle der bedingten Anweisung zu verwenden. Der Code lautet wie folgt:

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 
} 

Anmerkungen:

1. cloud.shadow.mask ein Raster ist, das TRUE ist, wenn es Wolke oder Schatten an dem Pixel ist, und sonst falsch.
2. Verwenden mask auf dem Raster, der die Mitte zwischen t2.stack[[i-1]] und t2.stack[[i+1]] ist jene Pixel zu setzen, für die cloud.shadow.maskFALSE Null ist.
3. Umgekehrt, mask das Raster t2.stack[[i]], um die Pixel, für die cloud.shadow.mask ist TRUE auf Null setzen.
4. Dann fügen Sie diese zwei Raster zu aktualisieren t2.stack[[i]].Hier

ist das Ergebnis für nT=3 wo nur t2.stack[[2]] berechnet:

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 

Für große Bilder, die nicht in den Speicher passen kann, verwenden overlay statt calc für Effizienz. Hier haben wir die Berechnung wiederholen Sie den ursprünglichen t2.stack verwenden, der 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) 
} 

Die Ergebnisse sind die gleichen kopiert wurde.

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