Implement rasterization algorithms for other shape types

[noamross]

• [ ] Lines rasterization
• [ ] Points rasterization
• [ ] Auto detection of shape type
• [ ] Method or detection for GEOMETRYCOLLECTIONS of mutiple shape types

[mikoontz]

+1 For Lines rasterization! It would complement the polygon rasterization because we could rasterize by lines, and then by polygons if we want to choose to be conservative or anti-conservative with respect to the default behavior of rasterize/fasterize on polygons.

library(raster)
library(dplyr)
library(viridis)

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
p2[, 1] <- p2[, 1] + 200
p2 <- list(p2)
pols <- st_sf(value = 1:2,
              geometry = st_sfc(lapply(list(p1, p2), st_polygon)))
pols_sp <- as(pols, "Spatial")

r <- raster(resolution = 15.5)
extent(r) <- extent(pols_sp)

r <- 
  rasterize(as(pols_sp, "SpatialLines"), r, field = rep(1, times = length(pols_sp))) %>%
  rasterize(pols_sp, ., update = TRUE, field = rep(1, times = length(pols_sp)))

plot(r, col = viridis(2))
plot(pols_sp, add = TRUE)
plot(as(r, "SpatialGrid"), add = TRUE)

Fig. 1. Anti-conservative implementation of rasterize, where a polygon value is transferred to pixels with any overlap to that polygon.

tmp <- r
tmp[] <- 0

r2 <- 
  rasterize(pols_sp, tmp, update = TRUE, field = rep(1, times = length(pols_sp))) %>%
  rasterize(as(pols_sp, "SpatialLines"), ., update = TRUE, field = rep(0, times = length(pols_sp)))

r2[r2[] == 0] <- NA
plot(r2, col = viridis(2))
plot(pols_sp, add = TRUE)
plot(as(r, "SpatialGrid"), add = TRUE)

Figure 2. Conservative implementation of rasterize, where a polygon value is only transferred to pixels that are entirely overlapped by the polygon.

[mdsumner]

Using the geometric Z is an Interesting option, nice example data set here: https://github.com/edzer/sfr/issues/355

[jhollist]

Just adding a +1 to the lines feature request!

[mattbk]

+1 for points from me. Would like to be able to go from sf points (representing several layers) to a stack of rasters more quickly.

[ctlamb]

+1 lines and points

[Martin-Jung]

+1 also for small polygons! Currently fasterize is not working for small polygons that are completely contained within a grid cell. This has been an issue for the velox package as well which has now added a "small = TRUE" parameter to the rasterization. See here. Also for testing data https://github.com/hunzikp/velox/issues/17

[LLeiSong]

Up for lines and points!

[mhweber]

I would also love to see implementation for small polygons as @Martin-Jung mentions in similar way to 'small=TRUE' parameter in velox. I'd like to use this functionality in another R package, but velox was dropped from CRAN in April 2020 and only available via GitHub install, which I'd like to avoid as a package dependency.

[mdsumner]

just a note, folks discussing features here should try stars and terra as two distinct approaches, the landscape has changed quite a bit (from one where sp and raster just weren't fast enough) and maybe those can cover some use cases now

[felippehsk]

+1 for the points feature.

[mdsumner]

IMO points is not needed, can it really be faster than point lookup?

library(terra)
r <- terra::rast(nrows = 50000, ncols = 80000)
pts <- geosphere::randomCoordinates(1e6)
system.time(cell <- cellFromXY(r, pts))
   user  system elapsed 
  0.064   0.020   0.084 

You don't even need spatial libs for this, happy to unpack how it works but terra is v fast. There'll be time taken to materialize the raster data and populate those cells, but I think it's not worth complicating the landscape more than it is. (Happy to discuss though!).

library(vaster) ## remotes::install_github("hypertidy/vaster")
ex <- c(-180, 180, -90, 90)
dm <- c(80000, 50000)
system.time(cell <- cell_from_xy(dm, extent = ex, pts))
 user  system elapsed 
  0.382   0.085   0.466 

You have to tabulate those cells to a count (say) but this is all it takes. A grouping in a dplyr call can aggreate other values onto those cell instances.

This will only work for smaller raster objects, but the same applies to fasterize with its materialized matrix, on computers I use fasterize won't work at 80000x50000 values. I think at 64Gb of RAM the largest you could do is about 10Kx10K.

##values(r) <- tabulate(cells, ncell(r))