Push-pull_tech_codes/3. maxent_model_using_sdm_package.R at main · icipe-official/Push-pull_tech_codes · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
# =========================================
# Generic Species Distribution Modeling Pipeline
# Example species: Napier grass
# Author: Quinto Juma Meltus
# Purpose: Process species occurrence data and environmental predictors
#          and fit MaxEnt models with configurable parameters
# =========================================

# ------------------------------
# 1. Install required packages (only if not installed)
# ------------------------------
packages <- c('terra', 'sp', 'raster', 'rgdal', 'dismo', 'rJava', 'sdm', 'gbm', 'vip', 'pdp', 'boot', 'dplyr', 'PresenceAbsence', 'ggpubr', 'RColorBrewer', 'quickPlot', 'foreign')
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if(length(new_packages)) install.packages(new_packages, dependencies = TRUE)

# ------------------------------
# 2. Load libraries
# ------------------------------
library(terra)
library(sp)
library(raster)
library(rgdal)
library(dismo)
library(rJava)
library(sdm)
library(dplyr)

# ------------------------------
# 3. Load species occurrence data
# ------------------------------
species_file <- "path_to/Napier_grass_occurrences.shp"  # Replace with your shapefile
species <- shapefile(species_file)

# Define CRS (WGS84)
sr <- CRS("+proj=longlat +datum=WGS84 +no_defs")
species <- spTransform(species, sr)

# Plot species occurrences
plot(species, main = "Napier Grass Occurrences")

# ------------------------------
# 4. Load environmental predictors
# ------------------------------
raster_file <- "path_to/environmental_predictors.tif"  # Multi-band raster
preds <- brick(raster_file)  # If multiple bands stacked
plot(preds[[1]], main = "First Predictor Layer")
plot(species, add = TRUE)

# ------------------------------
# 5. Prepare SDM data
# ------------------------------
# Background points: 1500 random points; remove duplicates
sdm_data <- sdmData(
  formula = Occ ~ .,
  train = species,
  predictors = preds,
  bg = list(n = 1500, method = 'gRandom', remove = TRUE)
)

# ------------------------------
# 6. Check available methods in sdm
# ------------------------------
getmethodNames()

# ------------------------------
# 7. Fit MaxEnt model with options
# ------------------------------
# Configure regularization and feature types via maxent.args
maxent_args <- c(
  "betamultiplier=1.5",     # Regularization multiplier (change as needed)
  "linear=true",            # Use linear features
  "quadratic=true",         # Use quadratic features
  "product=false",
  "threshold=false",
  "hinge=true"
)

start_time <- Sys.time()
maxent_model <- sdm(
  Occ ~ .,
  data = sdm_data,
  methods = c('maxent'),
  replication = "cv",
  test.percent = 30,
  cv.folds = 5,
  maxent.args = maxent_args
)
end_time <- Sys.time()
cat("Modeling time:", end_time - start_time, "\n")

# ------------------------------
# 8. Predict across environmental layers
# ------------------------------
predicted_raster <- predict(
  maxent_model,
  newdata = preds,
  filename = "path_to/output/NapierGrass_prediction.tif",
  overwrite = TRUE
)
plot(predicted_raster, main = "Predicted Suitability for Napier Grass")

# ------------------------------
# 9. Parallel prediction (optional)
# ------------------------------
pred_parallel <- predict(
  maxent_model,
  newdata = preds,
  filename = "path_to/output/NapierGrass_prediction_parallel.tif",
  parallelSetting = list(ncore = 2),
  overwrite = TRUE
)
plot(pred_parallel)

# ------------------------------
# 10. Ensemble modeling (optional)
# ------------------------------
ensemble_pred <- ensemble(
  maxent_model,
  newdata = preds,
  filename = "path_to/output/NapierGrass_ensemble.tif",
  setting = list(method = 'weighted', stat = 'AUC'),
  overwrite = TRUE
)
plot(ensemble_pred, main = "Ensemble Prediction")

# ------------------------------
# 11. Variable importance
# ------------------------------
varImp <- getVarImp(maxent_model)
print(varImp)
plot(varImp)

# ------------------------------
# 12. Model evaluation
# ------------------------------
rmse_vals <- getEvaluation(maxent_model, stat = 'RMSE')
print(rmse_vals)

# ------------------------------
# 13. Save evaluation metrics
# ------------------------------
write.csv(rmse_vals, "path_to/output/NapierGrass_model_evaluation.csv", row.names = FALSE)
cat("SDM workflow complete. Predictions and evaluations saved.\n")