R/cellnumbers.R

Defines functions line_cellnumbers pix as.owin.BasicRaster cellnumbers.sf cellnumbers.sfc cellnumbers.SpatialLines cellnumbers.default cellnumbers

Documented in cellnumbers cellnumbers.default cellnumbers.sf cellnumbers.sfc cellnumbers.SpatialLines

#' Extract cell numbers from a Raster object.
#'
#' Provide the 'cellnumbers' capability of [raster::extract] and friends
#' directly, returning a data frame of query-object identifiers 'object_' 
#' and the cell number. 
#' 
#' Raster data is inherently 2-dimensional, with a time or 'level' dimension
#' treated like a layers of these 2D forms. The 'raster' package cell number
#' is counted from 1 at the top-left, across the rows and down. This corresponds
#' the the standard "raster graphics" convention used by 'GDAL' and the 'sp' 
#' package, and many other implementations. Note that this is different to the
#' convention used by the [graphics::image] function. 
#' 
#' Currently this function only operates as if the input is a single layer objects, it's not clear if adding
#' an extra level of grouping for layers would be sensible.   
#' 
#' @param x Raster object
#' @param ... unused
#' @param query Spatial object or matrix of coordinates
#' @return a data frame (tibble) with columns
#' 
#' * `object_` - the object ID (what row is it from the spatial object)
#' * `cell_`   - the cell number of the raster
#' 
#' @export
#' @importFrom dplyr bind_rows
#' @importFrom raster cellFromPolygon cellFromLine cellFromXY projection
#' @examples
#' library(raster)
#' library(dplyr)
#' r <- raster(volcano) %>% aggregate(fact = 4)
#' cellnumbers(r, rasterToContour(r, level = 120))
#' library(dplyr)
#' 
#' cr <- cut(r,  pretty(values(r)))
#' 
#' suppressWarnings(tt <- cellnumbers(cr, polycano))
#' library(dplyr)
#' tt %>% mutate(v = extract(r, cell_)) %>% 
#' group_by(object_) %>% 
#' summarize(mean(v)) 
#' head(pretty(values(r)), -1)
#' @export
cellnumbers <- function(x, query, ...) {
  UseMethod("cellnumbers", object = query)
}
#' @name cellnumbers
#' @export
cellnumbers.default <- function(x, query, ...) {
  a <- NULL ## allow a basic check for not understanding "query"
  if (inherits(query, "sf")) {
    ## we need this for points, mpoints
    pth <- silicate::sc_path(query)
 
    tab <- tibble(object_ = rep(as.integer(ordered(pth$object_, unique(pth$object_))), pth$ncoords_), 
                  cell_ = raster::cellFromXY(x, as.matrix(silicate::sc_coord(query)[c("x_" , "y_")])))
    return(tab)
                  
 }
  if (is.na(projection(x)) || is.na(projection(query)) || projection(x) != projection(query)) {
    message(sprintf("projections not the same \n    x: %s\nquery: %s", projection(x), projection(query)), call. = FALSE)
  }
  if (inherits(query, "SpatialPolygons")) {
    message("cellnumbers is very slow for SpatialPolygons, consider conversion with 'sf::st_as_sf'")
    a <- cellFromPolygon(x, query)
  }
  if (is.matrix(query) | inherits(query, "SpatialPoints")) {
    if (is.matrix(query)) stopifnot(ncol(query) >= 2)
    a <- as.list(cellFromXY(x, query))
  }
  if (inherits(query, "SpatialMultiPoints")) {
    a <- lapply(query@coords, function(xymat) cellFromXY(x, xymat))
  }
  if (is.null(a)) stop(sprintf("no method for 'query' of type %s", class(query)[1L]))
  d <- dplyr::bind_rows(lapply(a, mat2d_f), .id = "object_")
  d[["object_"]] <- as.integer(d[["object_"]])
  if (ncol(d) == 2L) names(d) <- c("object_", "cell_")
  if (ncol(d) == 3L) names(d) <- c("object_", "cell_", "weight_")
  d
  
}
#' @name cellnumbers
#' @export
cellnumbers.SpatialLines <- function(x, query, ...) {
  line_cellnumbers(query, x)
}
#' @name cellnumbers
#' @export
cellnumbers.sfc <- function(x,  query, ...) {
  if (!"list" %in% class(query)) {
    class(query) <- c(class(query), "list")
  }

   sf1 <- tibble::tibble(geometry = query)
   cellnumbers(x, structure(sf1, sf_column = "geometry", agr = NULL, class = c("sf", "data.frame")))
}
#' @name cellnumbers
#' @export
cellnumbers.sf <- function(x, query, ...) {
  g <- query[[attr(query, "sf_column")]]
  if (inherits(g, "sfc_GEOMETRY")) stop("GEOMETRY and GEOMETRYCOLLECTION not supported")
  if (inherits(g, "sfc_LINESTRING") || inherits(g, "sfc_MULTILINESTRING")) {
    return(line_cellnumbers(query, x))
  }
  
  if (inherits(g, "sfc_POLYGON") || inherits(g, "sfc_MULTIPOLYGON")) {
    query$polygon <- seq_len(nrow(query))
    rast <- fasterize::fasterize(query, x, field = "polygon")
    v <- values(rast)
    ok <- !is.na(v)
    return(tibble::tibble(object_ = v[ok], cell_ = which(ok)))
  }

  cellnumbers.default(x, query)  ## needed for points to pass through
  #   stop(sprintf("%x not supported", class(x)))
}



#' @importFrom spatstat.geom owin as.owin pixellate
as.owin.BasicRaster <- function(W, ...) {
  msk <- matrix(TRUE, nrow(W), ncol(W))
  spatstat.geom::owin(c(raster::xmin(W), raster::xmax(W)), c(raster::ymin(W), raster::ymax(W)), mask = msk)
}
pix <- function(psp, ras) {
  spatstat.geom::pixellate(psp, as.owin(ras), weights = 1)   
}

line_cellnumbers <- function(x, r) {
  sc <- silicate::SC0(x)
  xy <- as.matrix(silicate::sc_vertex(sc)[c("x_", "y_")])
  l <- vector("list", nrow(silicate::sc_object(sc)))
  ow <- spatstat.geom::owin(range(xy[,1]), range(xy[,2]))
  for (i in seq_along(l)) {
    segs <- as.matrix(do.call(rbind, sc$object$topology_[i])[c(".vx0", ".vx1")])
    pspii <- spatstat.geom::psp(xy[segs[,1L],1L], xy[segs[,1L],2L], 
                           xy[segs[,2L],1L], xy[segs[,2L],2L], window = ow)
    
      im <- (raster(pix(pspii, r)) > 0)
    l[[i]] <- tibble(object_ = i, cell_ =   which(raster::values(im) > 0))    
  }
  dplyr::bind_rows(l)
} 
r-gris/tabularaster documentation built on Jan. 19, 2024, 8:57 a.m.