grow

# GROW A RASTER SURFACE TOWARDS 'NA' AREAS USING THE VALUE OF THE CLOSEST PIXEL

# version 1.3 (2 Jan 2022)

# this function is based on a script provided by Lennert at https://stackoverflow.com/questions/27562076/if-raster-value-na-search-and-extract-the-nearest-non-na-pixel/37392542#37392542

# it just converts the script to a function, and to 'terra' (rather than 'raster') package functions wherever possible
# requires 'terra' package development version, installed after 2 Jan 2022 with:
# install.packages("terra", repos='https://rspatial.r-universe.dev')


grow <- function(rst,  # a SpatRaster or Raster layer in a projected CRS
                 dist = 10000  # distance up to which to grow pixels; if NA, the entire raster is filled with values
) {

  if (is.lonlat(rst)) stop("'rst' is in unprojected lon-lat coordinates; this function requires a projected 'rst'.")
  else if (is.lonlat(rst, perhaps = TRUE)) warning("'rst' does not have a defined CRS, but the coordinates look like they could be unprojected lon-lat. This function requires a projected 'rst'.")

  if(is(rst, "Raster"))  rst <- rast(rst)

  dis <- distance(rst)
  if (is.finite(dist))  dis[dis > dist] <- NA

  dir <- direction(rst, from = FALSE)

  rna <- is.na(rst)  # returns NA raster

  # store coordinates in new raster: https://stackoverflow.com/a/35592230/3752258
  na.x <- init(rna, "x")
  na.y <- init(rna, "y")

  # calculate coordinates of the nearest non-NA pixel
  # if we have an orthogonal, projected CRS, we can use (Pythagorean) calculations:
  co.x <- na.x + dis * sin(dir)
  co.y <- na.y + dis * cos(dir)

  # matrix with point coordinates of nearest non-NA pixel
  co <- cbind(co.x[], co.y[])

  # extract values of nearest non-NA cell with coordinates 'co'
  NAVals <- extract(rst, co)
  r.NAVals <- rna  # initiate new raster
  r.NAVals[] <- NAVals[ , 1]  # store values in raster

  # cover nearest non-NA value at NA locations of original raster
  r.filled <- cover(x = rst, y = r.NAVals)

  return(r.filled)
}