Rasterize SpatialPolygons with holes in R

2
votes

I'm trying to rasterize a SpatialPolygons object in a way that I get the number of polygons that overlap with each cell of the raster. My polygons have holes. The issue I have is that if two holes overlap it includes this overlap in the count.

I've created a toy example:

library(sp)
library(raster)

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(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-180,-20), c(-130,55), c(0,60), c(40,5), c(15,-45), c(-180,-20))
p3 <- list(p3, hole)

pols <- spPolygons(p1, p2, p3)

pols@polygons[[1]]@Polygons[[2]]@hole
pols@polygons[[3]]@Polygons[[2]]@hole


r1 <- raster(ncol=180, nrow=180)
r <- rasterize(pols, r1, fun="count")
plot(pols, col=rgb(1,0,0,0.5))
plot(r)

I think this might be what was causing the problem associated to question use R to rasterize shapefile with holes? , but there is no final solutions. In question rasterize ESRI shapefile with holes but FALSE hole slots the problem is that the holes are not properly formated, but that's not my case. If you look at the holes, they are both TRUE, see:

pols@polygons[[1]]@Polygons[[2]]@hole
pols@polygons[[3]]@Polygons[[2]]@hole

Any help would be appreciated!

2
rr-rastersp

2 Answers

2
votes

Using package fasterize to do the rasterization (besides being faster) seems to avoid the problem: 

library(fasterize)
library(sf)

r1 <- raster(ncol=180, nrow=180)
r <- fasterize::fasterize(sf::st_as_sf(pols),
                          r1, fun = "count")
plot(r, col = rainbow(3))

HTH!

1
votes

Here is a workaround. Rasterize each polygon first and then sum them up. By doing this the hole would be zero.

library(sp)
library(raster)

# Create a list of Spatial Polygons
pols_list <- lapply(list(p1, p2, p3), spPolygons)

# Create a blank raster
r1 <- raster(ncol = 180, nrow = 180)

# Rasterize each polygon
r_list <- lapply(pols_list, rasterize, y = r1, fun = "count")

# Create a raster stack
s <- stack(r_list)

# Calculate the sum
r2 <- sum(s, na.rm = TRUE)

# Plot r2
plot(r2)

enter image description here