R/landsat_transform.R

In paleolimbot/landsatutils: Implements Helpers for Common Landsat Workflows

#' Propogate Landsat attributes through transformations
#'
#' This function applies a raster transformation function (a function that takes
#' a RasterStack/RasterBrick and returns a RasterStack/RasterBrick) to a
#' \link{landsat_scene} or landsat_scene_df (from \link{landsat_load_scenes}).
#'
#' @param x A \link{landsat_scene} or landsat_scene_df (from \link{landsat_load_scenes})
#' @param fun A funcion takes a RasterStack/RasterBrick and returns a
#'   RasterStack/RasterBrick
#' @param ... Arguments to be passed to \code{fun}
#'
#' @return A  \link{landsat_scene}, RasterStack, RasterBrick, or landsat_scene_df
#' @export
#'
landsat_transform <- function(x, ...) UseMethod("landsat_stransform")

#' @rdname landsat_transform
#' @export
landsat_transform.default <- function(x, fun, ...) {
  # this is a base method that applies a function to a landsat_scene and propogates
  # the landsat_attrs

  # make sure input is a landsat scene
  if(!is.landsat_scene(x)) stop("Use landsat_scene() to create a landsat scene")

  # apply method
  transformed <- fun(x, ...)

  # propogate attributes
  attr(transformed, "landsat_attrs") <- attr(x, "landsat_attrs")

  # return transformed version
  transformed
}

#' @rdname landsat_transform
#' @export
landsat_transform.landsat_scene_list <- function(x, fun, ...) {
  structure(
    lapply(x, function(scene) landsat_transform.default(scene, fun, ...)),
    class = "landsat_scene_list"
  )
}

#' @rdname landsat_transform
#' @export
landsat_transform.landsat_scene_df <- function(x, fun, ...) {
  x$scene <- landsat_transform(x$scene, fun, ...)
  x
}

#' @rdname landsat_transform
#' @export
landsat_transform.data.frame <- function(x, fun, ...) {
  if("scene" %in% colnames(x)) {
    x$scene <- landsat_transform(x$scene, fun, ...)
    class(x) <- c("landsat_scene_df", class(x))
    x
  } else {
    stop("x must have a scene column to be passed to landsat_transform")
  }
}


#' Crop Landsat objects to a boundary
#'
#' This function crops a \link{landsat_scene} or landsat_scene_df
#' (as generated by \link{landsat_load_scenes}) to a spatial object
#' describing a boundary. This is usually an object of type
#' \link[sp]{SpatialPolygons-class}, but could also be from the
#' sf package, a \link[sp]{bbox}, or an \link[raster]{extent}.
#'
#' @param x A \link{landsat_scene} or landsat_scene_df (from \link{landsat_load_scenes})
#' @param boundary Usually an object of type
#'   \link[sp]{SpatialPolygons-class}, but could also be from the
#'   sf package, a \link[sp]{bbox}, or an \link[raster]{extent}.
#'
#' @return A landsat_scene or landsat_scene_df
#' @export
#'
landsat_crop <- function(x, boundary) {
  # sanitize boundary
  boundary_proj <- sanitize_boundary(x, boundary)

  # use landsat_transform to crop and propogate attributres
  landsat_transform(x, raster::crop, boundary_proj)
}

#' @rdname landsat_crop
#' @export
landsat_mask <- function(x, boundary) {
  # sanitize boundary
  boundary_proj <- sanitize_boundary(x, boundary)

  # use landsat_transform to mask and propogate attributres
  landsat_transform(x, raster::mask, boundary_proj)
}


#' Project Landsat objects
#'
#' This function projects a \link{landsat_scene} or landsat_scene_df
#' (as generated by \link{landsat_load_scenes}) using an object
#' describing a CRS This is usually an object of type
#' \link[sp]{CRS}, but could also be from the
#' sf package (\link[sf]{st_bbox}), an integer describing an EPSG code,
#' or a character string describing proj4 arguments.
#'
#' @param x A \link{landsat_scene} or landsat_scene_df (from \link{landsat_load_scenes})
#' @param crs_obj Usually an object of type
#'   \link[sp]{CRS}, but could also be from the
#'   sf package (\link[sf]{st_bbox}), an integer describing an EPSG code,
#'   or a character string describing proj4 arguments.
#'
#' @return A landsat_scene or landsat_scene_df
#' @export
#'
landsat_project <- function(x, crs_obj) {
  # use landsat_transform to call projectRaster and propogate attributes
  landsat_transform(x, raster::projectRaster, sanitize_CRS(crs_obj))
}

#' Overlay a function on a Landsat object
#'
#' This function makes it easy to calculate indicies such as NDWI, NDVI,
#' etc. on \link{landsat_scene} objects.
#'
#' @param x A \link{landsat_scene} or landsat_scene_df (from \link{landsat_load_scenes})
#' @param fun Function with formal arguments B1, B2, B3, etc.
#' @param ... Passed to raster::\link[raster]{overlay}
#'
#' @return A landsat_scene or landsat_scene_df
#' @export
#'
landsat_overlay <- function(x, fun, ...) {
  # use landsat_transform to handle multiple classes and
  transformed <- landsat_transform(x, landsat_overlay_base, fun, ...)
  # custom landsat_scene classes aren't really applicable anymore
  attr(transformed, "landsat_attrs") <- NULL
  class(transformed) <- class(transformed)[!grepl("landsat", class(transformed))]
  transformed
}


landsat_overlay_base <- function(x, fun, ...) {
  # ensure scene is a landsat_scene
  if(!is.landsat_scene(x)) stop("Cannot use landsat_overlay without a landsat_scene")

  # ensure fun is a function
  fun <- match.fun(fun)

  # extract arguments of the function
  fun_args <- names(formals(fun))

  # extract bands available from landsat_attrs
  band_names <- attr(x, "landsat_attrs")$.band_names

  # check that all band_names are in fun_ags
  missing_args <- fun_args[!(fun_args %in% band_names)]
  if(any(missing_args)) stop("The following bands are missing in scene that are required by fun: ",
                             paste(missing_args, collapse = ", "))

  # modify function arguments to include all the band names so that raster::calc can be used
  formals(fun) <- stats::setNames(rep(list(rlang::missing_arg()), length(band_names)), band_names)

  # return result of raster::overlay
  raster::overlay(x, fun = fun, ...)
}



# private method to sanitize various objects that get passed as the boundary
sanitize_boundary <- function(x, boundary) {
  # make sure boundary is an sp object
  if(inherits(boundary, "sf") || inherits(boundary, "sfc")) {
    # sf objects can be coerced to sp objects
    boundary <- methods::as(boundary, "Spatial")
  } else if(inherits(boundary, "matrix") &&
            (nrow(boundary) == 2) && (ncol(boundary) == 2)) {
    # bbox objects get turned into polygons with same CRS as x
    boundary <- bbox_as_spatial(
      xmin = box[1, 1],
      xmax = box[1, 2],
      ymin = box[2, 1],
      ymax = box[2, 2],
      n = 10,
      crs = x@crs
    )
  } else if(inherits(boundary, "bbox")) {
    boundary <- bbox_as_spatial(
      xmin = boundary["xmin"],
      xmax = boundary["xmax"],
      ymin = boundary["ymin"],
      ymax = boundary["ymax"],
      n = 10,
      crs = x@crs
    )
  } else if(methods::is(boundary, "Extent")) {
    boundary <- bbox_as_spatial(
      xmin = boundary@xmin,
      xmax = boundary@xmax,
      ymin = boundary@ymin,
      ymax = boundary@ymax,
      n = 10,
      crs = x@crs
    )
  }

  # make sure boundary is in the same crs as the scene object
  boundary_proj <- sp::spTransform(boundary, x@crs)

  # return projected boundary
  boundary_proj
}

# private method to turn bounding boxes into spatial objects
bbox_as_spatial <- function(xmin, xmax, ymin, ymax, n = 10, crs = sp::CRS(NA_character_)) {
  # generate a polygon with the corners, and the dots filled in
  coords <- unique(rbind(
    data.frame(x = xmin, y = seq(ymin, ymax, length.out = n)),
    data.frame(x = seq(xmin, xmax, length.out = n), y = ymax),
    data.frame(x = xmax, y = seq(ymax, ymin, length.out = n)),
    data.frame(x = seq(xmax, xmin, length.out = n), y = ymin)
  ))

  # convert to spatial polygons object
  sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(coords)), 1)), proj4string = crs)
}

# private method to sanitize CRS objects
sanitize_CRS <- function(x) UseMethod("sanitize_CRS")
sanitize_CRS.CRS <- function(x) x
sanitize_CRS.numeric <- function(x) sp::CRS(sprintf("+init=epsg:%d", x))
sanitize_CRS.character <- function(x) sp::CRS(x)
sanitize_CRS.crs <- function(x) sp:CRS(x$proj4string)