############################ # a) ecospat.occupied.patch ############################ ### calculates the occupied patches of a binary SDM suitability map ## FUNCTION'S ARGUMENTS ## bin.map: Binary map (single layer or raster stack) from a Species Distribution Model. ## Sp.occ.xy: xy-coordinates of the species presence ## buffer: numeric. Calculate occupied patch models from the binary map using a buffer (predicted area occupied by the species or within a buffer around the species). ## Details: ## Value # a RasterLayer with value 1 () ## Author(s) # Frank Breiner ecospat.occupied.patch <- function(bin.map, Sp.occ.xy, buffer = 0){ cl<-clump(bin.map) d<-extract(cl,Sp.occ.xy,buffer=buffer) d <- unique(na.omit(unlist(d))) b.map<-(cl %in% d) + bin.map names(b.map) <- names(bin.map) return(b.map)} ## EXAMPLE library(dismo) # only run if the maxent.jar file is available, in the right folder jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep='') # checking if maxent can be run (normally not part of your script) file.exists(jar) require(rJava) # get predictor variables fnames <- list.files(path=paste(system.file(package="dismo"), '/ex', sep=''), pattern='grd', full.names=TRUE ) predictors <- stack(fnames) #plot(predictors) # file with presence points occurence <- paste(system.file(package="dismo"), '/ex/bradypus.csv', sep='') occ <- read.table(occurence, header=TRUE, sep=',')[,-1] colnames(occ) <- c("x","y") occ <- ecospat.mindist(occ,mindist=1) # fit model, biome is a categorical variable me <- maxent(predictors, occ, factors='biome') # predict to entire dataset pred <- predict(me, predictors) plot(pred) points(occ) # use MPA to convert suitability to binary map mpa.cutoff <- ecospat.mpa(pred,occ) # use Boyce index to convert suitability to binary map boyce <- ecospat.boyce(pred, occ) ### use the boyce index to find a threshold boyce.cutoff <- ecospat.boyce.cutoff(boyce) pred.bin.mpa <- ecospat.binary.model.mod(pred,mpa.cutoff) names(pred.bin.mpa) <- "me.mpa" pred.bin.boyce <- ecospat.binary.model.mod(pred,max(boyce.cutoff)) names(pred.bin.boyce) <- "me.boyce" mpa.ocp <- ecospat.occupied.patch(pred.bin.mpa,occ) boyce.ocp <- ecospat.occupied.patch(pred.bin.boyce,occ) par(mfrow=c(1,2)) plot(mpa.ocp) ## occupied patches: green area points(occ,col="red",cex=0.5,pch=19) plot(boyce.ocp) points(occ,col="red",cex=0.5,pch=19) ## with buffer: mpa.ocp <- ecospat.occupied.patch(pred.bin.mpa,occ, buffer=500000) boyce.ocp <- ecospat.occupied.patch(pred.bin.boyce,occ, buffer=500000) plot(mpa.ocp) ## occupied patches: green area points(occ,col="red",cex=0.5,pch=19) plot(boyce.ocp) points(occ,col="red",cex=0.5,pch=19) ############################ # b) ecospat.rangesize ############################ ## This function estimates the range size of a species using standard IUCN criteria (AOO, EOO) ## or binary predictions of Species Distribution Models ## FUNCTION'S ARGUMENTS ## bin.map: Binary map (single layer or raster stack) from a Species Distribution Model. ## ocp: logical. Calculate occupied patch models from the binary map (predicted area occupied by the species). ## buffer: numeric. Calculate occupied patch models from the binary map using a buffer (predicted area occupied by the species or within a buffer around the species). ## eoo.around.model: logical. The EOO around all positive predicted values from the binary map ## eoo.around.modelocp: logical. EOO around all positive predicted values of occupied patches. ## xy: xy-coordinates of the species presence ## EOO: logical. Extent of Occurrence. Convex Polygon around occurrences. ## Model.within.eoo: logical. Postitive predicted Area within EOO. ## AOO: logical. Area of Occurrences. ## resol: Resolution of the grid frame at which AOO should be calculated. ## AOO.circles: logical. AOO calculated by circles around the occurrences instead of using a grid frame. ## d: Radius of the AOO.circles around the occurrences. ## lonlat: Are these longitude/latidue coordinates? (Default FALSE) ## alphahull: logical. Calculates the alpha-hull for the occurrences ## alpha: numeric. ## return.obj: logical. should the objects created to estimate range size be returned (rasterfiles and spatial polygons) ## save.obj: logical. should opjects be saved on hard drive? ## save.rangesize: logical. should range size estimations be saved on hard drive ## directory: directory toin which objects should be saved (Default: getwd()) ## Details: ## Value # ## Author(s) # Frank Breiner ## See Also # convHull, circles, ahull, ecospat.occupied.patch ecospat.rangesize <- function(bin.map = NULL, ocp = T, buffer = 0, eoo.around.model = T, eoo.around.modelocp = F, xy = NULL, EOO = T, Model.within.eoo = T, AOO = T, resol = c(2000, 2000), AOO.circles = F, d = sqrt(2000^2/pi), lonlat = F, alphahull= F, alpha = 2, return.obj = T, save.obj = F, save.rangesize = F, directory = getwd()){ #if(ocp | EOO | alphahull | AOO & is.null(xy)){stop("need xy-coordinates of species occurrence")} if(!require(dismo)){stop("dismo package required!")} ### IUCN methods: if(Model.within.eoo | EOO){ xy.eoo<-convHull(xy) eoo <- round(xy.eoo@polygons@polygons[[1]]@area) #quadrat km }else{xy.eoo <- eoo <- NULL} if(AOO){ aoo.obj <- bin.map[[1]] res(aoo.obj)<- resol aoo.obj[] <- 0 aoo.obj <- rasterize(xy, aoo.obj, field=1) aoo.rs <- as.numeric(freq( aoo.obj==1)[1,2]*prod(resol)) ### km2 AOO2 }else{aoo.obj <- aoo.rs <- NULL} if(AOO.circles){ circ <- dismo::circles(xy,d=d,lonlat=lonlat) circ.rs <- round(raster::area(circ@polygons)) }else{circ.rs <- circ <- NULL} if(alphahull){ if(!require(alphahull)){stop("alphahull package required!")} del<-delvor(xy) dv<-del$mesh mn <- mean(sqrt(abs(del$mesh[,3]-del$mesh[,5])^2+abs(del$mesh[,4]-del$mesh[,6])^2))*alpha h<-ahull(del,alpha=mn) alpha.hull <- round(areaahull(h)) }else{h <- alpha.hull <- NULL} ## Model range size estimation bin.map.rs <-NULL bin.map.ocp.rs <- bin.map.ocp <- NULL eoo.around.mo <- list();eoo.around.mo.rs <- NULL eoo.around.ocp <- NULL eoo.around.mo.ocp <- list();eoo.around.mo.ocp.rs <- NULL map.mo.within.eoo <- mo.within.eoo <- NULL if(!is.null(bin.map)){ for(i in 1:nlayers(bin.map)){ bin.m <- bin.map[[i]] bin.map.rs <- c(bin.map.rs,sum(bin.m[bin.m==1])*prod(res(bin.m))) names(bin.map.rs) <- names(bin.map[[1:i]]) if(ocp | eoo.around.modelocp){ if(i==1){ bin.map.ocp <- ecospat.occupied.patch(bin.map[[i]],xy, buffer=buffer) }else{ bin.map.ocp <- addLayer(bin.map.ocp, ecospat.occupied.patch(bin.map[[i]],xy, buffer=buffer)) } bin.map.ocp.rs <- c(bin.map.ocp.rs,(sum(bin.map.ocp[[i]][bin.map.ocp[[i]]==2])/2)*prod(res(bin.map.ocp))) names(bin.map.ocp.rs) <- paste(names(bin.map[[1:i]]),"ocp",sep="_") } if(eoo.around.model | eoo.around.modelocp){ ids <- init(bin.map[[i]], v='cell') # extract these cells <- extract(ids, ids[bin.map[[i]]==1]) if(eoo.around.model){ xy.model.eoo <- xyFromCell(ids,ids[bin.map[[i]]==1]) eoo.around.mo[[i]] <- convHull(xy.model.eoo) eoo.around.mo.rs <- c(eoo.around.mo.rs,round(eoo.around.mo[[i]]@polygons@polygons[[1]]@area)) names(eoo.around.mo.rs) <- paste(names(bin.map[[1:i]]),"eoo.around.model",sep="_") names(eoo.around.mo) <- paste(names(bin.map[[1:i]]),"eoo.around.model",sep="_") } } if(eoo.around.modelocp){ xy.model.eoo.ocp <- xyFromCell(ids,ids[bin.map.ocp[[i]]==2]) eoo.around.mo.ocp[[i]] <- convHull(xy.model.eoo.ocp) eoo.around.mo.ocp.rs <- c(eoo.around.mo.ocp.rs,round(eoo.around.mo.ocp[[i]]@polygons@polygons[[1]]@area)) names(eoo.around.mo.ocp.rs) <- paste(names(bin.map[[1:i]]),"eoo.around.model.ocp",sep="_") names(eoo.around.mo.ocp) <- paste(names(bin.map[[1:i]]),"eoo.around.model.ocp",sep="_") } if(Model.within.eoo){ d<-extract(bin.map[[i]], xy.eoo@polygons,cellnumbers=T) mo.within.eoo <- c(mo.within.eoo,round(sum(d[[1]][,2],na.rm=T)*prod(res(bin.map[[i]])))) cells <- d[[1]][,1][d[[1]][,2]==1 & !is.na(d[[1]][,2])] if(return.obj){ if(i==1){ map.mo.within.eoo <- bin.map[[i]] map.mo.within.eoo[cells] <-2 }else{ map.mo.within.eoo1 <- bin.map[[i]] map.mo.within.eoo1[cells] <-2 map.mo.within.eoo <- addLayer(map.mo.within.eoo, map.mo.within.eoo1) } } names(mo.within.eoo) <- paste(names(bin.map[[1:i]]),"mo.within.eoo",sep="_") } } } range.sizes <-list(AOO = aoo.rs, AOO.circle = circ.rs, alpha.hull = alpha.hull, EOO = eoo, model=bin.map.rs, models.ocp= bin.map.ocp.rs, eoo.around.model = eoo.around.mo.rs, eoo.around.mo.ocp =eoo.around.mo.ocp.rs, model.within.eoo = mo.within.eoo) objects <- list(AOO = aoo.obj, AOO.circle = circ, alpha.hull = h, EOO = xy.eoo, models.ocp = bin.map.ocp, eoo.around.model = eoo.around.mo, eoo.around.mo.ocp =eoo.around.mo.ocp, model.within.eoo = map.mo.within.eoo) if(return.obj){ return(list(RangeSize= range.sizes, RangeObjects=objects)) }else{ return(list(RangeSize= range.sizes)) } if(save.obj){ save(objects, file=paste(directory,"rangesize_objects",sep="/")) } if(save.obj){ save(save.rangesize, file=paste(directory,"rangesize",sep="/")) } } ## EXAMPLE library(dismo) # only run if the maxent.jar file is available, in the right folder jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep='') # checking if maxent can be run (normally not part of your script) file.exists(jar) require(rJava) # get predictor variables fnames <- list.files(path=paste(system.file(package="dismo"), '/ex', sep=''), pattern='grd', full.names=TRUE ) predictors <- stack(fnames) #plot(predictors) # file with presence points occurence <- paste(system.file(package="dismo"), '/ex/bradypus.csv', sep='') occ <- read.table(occurence, header=TRUE, sep=',')[,-1] colnames(occ) <- c("x","y") occ <- ecospat.mindist(occ,mindist=1) # fit model, biome is a categorical variable me <- maxent(predictors, occ, factors='biome') # predict to entire dataset pred <- predict(me, predictors) plot(pred) points(occ) # use MPA to convert suitability to binary map mpa.cutoff <- ecospat.mpa(pred,occ) # use Boyce index to convert suitability to binary map boyce <- ecospat.boyce(pred, occ) ### use the boyce index to find a threshold boyce.cutoff <- ecospat.boyce.cutoff(boyce,T) pred.bin.mpa <- ecospat.binary.model.mod(pred,mpa.cutoff) names(pred.bin.mpa) <- "me.mpa" pred.bin.boyce <- ecospat.binary.model.mod(pred,max(boyce.cutoff)) names(pred.bin.boyce) <- "me.boyce" rangesize <- ecospat.rangesize(stack(pred.bin.mpa,pred.bin.boyce), xy=occ, resol=c(1,1), eoo.around.modelocp =T, AOO.circles = F, d=200000, lonlat =T, alphahull=T) rangesize$RangeSize names(rangesize$RangeObjects) par(mfrow=c(1,3)) plot(ecospat.binary.model.mod(pred,0),legend=F, main="IUCN criteria") ## IUCN criteria & derivates # plot AOO plot(rangesize$RangeObjects$AOO,add=T, col="red",legend=F) # plot EOO plot(rangesize$RangeObjects$EOO@polygons,add=T, border="red", lwd=2) # plot circles around occurrences plot(rangesize$RangeObjects$AOO.circle@polygons,add=T,border="blue") # plot alphahulls plot(rangesize$RangeObjects$alpha.hull,add=T,lwd=1, wpoints=F) for(i in 1:2){ ## plot the occupied patches of the model plot(rangesize$RangeObjects$models.ocp[[i]],col=c("grey","blue","darkgreen"),main=names(rangesize$RangeObjects$models.ocp[[i]]),legend=F) points(occ,col="red",cex=0.5,pch=19) ## plot EOO around model plot(rangesize$RangeObjects$eoo.around.model[[i]]@polygons,add=T,border="blue",lwd=2) ## plot EOO around occupied patches plot(rangesize$RangeObjects$eoo.around.mo.ocp[[i]]@polygons,add=T,border="darkgreen",lwd=2) ## plot the modeled area within EOO #plot(rangesize$RangeObjects$model.within.eoo[[i]],col=c("grey","blue","darkgreen"),legend=F) #points(occ,col="red",cex=0.5,pch=19) #plot(rangesize$RangeObjects$EOO@polygons,add=T, border="red", lwd=2) }