1

I am working with series of raster data.

I am checking temporal outlier of each pixel by using image before and after, then replace current pixel value (t) by average of image before (t-1) and after (t+1). I am able to run with test of 2 rows x 2 columns, but with the real data (NA value in background of image). It can be run but very slow. It cost more than one week to complete one image (979x526 pixels)

I am a beginner R user.

Following is the script

# rm(list=ls())
# setwd ("D:\\5Validation\\92NDVIValid\\GapTestNDVI")
# I've commented these out because these could catch people out

Files <- list.files(pattern="*.tif") 

library(raster)
library(stats)

#### Need to difine variables in loop
Ndvit       <- as.matrix(raster(Files[2]))                # The second image

Ndvit_1     <- as.matrix(raster(Files[1]))                 # The first image

Ndvit_3     <- as.matrix(raster(Files[3]))                 # The third image

Con1        <- (Ndvit-Ndvit_1) <(-0.01*Ndvit_1)           # Condition 1  
Con2        <- (Ndvit-Ndvit_3) <(-0.01*Ndvit_3)           # Condition 2

Con         <-  Con1 & Con2                               # Conbine both Condition 1&2 (T==1, F==0)
summary(Con1)

## Dimension of image in series
NRow       <- nrow(raster(Files[1]))

NCol       <- ncol(raster(Files[1]))

#Create blank Table to merge all output  
NdviReli  <- matrix(nrow = 979, ncol = 526)

Out       <- matrix(nrow = 979, ncol = 526)


## Need image before and after to fill, 
## So it will process from second image to lenghth of image - 1)

for (t in 2:(length(Files)-1)) {   
  for (i in 1:NRow) { 
    for (j in 1:NCol) { 
        x <- t-1         ## Set variabe t-1
        y <- t+1         ## ## Set variabe t+1

        ## Assigned input image as matrix
        Ndvit       <- as.matrix(raster(Files[t]))
        Ndvit_x     <- as.matrix(raster(Files[x]))
        Ndvit_y     <- as.matrix(raster(Files[y]))

        ## Condition 
        Con1     <- (Ndvit-Ndvit_x) <(-0.01*Ndvit_x) 
        Con2     <- (Ndvit-Ndvit_y) <(-0.01*Ndvit_y)  
        Con      <- Con1 & Con2
        if(is.na(Con[i,j])) {

           NdviReli[i,j] <- NA

          } else if  (Con[i,j] == 1) {
          ## If Condition satify
          NdviReli[i,j] <- (Ndvit_x[i,j] + Ndvit_y[i,j])/2

          } else  {
          ## If Condition satify
          NdviReli[i,j] <- Ndvit[i,j]
          }

          Reli <- print(paste(NdviReli[i,j])) ## Using this if want to check the value one by one
          Out <- raster(NdviReli)                                       ## Assign matrix as raster
          Ex       <- extent(c(487124.6, 618624.6, 1902496, 2147246))   ## Assign coordination, projection base on image
          extent(Out) <- Ex    
          projection(Out) <- CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
          writeRaster(Out,filename= paste("Reli",Files[t]),format="GTiff", overwrite=T)
    } 
  }
}
6
  • It's the three nested for loops. Can you provide a minimal reproducible example (i.e. an example image) that we can test alternatives with?
    – Phil
    Commented Sep 23, 2016 at 9:37
  • @ Midavalo: Thank you very much. Please, help me on this issue, I am new R user. I am sharing link of sample data as below. drive.google.com/drive/folders/…
    – Bang Tam
    Commented Sep 24, 2016 at 14:55
  • The Google Drive link has restricted permissions; I can't download any files.
    – Phil
    Commented Sep 26, 2016 at 20:50
  • @ Midavalo: Could you please to check the the link as below 4shared.com/folder/O6XyaGoE/TestingData.html
    – Bang Tam
    Commented Sep 28, 2016 at 6:59
  • It says I need to sign up for an account to download the files, which I'm not going to do. You could just change the permissions on your Google Drive version to public.
    – Phil
    Commented Sep 28, 2016 at 8:59

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.