Create a gist now

Instantly share code, notes, and snippets.

@amsantac /r-arcgis-spdistmod.R Secret
Last active Nov 14, 2016

What would you like to do?
R ArcGIS species distribution modeling script
tool_exec <- function(in_params, out_params)
{
# load the required packages
cat(paste0("\n", "Loading packages...", "\n"))
if (!requireNamespace("dismo", quietly = TRUE))
install.packages("dismo")
if (!requireNamespace("raster", quietly = TRUE))
install.packages("raster")
require(dismo)
require(raster)
# read input and output parameters
occurrence_dataset = in_params[[1]]
continuous_rasters_folder = in_params[[2]]
biome_raster = in_params[[3]]
model = in_params[[4]]
out_raster = out_params[[1]]
out_table = out_params[[2]]
out_shp = out_params[[3]]
cat(paste0("\n", "Loading datasets...", "\n"))
# open the input shapefile dataset
d <- arc.open(occurrence_dataset)
occurrence <- arc.data2sp(arc.select(d))
# read the continuous raster files from a folder
rfiles1 <- list.files(path = continuous_rasters_folder, full.names = TRUE)
rasters1 <- stack(rfiles1[-grep(".aux", rfiles1)])
# read the categorical raster (e.g., biome)
raster2 <- raster(gsub("/", "\\\\", biome_raster))
# create a RasterStack object for the model predictors
predictors <- stack(rasters1, raster2)
# extract values from the RasterStack at the locations of the species occurrence
presvals <- as.data.frame(extract(predictors, occurrence))
cat(paste0("\n", "Adjusting model...", "\n"))
# adjust the chosen model: bioclim, domain, Mahalanobis or a generalized linear model
if(model == "bioclim"){
fitmodel <- bioclim(subset(presvals, select = -c(biome)))
p <- predict(predictors, fitmodel)
}
if(model == "domain"){
fitmodel <- domain(dropLayer(predictors, 'biome'), as.data.frame(occurrence)[, c('coords.x1', 'coords.x2')])
p <- predict(dropLayer(predictors, 'biome'), fitmodel)
}
if(model == "mahal"){
fitmodel <- mahal(dropLayer(predictors, 'biome'), as.data.frame(occurrence)[, c('coords.x1', 'coords.x2')])
p <- predict(dropLayer(predictors, 'biome'), fitmodel)
}
if(model == "glm"){
set.seed(0)
backgr <- randomPoints(predictors, 500)
absvals <- extract(predictors, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
sdmdata[, 'biome'] = as.factor(sdmdata[, 'biome'])
fitmodel = glm(pb ~ ., data = sdmdata)
p <- predict(predictors, fitmodel)
}
# export the predicted species distribution raster object and the output table and shapefile
cat(paste0("\n", "Writing result datasets...", "\n"))
if (!is.null(out_raster) && out_raster != "NA")
writeRaster(p, out_raster)
if (!is.null(out_table) && out_table != "NA")
arc.write(out_table, presvals)
if (!is.null(out_shp) && out_shp != "NA")
arc.write(out_shp, presvals, coords = coordinates(occurrence), shape_info = arc.shapeinfo(d))
cat(paste0("\n", "Done.", "\n"))
return(out_params)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment