Implement rasterization algorithms for other shape types
noamross opened this issue · 11 comments
- Lines rasterization
- Points rasterization
- Auto detection of shape type
- Method or detection for GEOMETRYCOLLECTIONS of mutiple shape types
+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.
Using the geometric Z is an Interesting option, nice example data set here: r-spatial/sf#355
Just adding a +1 to the lines feature request!
+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.
+1 lines and points
+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
hunzikp/velox#17
Up for lines and points!
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.
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
+1 for the points feature.
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))