Maps and Matrices

Purpose: Create species richness maps of the rodents of Chile for all rodent species and for only rodent species that have ectoparasite associations.

Load packages

library(raster)
library(RColorBrewer)
library(sf)
library(phyloraster)
library(terra)
library(dplyr)

Read in shape file with the polygons of ranges for rodent species in Chile. The file being read in as rodrange is a shapefile with polygons of ranges for 65 rodent species in Chile from IUCN, and from gbif for Rattus rattus and Rattus norvegicus. From this, the rodents which do not have ectoparasites are being excluded and named “ectorodents”.And then the most central rodents are also sorted out.

rodrange <- read_sf("./mergedChile/mergedChile.shp", stringsAsFactors = F) # all rodents

ectorodents <- subset(rodrange, SCI_NAME !="Akodon iniscatus" & SCI_NAME != "Abrothrix jelskii" 
                  & SCI_NAME !="Auliscomys pictus" & SCI_NAME != "Calomys lepidus" 
                  & SCI_NAME !="Cavia tschudii" & SCI_NAME != "Chelemys megalonyx" 
                  & SCI_NAME !="Chinchilla chinchilla" & SCI_NAME != "Chinchilla lanigera"
                  & SCI_NAME !="Ctenomys coyhaiquensis" & SCI_NAME != "Eligmodontia morgani" 
                  & SCI_NAME !="Euneomys petersoni" & SCI_NAME != "Galea musteloides" 
                  & SCI_NAME !="Galenomys garleppi" & SCI_NAME != "Geoxus annectens"
                  & SCI_NAME !="Lagidium wolffsohni" & SCI_NAME != "Microcavia australis" 
                  & SCI_NAME !="Microcavia niata" & SCI_NAME != "Microtus juldaschi"
                  & SCI_NAME != "Neotomys ebriosus"
                  & SCI_NAME !="Octodon pacificus" & SCI_NAME != "Oligoryzomys magellanicus" 
                  & SCI_NAME !="Ondatra zibethicus" & SCI_NAME != "Phyllotis limatus") #exclude rodents without ectoparasites


badrodents <- subset(rodrange, SCI_NAME == "Rattus norvegicus" |
                       SCI_NAME == "Chelemys macronyx" |
                       SCI_NAME == "Octodon degus" |
                       SCI_NAME == "Phyllotis xanthopygus" |
                       SCI_NAME == "Aconaemys porteri"| 
                        SCI_NAME == "Loxodontomys micropus" |
                        SCI_NAME == "Reithrodon auritus" |
                        SCI_NAME == "Oligoryzomys longicaudatus" |
                        SCI_NAME == "Rattus rattus" |
                        SCI_NAME == "Abrothrix longipilis" |
                       SCI_NAME == "Phyllotis darwini" |
                      SCI_NAME == "Abrothrix olivaceus" ) # central rodents

Here, we use the phyloraster package to map species richness maps for the all the rodents, rodents with ectoparasites, and the most central rodents.

v_all <-vect(rodrange) #make spatvectors 
v_ecto <- vect(ectorodents)
v_bad <- vect(badrodents)

r_all <- shp2rast(v_all, sps.col = "SCI_NAME", ymask = FALSE, background = 0, 
             resolution = 0.1) #make rasters
r_ecto <- shp2rast(v_ecto, sps.col = "SCI_NAME", ymask = FALSE, background = 0, 
             resolution = 0.1)
r_bad <- shp2rast(v_bad, sps.col = "SCI_NAME", ymask = FALSE, background = 0, 
             resolution = 0.1)

sr_all<- rast.sr(x=r_all) #use phyloraster function to get species richness
sr_ecto <- rast.sr(x=r_ecto)
sr_bad <- rast.sr(x=r_bad)

chile <- read_sf("./CHL_adm/CHL_adm0.shp") #CHile from DIVA gis
sr_allm<-mask(sr_all, chile) #limit to geographical extent of Chile, this mostly just helps with getting the background, when plotted, to not have a raster value
sr_ectom<- mask(sr_ecto, chile)
sr_badm<- mask(sr_bad, chile)

Plot the maps

color<-colorRampPalette(c("lightgray", "#090387ff","lightblue", "#fada5fff","orange", "#aa0000")) #make ramp palette

par(mfrow=c(1,1))
plot(sr_allm, col=color(22),
     ext=c(-81, -59, -59, -15), main="Total Rodent Species Richness",
     plg=list(cex = 1.3 ))

plot(sr_ectom,col=color(18), 
     ext=c(-81, -59, -59, -15), main="Rodent Species Richness",
    plg=list(cex = 1.3 ))

plot(sr_badm,col=color(10), 
     ext=c(-81, -59, -59, -15), main="Central Rodent Species Richness",
       plg=list(cex = 1.3 ))

The next step is to calculate pairwise range overlap between rodent species. We will do this only for the rodents with ectoparasites

sf::sf_use_s2(FALSE) #turn off spherical geometry

## Spherical geometry (s2) switched off

ectogeom<- st_sfc(ectorodents$geometry)

overlap<- lapply(ectogeom, function(x){
  lapply(ectogeom, function(y) st_intersection( x, y ) %>% st_area(x,y)/(st_area(x)+st_area(y)-st_area(x, y)) )
  })
geolap<- matrix(unlist(overlap), ncol = length(ectorodents$geometry), byrow = TRUE)
colnames(geolap)<- ectorodents$SCI_NAME
rownames(geolap)<- ectorodents$SCI_NAME

write.csv(geolap, "OverlapMatrix.csv")

Now you have a matrix with pairwise range overlap between rodents in a csv file, “OverlapMatrix.csv”.

Calculate the distance between range centroids for pairs of rodents and write a csv file with the distance matrix.

ectorodents <- ectorodents %>% mutate(centroids = st_centroid(st_geometry(.)))
geodist<- ectorodents %>% st_set_geometry('centroids') %>% st_distance()
colnames(geodist) <- ectorodents$SCI_NAME
row.names(geodist)<- ectorodents$SCI_NAME
df<- as.data.frame(geodist)

write.csv(df, "CentroidDist.csv")