## import packages
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(terra)
## terra 1.8.29
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(mapview)
library(rJava)
library(corrplot)
## corrplot 0.95 loaded
library(predicts)
## load data
setwd('/Users/brunetca/Documents/MEGA/phython exercice/data_python/')
python_presence <- read.table("python_presence.txt", header = TRUE)
head(python_presence)
## lon lat
## 1 90.33333 24.01667
## 2 104.04087 13.46905
## 3 118.28360 25.93977
## 4 113.08333 24.36667
## 5 112.43667 23.13833
## 6 106.59409 26.92268
str(python_presence)
## 'data.frame': 91 obs. of 2 variables:
## $ lon: num 90.3 104 118.3 113.1 112.4 ...
## $ lat: num 24 13.5 25.9 24.4 23.1 ...
ncol(python_presence)
## [1] 2
colnames(python_presence)
## [1] "lon" "lat"
python_presence_sf <-st_as_sf(python_presence, coords=c("lon","lat"),crs=4326)
env_asia_c <- raster::stack('current_asia.grd')
names(env_asia_c)
## [1] "annual_mean_temp" "mean_diurn_range" "isothermality"
## [4] "temp_seasonality" "temp_max_warm" "temp_min_cold"
## [7] "temp_annual_range" "mean_temp_wet" "mean_temp_dry"
## [10] "mean_temp_warm" "mean_temp_cold" "annual_percip"
## [13] "percip_wet" "percip_dry" "percip_seasonality"
## [16] "percip_wet_quart" "percip_dry_quart" "percip_warm_quart"
## [19] "percip_cold_quart"
plot(env_asia_c)
env_america_45 <- raster::stack('future45_america.grd')
names(env_america_45)
## [1] "annual_mean_temp" "mean_diurn_range" "isothermality"
## [4] "temp_seasonality" "temp_max_warm" "temp_min_cold"
## [7] "temp_annual_range" "mean_temp_wet" "mean_temp_dry"
## [10] "mean_temp_warm" "mean_temp_cold" "annual_percip"
## [13] "percip_wet" "percip_dry" "percip_seasonality"
## [16] "percip_wet_quart" "percip_dry_quart" "percip_warm_quart"
## [19] "percip_cold_quart"
plot(env_america_45)
env_america_85 <- raster::stack('future85_america.grd')
names(env_america_85)
## [1] "annual_mean_temp" "mean_diurn_range" "isothermality"
## [4] "temp_seasonality" "temp_max_warm" "temp_min_cold"
## [7] "temp_annual_range" "mean_temp_wet" "mean_temp_dry"
## [10] "mean_temp_warm" "mean_temp_cold" "annual_percip"
## [13] "percip_wet" "percip_dry" "percip_seasonality"
## [16] "percip_wet_quart" "percip_dry_quart" "percip_warm_quart"
## [19] "percip_cold_quart"
plot(env_america_85)
## estimate asia
env_asia_c_terra <- rast(env_asia_c) ## wandelt RasterStack in SpatRaster um
class(env_asia_c_terra)
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
##
asia_background <- nrow(python_presence) * 10
background_points_asia <- predicts::backgroundSample(
env_asia_c_terra,
n = asia_background,
p = python_presence[, c("lon", "lat")]
)
## Warning in predicts::backgroundSample(env_asia_c_terra, n = asia_background, :
## generated random points = 0.959340659340659 times requested number
background_points_asia_sf <- st_as_sf(
as.data.frame(background_points_asia),
coords = c("x", "y"),
crs = 4326
)
python_presence_sf$occ <- 1
background_points_asia_sf$occ <- 0
colnames(python_presence_sf)
## [1] "geometry" "occ"
colnames(background_points_asia_sf)
## [1] "geometry" "occ"
python_pa_sf <- rbind(python_presence_sf, background_points_asia_sf)
env_asia_extract <- terra::extract(env_asia_c, python_pa_sf, sp = TRUE)
model_df <- as.data.frame(env_asia_extract)
model_df$occ <- python_pa_sf$occ
model_df <- model_df[complete.cases(model_df), ]
colnames(model_df)
## [1] "occ" "annual_mean_temp" "mean_diurn_range"
## [4] "isothermality" "temp_seasonality" "temp_max_warm"
## [7] "temp_min_cold" "temp_annual_range" "mean_temp_wet"
## [10] "mean_temp_dry" "mean_temp_warm" "mean_temp_cold"
## [13] "annual_percip" "percip_wet" "percip_dry"
## [16] "percip_seasonality" "percip_wet_quart" "percip_dry_quart"
## [19] "percip_warm_quart" "percip_cold_quart" "coords.x1"
## [22] "coords.x2"
table(model_df$occ)
##
## 0 1
## 873 91
head(model_df)
## occ annual_mean_temp mean_diurn_range isothermality temp_seasonality
## 1 1 253.32 75.00 345.96 3894.04
## 2 1 277.88 59.00 475.84 1213.92
## 3 1 181.08 74.00 279.72 6299.56
## 4 1 210.84 66.12 259.20 6446.28
## 5 1 216.76 67.00 288.80 5573.44
## 6 1 140.48 67.00 262.08 6361.80
## temp_max_warm temp_min_cold temp_annual_range mean_temp_wet mean_temp_dry
## 1 331.48 115.00 216.72 286.00 182.00
## 2 335.60 211.28 124.00 275.64 259.88
## 3 305.96 41.84 263.92 234.04 118.40
## 4 320.76 64.16 256.56 268.76 145.12
## 5 311.36 79.80 231.68 267.32 160.36
## 6 257.84 0.84 256.84 201.04 54.04
## mean_temp_warm mean_temp_cold annual_percip percip_wet percip_dry
## 1 287.00 182.00 2007.84 381.52 8.00
## 2 295.12 258.44 1368.48 255.80 1.00
## 3 263.64 92.44 1477.68 207.36 44.68
## 4 290.28 115.80 1871.72 350.00 33.72
## 5 283.20 133.84 1592.28 260.32 28.92
## 6 222.88 45.28 1143.80 222.12 20.12
## percip_seasonality percip_wet_quart percip_dry_quart percip_warm_quart
## 1 83.16 1116.40 24.00 979.16
## 2 82.76 744.12 10.84 129.80
## 3 44.80 617.08 141.52 488.24
## 4 70.76 1013.12 105.44 599.88
## 5 65.56 773.28 96.68 686.36
## 6 71.72 613.04 60.56 474.48
## percip_cold_quart coords.x1 coords.x2
## 1 24.00 90.33333 24.01667
## 2 73.72 104.04087 13.46905
## 3 162.60 118.28360 25.93977
## 4 136.76 113.08333 24.36667
## 5 113.56 112.43667 23.13833
## 6 61.44 106.59409 26.92268
variable_names <- names(env_asia_c)
cor_results <- cor(model_df[,variable_names])
corrplot::corrplot(cor_results, addCoef.col = "black", number.digits = 1, diag = FALSE)
## validating data
train_ids <- sample(1:nrow(model_df), size = 0.8*nrow(model_df))
training_data <- model_df[train_ids,]
testing_data <- model_df[-train_ids,]
nrow(training_data)
## [1] 771
nrow(testing_data)
## [1] 193
glm1_val <- glm(occ ~ annual_mean_temp+annual_percip, family=binomial(), data = training_data)
glm2_val <- glm(occ ~ mean_temp_cold+percip_dry, family=binomial(), data = training_data)
eT <- predicts::pa_evaluate(p = testing_data[testing_data$occ == 1,],
a = testing_data[testing_data$occ == 0,],
glm1_val,
type = "response")
eT
## @stats
## np na prevalence auc cor pcor ODP
## 1 17 176 0.088 0.836 0.337 0 0.912
##
## @thresholds
## max_kappa max_spec_sens no_omission equal_prevalence equal_sens_spec
## 1 0.144 0.116 0.032 0.088 0.149
##
## @tr_stats
## treshold kappa CCR TPR TNR FPR FNR PPP NPP MCR OR
## 1 0 0 0.09 1 0 1 0 0.09 NaN 0.91 NaN
## 2 0 0 0.09 1 0.01 0.99 0 0.09 1 0.91 Inf
## 3 0 0 0.09 1 0.01 0.99 0 0.09 1 0.91 Inf
## 4 ... ... ... ... ... ... ... ... ... ... ...
## 193 0.33 -0.01 0.91 0 0.99 0.01 1 0 0.91 0.09 0
## 194 0.33 -0.01 0.91 0 0.99 0.01 1 0 0.91 0.09 0
## 195 0.33 0 0.91 0 1 0 1 NaN 0.91 0.09 NaN
eP <- predicts::pa_evaluate(p = testing_data[testing_data$occ == 1,],
a = testing_data[testing_data$occ == 0,],
glm2_val,
type = "response")
eP
## @stats
## np na prevalence auc cor pcor ODP
## 1 17 176 0.088 0.846 0.374 0 0.912
##
## @thresholds
## max_kappa max_spec_sens no_omission equal_prevalence equal_sens_spec
## 1 0.363 0.11 0.053 0.088 0.134
##
## @tr_stats
## treshold kappa CCR TPR TNR FPR FNR PPP NPP MCR OR
## 1 0 0 0.09 1 0 1 0 0.09 NaN 0.91 NaN
## 2 0 0 0.09 1 0 1 0 0.09 NaN 0.91 NaN
## 3 0 0 0.09 1 0 1 0 0.09 NaN 0.91 NaN
## 4 ... ... ... ... ... ... ... ... ... ... ...
## 193 0.43 0.1 0.92 0.06 1 0 0.94 1 0.92 0.08 Inf
## 194 0.43 0 0.91 0 1 0 1 NaN 0.91 0.09 NaN
## 195 0.43 0 0.91 0 1 0 1 NaN 0.91 0.09 NaN
plot(eT, "ROC")
plot(eP, "ROC")
AIC(glm1_val) #
## [1] 403.3527
AIC(glm2_val) #
## [1] 386.1748
### estimate maxent asia
predictors <- c("percip_dry", "temp_min_cold")
maxent_args <-c('linear=true',
'quadratic=true',
'product=true',
'hinge=true',
'threshold=false',
'betamultiplier=1')
maxent_model <- MaxEnt(x = model_df[,predictors],
p = model_df$occ,
args = maxent_args)
predicts::partialResponse(maxent_model)
pred_maxent_asia <- terra::predict(env_asia_c, maxent_model)
plot(pred_maxent_asia)
train_ids <- sample(1:nrow(model_df), size = 0.8 * nrow(model_df))
training_data <- model_df[train_ids, ]
testing_data <- model_df[-train_ids, ]
e <- predicts::pa_evaluate(
p = testing_data[testing_data$occ == 1, ],
a = testing_data[testing_data$occ == 0, ],
maxent_model
)
plot(e, "ROC")
habitat_suitability <- predict(env_asia_c, glm2_val, type ="response")
mapview(habitat_suitability)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 5764225
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 5764225 '