Raster cell numbers

It bugs me how raster::cellFromPolygon and co give a list, it’s a hassle to use this for further extractions, and if it was a tibble of object (polygon, line) ID and cell number it’s easy to drive extractions and summarize with the tidyverse.

Enter cellnumbers()

#' Extract cell numbers from a Raster object.
#'
#' Currently only for single layer objects.
#' @param x Raster object
#' @param ... arguments passed on
#' @param query Spatial object or matrix of coordinates
#' @return tbl_df data frame
#' @export
#' @importFrom dplyr bind_rows
#' @importFrom raster cellFromPolygon cellFromLine cellFromXY projection
#' @examples 
#' library(raster)
#' r <- raster(volcano)
#' cr <- cut(r,  pretty(values(r)))
#' p <- raster::rasterToPolygons(cr, dissolve = TRUE)
#' tt <- cellnumbers(cr, p)
#' library(dplyr)
#' tt %>% mutate(v = extract(r, cell_)) %>% 
#' group_by(object_) %>% 
#' summarize(mean(v)) %>% mutate(head(pretty(values(r)), -1))
cellnumbers <- function(x, query, ...) {
  if (inherits(query, "sf")) query <- sf::as(query, "Spatial")
  if (is.na(projection(x)) || is.na(projection(query)) || projection(x) != projection(query)) {
    warning(sprintf("projections not the same \n    x: %s\nquery: %s", projection(x), projection(query)))
  }
  if (inherits(query, "SpatialPolygons")) {
    a <- cellFromPolygon(x, query, ...)
  }
  if (inherits(query, "SpatialLines")) {
    a <- cellFromLine(x, query, ...)
  }
  if (is.matrix(query) | inherits(query, "SpatialPoints")) {
    a <- cellFromXY(x, query)
  }
  d <- dplyr::bind_rows(lapply(a, mat2d_f), .id = "object_")
  if (ncol(d) == 2L) names(d) <- c("object_", "cell_")
  if (ncol(d) == 3L) names(d) <- c("object_", "cell_", "weight_")
  d
  
}

#' @importFrom tibble as_tibble
mat2d_f <- function(x) {
  
  if (is.null(x)) {
    return(NULL)
  }
  tibble::as_tibble(x)
}


library(raster)

## Loading required package: sp

r <- raster(volcano)
cr <- cut(r,  pretty(values(r)))
p <- raster::rasterToPolygons(cr, dissolve = TRUE)

## Loading required namespace: rgeos

## get the cell numbers in tidy form, rather than a list of stuff
(tt <- cellnumbers(cr, p))

## Warning in cellnumbers(cr, p): projections not the same 
##     x: NA
## query: NA

## # A tibble: 5,307 × 2
##    object_ cell_
##      <chr> <dbl>
##  1       1     1
##  2       1     2
##  3       1     8
##  4       1     9
##  5       1    10
##  6       1  3843
##  7       1  3903
##  8       1  3904
##  9       1  3963
## 10       1  3964
## # ... with 5,297 more rows

## now apply the cell index to an extraction task
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

tt %>% mutate(v = extract(r, cell_)) %>%
group_by(object_) %>%
summarize(v = mean(v)) 

## # A tibble: 6 × 2
##   object_         v
##     <chr>     <dbl>
## 1       1  97.38693
## 2       2 110.20931
## 3       3 130.40073
## 4       4 149.72827
## 5       5 170.44589
## 6       6 186.32022

## "object_" is the character version of the row number of 'p' above (needs work to be a unique ID)

Vector object identity is a bit loose, but if you’re careful this is a really neat way to work with heavy extraction tasks.