This function fills NA with interpolation. Only NA locations with valid up, down, left, and right values will be interpolated.

fill-raster-NA.R

# "`-''-/").___..--''"`-._
#  (`6_ 6  )   `-.  (     ).`-.__.`)   WE ARE ...
#  (_Y_.)'  ._   )  `._ `. ``-..-'    PENN STATE!
#    _ ..`--'_..-_/  /--'_.' ,'
#  (il),-''  (li),'  ((!.-'
# 
# Author: Weiming Hu <wuh20@psu.edu>
#         Geoinformatics and Earth Observation Laboratory (http://geolab.psu.edu)
#         Department of Geography and Institute for CyberScience
#         The Pennsylvania State University

# This function fills NA vlaues with interpolation.
# Only NA locations with valid values from the up,
# down, left, and right locations will be interpolated.
#
fill.raster.NA <- function(rast.in, verbose = F) {
  require(raster)
  rast <- rast.in
  
  cells.na <- which(is.na(values(rast)))
  
  if (length(cells.na) != 0) {
    # only when both the grids above and below 
    # are valid values, we interpolate the middle
    #
    if (verbose) {
        print("Extracting NA grid points with valid up/down grid points ...")
    }
    rows <- rowFromCell(rast, cells.na)
    cols <- colFromCell(rast, cells.na)
    
    # ignore the first and the last rows
    index <- which(rows > 1)
    rows <- rows[index]
    cols <- cols[index]
    cells.na <- cells.na[index]
    index <- which(rows < nrow(rast))
    rows <- rows[index]
    cols <- cols[index]
    cells.na <- cells.na[index]
    
    values.up <- rast[matrix(c(rows - 1, cols), ncol = 2)]
    values.down <- rast[matrix(c(rows + 1, cols), ncol = 2)]
    
    # grid points with valid values in up and down grids
    index <- intersect(
      which(!is.na(values.up)),
      which(!is.na(values.down)))
    if (verbose) {
        print(paste(length(index), " gird points to fill.", sep = ''))
    }
    
    # compute and assign the mean
    if (verbose) {
        print("Assign grid points ...")
    }
    
    rast[cells.na[index]] <- 
      (values.up[index] + values.down[index]) / 2
    
    # only when both the grids left and right
    # are valid values, we interpolate the middle
    #
    if (verbose) {
        print("Extracting NA grid points with valid left/right grid points ...")
    }
    
    cells.na <- which(is.na(values(rast)))
    rows <- rowFromCell(rast, cells.na)
    cols <- colFromCell(rast, cells.na)
    
    # ignore the first and the last cols
    index <- which(cols > 1)
    rows <- rows[index]
    cols <- cols[index]
    cells.na <- cells.na[index]
    index <- which(cols < ncol(rast))
    rows <- rows[index]
    cols <- cols[index]
    cells.na <- cells.na[index]
    
    values.left <- rast[matrix(c(rows, cols - 1), ncol = 2)]
    values.right <- rast[matrix(c(rows, cols + 1), ncol = 2)]
    
    # grid points with valid values in up and down grids
    index <- intersect(
      which(!is.na(values.left)),
      which(!is.na(values.right)))
    if (verbose) {
        print(paste(length(index), " gird points to fill.", sep = ''))
    }
    
    # compute and assign the mean
    if (verbose) {
        print("Assign grid points ...")
    }
    rast[cells.na[index]] <- 
      (values.left[index] + values.right[index]) / 2
  }
  
  if (verbose) {
      print("Done!")
  }
  return(rast)
}