##### Supplementary code file 1/5  PREPARE SPATIAL DATA

##### M. Cardillo "Clarifying the relationship between body size and extinction risk in amphibians by complete mapping of model space"
##### marcel.cardillo@anu.edu.au
##### Apologies for the messiness and inelegance of the code!


library(rgdal)
library(rgeos)
library(sp)
library(maptools)
library(maps)
library(RColorBrewer)
library(raster)
library(rworldmap)


##################################################################################
######  Read in map of the world:

cea <- '+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30 +ellps=WGS84 +datum=WGS84 +units=m'

data(coastsCoarse) # from rworldmap
world.cont <- countriesCoarseLessIslands
proj4string(world.cont)<-CRS(cea)
world.cont <- spTransform(world.cont, CRS=CRS(cea))
world.cont.2 <- gUnaryUnion(world.cont,id=world.cont@data$continent)

##################################################################################
###### read in amphibians shapefile

path <- "C:/-" # directory for amphibian ranges shapefile from IUCN
amph_maps <- readOGR(dsn=path, layer='AMPHIBIANS',verbose=FALSE)
proj4string(amph_maps)<-CRS(cea)
amph_maps <- spTransform(amph_maps, CRS=CRS(cea))


binoms <- unique(amph_maps@data$binomial)
amphranges <- data.frame("binomial"=as.character(binoms)"rangesize"=NA,"latcentroid"=NA)
plot(world.cont.2,axes=TRUE)
for(i in 1:length(binoms)){
    curr.map <- amph_maps[amph_maps@data$binomial%in%binoms[i],]
    amphranges$latcentroid[i] <- gCentroid(curr.map,byid=TRUE)@coords[2]
    }

amphranges$binomial <- gsub(" ","_",amphranges$binomial)

##################################################################################
###### overlay amphibian ranges with continents and extract species list for each continent

amphranges$continent <-NA
plot(world.cont.2,axes=TRUE)
for(i in 1:length(amphranges$binomial)){
    tryCatch({
    print(i)
    curr.map <- amph_maps[amph_maps@data$binomial%in%binoms[i],]
    plot(curr.map,add=TRUE,col="orange")
    test<- intersect(curr.map,world.cont.2)
    if(!is.null(test))amphranges$continent[i] <- as.character(test@data$continent)
    }, error=function(e){})
    }
    

##################################################################################
###### Define a global raster grid with a specified resolution:

extent <- c(bbox(world.cont.2)[1,],bbox(world.cont.2)[2,])
extent.lat <- extent[2]-extent[1]
extent.long <- extent[4]-extent[3]

xy <- matrix(rnorm(10000),extent.long/14617.2,extent.lat/34735.06)  
image(xy)
rast <- raster(xy)
extent(rast) <- extent
projection(rast) <- CRS(cea)
amph_rast <- rasterize(amph_maps,rast)

##################################################################################
###### read in climate layers from worldclim:

w <- getData('worldclim', var='bio', res=10)  # the 19 bioclim layers   
mpdq <- w@layers[[17]]   # mean precipitation driest quarter(bio17)
mpdq <- projectRaster(from=mpdq, to=rast)

##################################################################################
###### human footprint & aridity data:

path <- "C:/-" # directory for HF dataset   
hf <- raster(path)
hf <- projectRaster(from=hf, to=rast)

path <- "C:/-" # directory for the aridity index dataset    
aridity <- raster(path)
aridity <- projectRaster(from=aridity, to=rast)


##################################################################################
###### extract climate layer values within species polygons:

amphranges$mpdq<-NA
amphranges$hf<-NA
amphranges$aridity<-NA

tryCatch({
for(i in 1:length(binoms)){
    print(i)
    curr.map <- amph_maps[amph_maps@data$binomial%in%binoms[i],]
    curr.map <- spTransform(curr.map, CRS=CRS(cea))
    curr.map <- gUnaryUnion(curr.map)
    plot(world.cont.2,axes=TRUE)
    plot(curr.map,add=TRUE,col="orange")
    e <- extract(aridity,curr.map)
    amphranges$aridity[i] <- mean(e[[1]])
    }
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})