R/create_points.R

In jnghiem/bfasttools: Tools for BFAST spatial analysis

#' Create an ESRI shapefile containing point features indexed to cell number
#'
#' This function converts a raster into a shapefile with point features
#' corresponding to each cell. Each point will have an attribute for its cell
#' number.
#'
#' If an output file is specified, this function will write an ESRI shapefile to
#' that file. The cell number will appear as the "no_cell" field in the
#' attribute table.
#'
#' @param template A Raster object to be convert into point features.
#'   Alternatively, an XML file generated by \code{create_raster_metadata()}.
#' @param boundary An optional SpatialPolygons object to subset the point
#'   features after they have been generated according to the the grid from
#'   \code{template}. It should be in the same projection as \code{template}.
#' @param output_file An optional path to a .shp file to which the shapefile
#'   will be written. If this argument is not supplied, the features will be
#'   returned as a SpatialPointsDataFrame object.
#' @return A SpatialPointsDataFrame object or nothing (with a side-effect of
#'   writing to a file).
#' @examples
#' \dontrun{
#' create_points(raster("C:/Desktop/extent.tif"), "C:/Desktop/shapefiles/all_points.shp")
#' }
#' @importFrom rgdal writeOGR
#' @importFrom rgeos gIntersection
#' @export
create_points <- function(template, boundary=NULL, output_file=NULL) {
  template <- create_template_raster(template)
  all_cells <- 1:ncell(template)
  df <- data.frame(no_cell=all_cells, index=all_cells)
  pts <- create_raster("no_cell", "index", df, template) %>%
    rasterToPoints(spatial=TRUE)
  if (!is.null(boundary)) {
    boundary@data <- data.frame()
    pts <- intersect(pts, boundary)
  }
  pts@data <- dplyr::rename(pts@data, no_cell=layer)
  if (is.null(output_file)) {
    return(pts)
  } else {
    writeOGR(pts, dsn=dirname(output_file), layer=tools::file_path_sans_ext(basename(output_file)), driver="ESRI Shapefile")
  }
}