[GIS] From a Grid map to a dataframe of presence/absence of each cell

Question

Working in R, I have a grid map of the world (as SpatialPolygonsDataFrame object). I have also many indipendent shapefiles (formed by several polygons) expressing species ranges. I grouped all these shapefiles representing distributions in a single list.

species <- list(sp1, sp2, sp3, sp4, sp5)

What I want to obtain as output is a dataframe where for each cell of my grid (put in the rows of the dataframe) I have in the columns the presence/absence (expressed by the values 1 and 0, respectively) of all species of my list. See the example I want to obtain:

Answer 1

The steps to achieve what you want could be:

• Intersect the geometries (world grid and species polygons home ranges)
• Merge the intersected data with world grid original data
• Make a table and transform to data frame

I'm using sf package to manage spatial objects, rworldmap to get example data and mapview to make interactive map plots.

Note: I'm simulating the species home range polygons as countries boundaries

Check the commented code below.

# Load libraries 
library(sf)
library(rworldmap)
library(mapview)

# Load data 
worldMap <- getMap(resolution = "less islands")
worldMap <- st_as_sf(worldMap)

# Create example grid
worldGrid <- st_make_grid(x = worldMap, what = "polygons", cellsize = 10)
worldGrid <- st_sf(idcell = 1:length(worldGrid), geom = worldGrid) %>% 
  st_cast("POLYGON")

# View map
mapview(worldMap, zcol = 'NAME') + 
  mapview(worldGrid, alpha.regions = 0)

# Create example polygons from countries
sp1 <- worldMap[which(worldMap$NAME == "Uruguay"), "NAME"]
sp1$NAME <- "sp1"

sp2 <- worldMap[which(worldMap$NAME == "Brazil"), "NAME"]
sp2$NAME <- "sp2"

sp3 <- worldMap[which(worldMap$NAME == "Niger"), "NAME"]
sp3$NAME <- "sp3"

sp4 <- worldMap[which(worldMap$NAME == "China"), "NAME"]
sp4$NAME <- "sp4"

sp5 <- worldMap[which(worldMap$NAME == "Australia"), "NAME"]
sp5$NAME <- "sp5"

# View map
mapview(worldMap, col.regions = 'gray20') + 
  mapview(worldGrid, alpha.regions = 0) + 
  mapview(sp1) +
  mapview(sp2) +
  mapview(sp3) +
  mapview(sp4) + 
  mapview(sp5)

# Bind species polygons
species <- rbind(sp1, sp2, sp3, sp4, sp5)

# Get intersection between geometries
intersection <- st_intersection(worldGrid, species)

# View map
mapview(worldMap, col.regions = 'gray20') + 
mapview(intersection, col.regions = 'green')

# Convert intersection to df
intersection_df <- as.data.frame(intersection)[,c(1,2)]

# Merge intersection df with worldgrid
merge_df <- merge.data.frame(x = data.frame(idcell = worldGrid$idcell), y = intersection_df, all.x = TRUE, by.x = TRUE, sort = TRUE)

# Get output as desired
output <- table(merge_df)

# As data frame
output_df <- as.data.frame.matrix(cbind(idcell = worldGrid$idcell, output)) 

print(head(output_df))

idcell sp1 sp2 sp3 sp4 sp5
1      1   0   0   0   0   0
2      2   0   0   0   0   0
3      3   0   0   0   0   0
4      4   0   0   0   0   0
5      5   0   0   0   0   0
6      6   0   0   0   0   0