# fitENM: Fit climatic niche models to focal palm lineages # Reference: Lim, J.Y.*, Huang, H,*, Fansworth, A., Lunt, D.J., Baker, W.J., Morley, R.J., Kissling, W.D. & Hoorn, C. (in review) Cenozoic tropical biome decline and the global diversification and biogeographic histories of palms # Author: Jun Ying Lim # Contact: junyinglim@gmail.com ## PACKAGES AND DIRECTORIES ==================== main.dir <- "~/Dropbox/Projects/2019/palms/projects/fossilPalms/" data.dir <- file.path(main.dir, "data") res.dir <- file.path(main.dir, "results") src.dir <- file.path(main.dir, "src") library(rgdal) library(spThin) library(sp) library(raster) library(maptools) library(ENMeval) source(file.path(src.dir, "utils.R")) ## IMPORT DATA =============== # Import palm occurrence data --------- calaminae_occ <- readRDS(file.path(data.dir, "calaminae_occ.rds")) eugeissona_occ <- readRDS(file.path(data.dir, "eugeissona_occ.rds")) mauritiinae_occ <- readRDS(file.path(data.dir, "mauritiinae_occ.rds")) nypa_occ <- readRDS(file.path(data.dir, "nypa_occ.rds")) # Get current climate data --------- paleoclim <- readRDS(file.path(data.dir, "paleoclim_stacks_masked.rds")) worldclimDat <- getData("worldclim", var="bio", res = 5) # 5 minute resolution worldclimPred <- stack(worldclimDat[[c(1, 4, 6, 12, 14, 15)]]) # 5, 6, 13, 14 saveRDS(worldclimPred, file.path(data.dir, "worldclimPred.rds")) # bioclim 1 = mean annual temperature # bioclim 4 = temperature seasonality # bioclim 5 = max temp of warmest month # bioclim 6 = min temp of coldest month # bioclim 12 = total annual temperature # bioclim 13 = ppt of wettest month # bioclim 14 = ppt of driest month # bioclim 15 = precipitation seasonality ## SPATIALLY THIN OCCURRENCE POINTS --------- calaminae_occ_thin <- spatialThinningByRaster(calaminae_occ, raster = worldclimPred, lat.col = "LAT", long.col = "LON") dim(calaminae_occ) dim(calaminae_occ_thin) saveRDS(calaminae_occ_thin, file.path(res.dir, "calaminae_occ_thin.rds")) nypa_occ_thin <- spatialThinningByRaster(nypa_occ, raster = worldclimPred, lat.col = "LAT", long.col = "LON") dim(nypa_occ) dim(nypa_occ_thin) saveRDS(nypa_occ_thin, file.path(res.dir, "nypa_occ_thin.rds")) eugeissona_occ_thin <- spatialThinningByRaster(eugeissona_occ, raster = worldclimPred, lat.col = "LAT", long.col = "LON") dim(eugeissona_occ) dim(eugeissona_occ_thin) saveRDS(eugeissona_occ_thin, file.path(res.dir, "eugeissona_occ_thin.rds")) mauritiinae_occ_thin <- spatialThinningByRaster(mauritiinae_occ, raster = worldclimPred, lat.col = "LAT", long.col = "LON") dim(mauritiinae_occ) dim(mauritiinae_occ_thin) saveRDS(mauritiinae_occ_thin, file.path(res.dir, "mauritiinae_occ_thin.rds")) ## CREATE BACKGROUND POINTS =============== # Sample from defined background (can also used spsample) # World clim is 0.08333 x 0.08333 degrees, ncells = 7,776,000, extent (-180, 180, -60, 90) # Paleoclim simulations are 3.75 x 2.75 degrees, ncells = 7,008, extent (-180, 180, -90, 90) # Resolution is 1350:1, 1350 = (3.75*12*2.5*12) set.seed(12345) bg_sample <- sampleRandom(worldclimPred, size = 100000, xy = TRUE, na.rm = T, cells = TRUE) saveRDS(bg_sample, file.path(res.dir, "bg_sample.rds")) ## FIT MAXENT MODELS =============== # To run with pre-generated files (caution: make sure there has not been any changes in code since their generation) worldclimPred <- readRDS(file.path(data.dir, "worldclimPred.rds")) bg_sample <- readRDS(file.path(res.dir, "bg_sample.rds")) # Notes: # Parallel functionality not used as it is buggy # By default, ENMevaluate tests RM between 0.5 to 4 # worldclim resolution is about 0.083333 degree which is approximately 10 km. # With an aggregation factor of c(2,2), adjacent points will be in different bins but points that are > 20 km apart will also be in different bins fc = c("L", "LQ", "LQH", "LQHP", "LQHPT") #L=linear, Q=quadratic, P=product, T=threshold, and H=hinge alg = "maxnet" parallel = FALSE # parallel not used as it is still buggy taxa <- c("nypa", "eugeissona", "calaminae", "mauritiinae") occ_thin_list <- list(nypa_occ_thin, eugeissona_occ_thin, calaminae_occ_thin, mauritiinae_occ_thin) # Run analysis with global background for(i in 1:length(taxa)){ print(paste("Fitting model to occurrence points of", taxa[i], sep = " ")) enmModel <- ENMevaluate(occ = occ_thin_list[[i]][,c("LON", "LAT")], fc = fc, env = worldclimPred, bg.coords = bg_sample[,c("x","y")], method = "randomkfold",#validation, kfolds = 10, algorithm = alg, rasterPreds = T, updateProgress = T) saveRDS(enmModel, file.path(res.dir, paste0(taxa[i], "_enm_kfolds_globalbg.rds"))) }