basic-use.R
In SDMtune: Species Distribution Model Selection

## ----knitr-options, include=FALSE---------------------------------------------
knitr::opts_chunk$set(comment = "#>",
                      collapse = TRUE,
                      eval = FALSE,
                      fig.align = "center")

## ----load-pkgs----------------------------------------------------------------
# library(ggplot2)    # To plot locations
# library(maps)       # To access useful maps
# library(rasterVis)  # To plot raster objects

## ----get-predictors-----------------------------------------------------------
# files <- list.files(path = file.path(system.file(package = "dismo"), "ex"),
#                     pattern = "grd",
#                     full.names = TRUE)

## ----create-raster-object-----------------------------------------------------
# predictors <- terra::rast(files)

## ----names-predicrtors--------------------------------------------------------
# names(predictors)

## ----plot-bio1----------------------------------------------------------------
# gplot(predictors$bio1) +
#     geom_tile(aes(fill = value)) +
#     coord_equal() +
#     scale_fill_gradientn(colours = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
#                          na.value = "transparent",
#                          name = "°C x 10") +
#     labs(title = "Annual Mean Temperature",
#          x = "longitude",
#          y = "latitude") +
#     scale_x_continuous(expand = c(0, 0)) +
#     scale_y_continuous(expand = c(0, 0)) +
#     theme_minimal() +
#     theme(plot.title = element_text(hjust = 0.5),
#           axis.ticks.x = element_blank(),
#           axis.ticks.y = element_blank())

## ----load-SDMtune-------------------------------------------------------------
# library(SDMtune)

## ----help-dataset-------------------------------------------------------------
# help(virtualSp)
# p_coords <- virtualSp$presence
# bg_coords <- virtualSp$background

## ----plot-presence------------------------------------------------------------
# ggplot(map_data("world"), aes(long, lat)) +
#     geom_polygon(aes(group = group), fill = "grey95", color = "gray40", size = 0.2) +
#     geom_jitter(data = p_coords, aes(x = x, y = y), color = "red",
#                 alpha = 0.4, size = 1) +
#     labs(x = "longitude", y = "latitude") +
#     theme_minimal() +
#     theme(legend.position = "none") +
#     coord_fixed() +
#     scale_x_continuous(limits = c(-125, -32)) +
#     scale_y_continuous(limits = c(-56, 40))

## ----plot-bg_model-locations--------------------------------------------------
# ggplot(map_data("world"), aes(long, lat)) +
#     geom_polygon(aes(group = group), fill = "grey95", color = "gray40", size = 0.2) +
#     geom_jitter(data = as.data.frame(bg_coords), aes(x = x, y = y),
#                 color = "blue", alpha = 0.4, size = 0.5) +
#     labs(x = "longitude", y = "latitude") +
#     theme_minimal() +
#     theme(legend.position = "none") +
#     coord_fixed() +
#     scale_x_continuous(limits = c(-125, -32)) +
#     scale_y_continuous(limits = c(-56, 40))

## ----prepare-SWD--------------------------------------------------------------
# data <- prepareSWD(species = "Virtual species",
#                    p = p_coords,
#                    a = bg_coords,
#                    env = predictors,
#                    categorical = "biome")
# 

## ----show-SWD-object----------------------------------------------------------
# data

## ----show-data----------------------------------------------------------------
# head(data@data)

## ----show-coords-data---------------------------------------------------------
# head(data@coords)

## ----show-species-------------------------------------------------------------
# data@species

## ----save-SWD-single, eval=FALSE----------------------------------------------
# swd2csv(data,
#         file_name = "data.csv")

## ----save-SWD-double, eval=FALSE----------------------------------------------
# swd2csv(data,
#         file_name = c("presence.csv", "background.csv"))

## ----train--------------------------------------------------------------------
# model <- train(method = "Maxent",
#                data = data)

## ----print-model--------------------------------------------------------------
# model

## ----get-slots----------------------------------------------------------------
# slotNames(model)

## ----model-slots--------------------------------------------------------------
# slotNames(model@model)

## ----retrain------------------------------------------------------------------
# model <- train(method = "Maxent",
#                data = data,
#                fc = "lqph",
#                reg = 1,
#                iter = 500)

## ----model-witout-default-arguments-------------------------------------------
# model <- train(method = "Maxent",
#                data = data,
#                fc = "lh",
#                reg = 0.5,
#                iter = 700)

## ----predict-train------------------------------------------------------------
# pred <- predict(model,
#                 data = data,
#                 type = "cloglog")

## ----print-pred---------------------------------------------------------------
# head(pred)

## ----predict-presence---------------------------------------------------------
# p <- data@data[data@pa == 1, ]
# pred <- predict(model,
#                 data = p,
#                 type = "cloglog")
# tail(pred)

## ----predict-raster-----------------------------------------------------------
# map <- predict(model,
#                data = predictors,
#                type = "cloglog")

## ----print-raster-output------------------------------------------------------
# map

## ----save-map, eval=FALSE-----------------------------------------------------
# map <- predict(model,
#                data = predictors,
#                type = "cloglog",
#                file = "my_map",
#                format = "GTiff")

## ----exercise-----------------------------------------------------------------
# # First create the extent that surrounds Chile
# e = terra::ext(c(-77, -60, -56, -15))
# # Now use the extent to make the prediction
# map_e <- predict(model,
#                  data = predictors,
#                  type = "cloglog",
#                  extent = e)

## ----plot-map-default---------------------------------------------------------
# plotPred(map)

## ----plot-map-custom----------------------------------------------------------
# plotPred(map,
#          lt = "Habitat\nsuitability",
#          colorramp = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"))

## ----thresholds---------------------------------------------------------------
# ths <- thresholds(model,
#                   type = "cloglog")
# ths

## ----plot-pa------------------------------------------------------------------
# plotPA(map,
#        th = ths[3, 2])

## ----plot-and-save-pa, eval=FALSE---------------------------------------------
# plotPA(map,
#        th = ths[3, 2],
#        filename = "my_pa_map",
#        format = "GTiff")

## ----auc----------------------------------------------------------------------
# auc(model)

## ----plot-roc-----------------------------------------------------------------
# plotROC(model)

## ----tss----------------------------------------------------------------------
# tss(model)

## ----aicc---------------------------------------------------------------------
# aicc(model,
#      env = predictors)

## ----train-test---------------------------------------------------------------
# library(zeallot)  # For unpacking assignment
# c(train, test) %<-% trainValTest(data,
#                                  test = 0.2,
#                                  only_presence = TRUE,
#                                  seed = 25)

## ----trai-with-train-dataset--------------------------------------------------
# model <- train("Maxent",
#                data = train)

## ----evaluate-test------------------------------------------------------------
# auc(model)
# auc(model, test = test)

## ----plot-AUC-train-test------------------------------------------------------
# plotROC(model, test = test)

## ----experiment---------------------------------------------------------------
# # Create an empty data.frame
# output <- data.frame(matrix(NA, nrow = 10, ncol = 3))
# colnames(output) <- c("seed", "trainAUC", "testAUC")
# 
# # Create 10 different random seeds
# set.seed(25)
# seeds <- sample.int(1000, 10)
# 
# # Loop through the seeds
# for (i in 1:length(seeds)) {
#   # Make the train/test split
#   c(train, test) %<-% trainValTest(data,
#                                    test = 0.2,
#                                    seed = seeds[i],
#                                    only_presence = TRUE)
# 
#   # train the model
#   m <- train("Maxent", data = train)
# 
#   # Populate the output data.frame
#   output[i, 1] <- seeds[i]
#   output[i, 2] <- auc(m)
#   output[i, 3] <- auc(m, test = test)
# }
# 
# # Print the output
# output
# # compute the range of the testing AUC
# range(output[, 3])

## ----random-folds-------------------------------------------------------------
# folds <- randomFolds(data,
#                      k = 4,
#                      only_presence = TRUE,
#                      seed = 25)

## ----cv, eval=FALSE-----------------------------------------------------------
# cv_model <- train("Maxent",
#                   data = data,
#                   folds = folds)
# cv_model

## ----cv-auc-------------------------------------------------------------------
# auc(cv_model)
# auc(cv_model, test = TRUE)

## ----enmeval-block------------------------------------------------------------
# library(ENMeval)
# block_folds <- get.block(occ = data@coords[data@pa == 1, ],
#                          bg.coords = data@coords[data@pa == 0, ])
# 
# model <- train(method = "Maxent",
#                data = data,
#                fc = "l",
#                reg = 0.8,
#                folds = block_folds)

## ----enmeval-checherboard1----------------------------------------------------
# cb_folds <- get.checkerboard1(occ = data@coords[data@pa == 1, ],
#                               env = predictors,
#                               bg.coords = data@coords[data@pa == 0, ],
#                               aggregation.factor = 4)
# 
# model <- train(method = "Maxent",
#                data = data,
#                fc = "l",
#                reg = 0.8,
#                folds = cb_folds)

## ----blockCV------------------------------------------------------------------
# library(blockCV)
# # Create sf object
# sf_df <- sf::st_as_sf(cbind(data@coords, pa = data@pa),
#                       coords = c("X", "Y"),
#                       crs = terra::crs(predictors,
#                                        proj = TRUE))
# 
# # Spatial blocks
# spatial_folds <- cv_spatial(x = sf_df,
#                             column = "pa",
#                             rows_cols = c(8, 10),
#                             k = 5,
#                             hexagon = FALSE,
#                             selection = "systematic")
# 
# model <- train(method = "Maxent",
#                data = data,
#                fc = "l",
#                reg = 0.8,
#                folds = spatial_folds)