2

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: The output I wolud like to obtain

1

You can use the lets.presab.grid function from the package letsR:

    library(letsR)

    # Grid 
    sp.r <- rasterToPolygons(raster(resol = 5))
    slot(sp.r, "data") <- cbind("ID" = 1:length(sp.r),
                                slot(sp.r, "data"))

    # Species polygons
    data(Phyllomedusa)
    projection(Phyllomedusa) <- projection(sp.r)

    # PAM
    resu <- lets.presab.grid(Phyllomedusa, sp.r, "ID")

    # Plot
    rich_plus1 <- rowSums(resu$PAM) + 1
    colfunc <- colorRampPalette(c("#fff5f0", "#fb6a4a", "#67000d"))
    colors <- c("white", colfunc(max(rich_plus1)))
    plot(resu$grid, border = "gray40",
         col = colors[rich_plus1])
    map(add = TRUE)

Here is a complete tutorial showing how you can do it: https://rmacroecology.netlify.com/2018/01/23/a-guide-to-transform-species-shapefiles-into-a-presence-absence-matrix-based-on-a-user-defined-grid-system/

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

2

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)

map1

# 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)

map2

# 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')

map3

# 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

table

  • Thank you for your very detailed reply. The shapefiles of my species (sp1, sp2, ...) are in WGS84, therefore the intersection of your code failed to work (Error message: Error: st_crs(x) == st_crs(y) is not TRUE). For several reasons, I would like to work in my native CRS (i.e. WGS84: 4326): is it possible with your suggested method? – Marco Jan 13 '18 at 22:21
  • 1
    @Marco hope the answer can help you. The crs of the sf objects used in my answer are in wgs84 or 4326 code. You can check it using st_crs(sfobject). So I think yes, it is possible. Check if your sf objects have the crs 4326 assigned. – Guzmán Jan 14 '18 at 0:15
  • You are right. I fixed the troubles concerning crs. Your code was perfectly functioning and pertinent to my aims. Thanks! – Marco Jan 14 '18 at 23:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.