---
title: "Script manuscrit"
author: "Bérénice AMBIELLE"
date: "2024-02-12"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

##Raster import 
```{r setup, include=FALSE}
e <- extent(741051.2, 808804, 7469297, 7523534)
p <- as(e, 'SpatialPolygons') 
h <-st_as_sf(p, coords=c("x","y"), crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
st_crs(h) <- "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"

##temperatures :
temperature <- raster("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/températures moyennes/moyenneRJ.tif")
temperature1 <- raster::projectRaster(temperature, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))

##precipitations :
precipitation<-raster("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/précipitations/moyprecipitations.tif")
precipitation1 <- raster::projectRaster(precipitation, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))

##vapor pressure :
vapeur<-raster("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vapeur/moyvapeur.tif")
vapeur1 <- raster::projectRaster(vapeur, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))

##radiation :
ra<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/radiation/radRJ.tif')
ra <- raster::projectRaster(ra, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))

##NDVI :
ndvi1<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/couverture/ndvi1.tif')
ndvi1 <- raster::projectRaster(ndvi1, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
ndvi1 <- projectRaster(ndvi1, precipitation1)

##population : 
pop<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/population/pop.tif')
pop <- raster::projectRaster(pop, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
pop <- projectRaster(pop, precipitation1)

##VCF Percent Non Tree Vegetation :
vcf_non_tree<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vcfok/PercentNonTreeVegetation_mean.tif')
vcf_non_tree<- raster::projectRaster(vcf_non_tree, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
vcf_non_tree <- projectRaster(vcf_non_tree, precipitation1)

##VCF Percent Non Vegetated :
vcf_non_vege<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vcfok/PercentNonVegetated_mean.tif')
vcf_non_vege<- raster::projectRaster(vcf_non_vege, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
vcf_non_vege <- projectRaster(vcf_non_vege, precipitation1)

##VCF Percent Tree Cover:
vcf_tree<-raster('C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vcfok/PercentTreeCover_mean.tif')
vcf_tree<- raster::projectRaster(vcf_tree, crs=CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))
vcf_tree<- projectRaster(vcf_tree, precipitation1)
```

###Changes to the study area :
```{r setup, include=FALSE}
temperature1<-crop(x = temperature1, y = h)
precipitation1 <- crop(x = precipitation1, y = h)
vapeur1<-crop(x=vapeur1, y=h)
ra1<-crop(x=ra, y=h)
ndvi2<-crop(x=ndvi1, y=h)
pop1<-crop(x=pop, y=h)
vcf_non_tree <- crop(x=vcf_non_tree, y=h)
vcf_non_vege <- crop(x=vcf_non_vege, y=h)
vcf_tree <- crop(x=vcf_tree, y=h)
land<-crop(x=land, y=h)
vege2<-crop(x=vege2, y=h)
```

###Multiple correlation test :
```{r setup, include=FALSE}
multicouche_sansfact <- stack(ndvi1, temperature1, precipitation1, vapeur1, ra1, pop1, vcf_non_vege, vcf_tree)
names(multicouche_sansfact)<- c("ndvi1", "moy_temp_ouis", "moy_prec_ouis", "moy_vap_ouis", "moy_ra_ouis", "moy_pop_ouis", "moy_vcf_non_vege_ouis", "moy_vcf_tree_ouis")

jnk=layerStats(multicouche_sansfact, 'pearson', na.rm=T)
corr_matrix=jnk$'pearson correlation coefficient'
corr_matrix
corrplot(corr_matrix, type = "upper")
```

### Vector layer import, rasterization, resolution : 
```{r setup, include=FALSE}
##vegetation :
vegetation<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vege", 'vegeRJ')


vegetation<-st_transform(vegetation, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
vegetation <- as_Spatial(vegetation)
raster_vegetation <- rasterize(vegetation, precipitation1, "leg_carga")
vege2<-raster_vegetation


##land use :
land_use<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/occupation de l'espace", 'Rio_de_Janeiro')
land_use<-st_transform(land_use, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
raster_land_use<-raster("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/occupation de l'espace/raster_land_use.tif")
raster_land_use <- projectRaster(raster_land_use, precipitation1, method = 'ngb')
land<-raster_land_use
```
###Roads
```{r setup, include=FALSE}
routeshap<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin", 'route')
route<-routeshap
route$geometry[st_is_empty(route$geometry)] <- NA
route$highway <- as.factor(route$highway)

niveaux <- c("motorway", "primary", "secondary", "tertiary","trunk")
shapefile_sub <- subset(route, highway == niveaux)
```

###Importing presence data :
```{r setup, include=FALSE}
data<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/données", 'master_pa')
data<-st_transform(data, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
data<-erase.point(data, h, inside = FALSE)
```

#Callithrix
```{r setup, include=FALSE}
data1<-data
data1<-data_article
ouis<-data1$Cllthrx

Cllthrx <- data1[ouis == "1",]
absence_ouis<-data1[data1$Cllthrx == "0",]

absence_ouis<-st_as_sf(absence_ouis, coords=data1$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
Cllthrx<-st_as_sf(Cllthrx, coords=data1$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

tampons_presence_ouis<- st_buffer(Cllthrx, dist = 177.29)


data1[5]<-data1[5][!is.na(data1[5]),]
points_non_contenus <- data1[5]
for(i in 1:nrow(tampons_presence_ouis)) {
  points_non_contenus <- points_non_contenus[!st_intersects(points_non_contenus, tampons_presence_ouis[i,], sparse = FALSE),]
}


data2 <- data.frame()
for(file in c("points_non_contenus", "Cllthrx")){
  temp_file <- get(file)
  data2 <- rbind(data2, temp_file[c("geometry", "Cllthrx")])
}


data2<-st_as_sf(data2, coords=data2$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
data2<-st_transform(data2, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
empty_geom <- st_is_empty(data2)
data_sf_nonempty <- data2[!empty_geom, ]
data2<-data_sf_nonempty
data2<-erase.point(data2, h, inside = FALSE)



data2_1<- st_as_sfc(Cllthrx)
empty_geoms <- st_is_empty(data2_1)
data2_1 <- st_drop_geometry(data2_1[!empty_geoms])
data_sp_1 <-as(data2_1, "Spatial")



cells <- cellFromXY(temperature1, data_sp_1)
dups <- duplicated(cells)
data_pixel_1 <- data_sp_1[!dups, ]
if (sum(dups) > 250) {
  dups_idx <- which(dups)
  data_pixel_1 <- rbind(data_pixel_1, data_sp_1[dups_idx[1:250],])
}
cat(nrow(data_sp_1) - nrow(data_pixel_1), "records are removed")


data_pixel_sf_obj <- as(data_pixel_1, "sf")
data_pixel_1<-st_as_sf(data_pixel_sf_obj, coords=data_pixel_sf_obj$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs") 


data2_0<- st_as_sfc(points_non_contenus)
empty_geoms <- st_is_empty(data2_0)
data2_0 <- st_drop_geometry(data2_0[!empty_geoms])
data_sp_0 <-as(data2_0, "Spatial")

cells <- cellFromXY(temperature1, data_sp_0)
dups <- duplicated(cells)
data_pixel_0 <- data_sp_0[!dups, ]
cat(nrow(data_sp_0) - nrow(data_pixel_0), "records are removed")


data_pixel_sf_obj <- as(data_pixel_0, "sf")
data_pixel_0<-st_as_sf(data_pixel_sf_obj, coords=data_pixel_sf_obj$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


absent<-rep(0,907)
data_pixel_0<-data.frame(data_pixel_0, absent)
names(data_pixel_0)<-c("geometry","Cllthrx")
data_pixel_0<-st_as_sf(data_pixel_0, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


pres<-rep(1,123)
data_pixel_1<-data.frame(data_pixel_1, pres)
names(data_pixel_1)<-c("geometry","Cllthrx")
data_pixel_1<-st_as_sf(data_pixel_1, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


data2<- data.frame()
for(file in c("data_pixel_0", "data_pixel_1")){
  temp_file <- get(file)
  data2 <- rbind(data2, temp_file[c( "Cllthrx","geometry")])
}
```

### Absence
```{r setup, include=FALSE}
ouis<-data2$Cllthrx
absence_ouis<-data2[ouis == "0",]

absence_ouis1<-st_cast(st_sample(absence_ouis, size=373), "POINT")
absence_ouis_rest<-st_cast(st_sample(absence_ouis1, size=300), "POINT")
st_crs(absence_ouis1)<-"+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
st_crs(absence_ouis_rest)<-"+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
absence_ouis1<-st_as_sf(absence_ouis1, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
absence_ouis_rest<-st_as_sf(absence_ouis_rest, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

 
absence_ouis1_test <- absence_ouis1[-which(absence_ouis1$x %in% absence_ouis_rest$x),]
absence_ouis1_test<-st_cast(st_sample(absence_ouis1_test, size=73), "POINT")

ab<-rep(0,300)
absence_ouis2<-data.frame(absence_ouis_rest, ab)
names(absence_ouis2)<-c("geometry", "Cllthrx")
absence_ouis2<-st_as_sf(absence_ouis2, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

abc<-rep(0,73)
absence_ouis2_test<-data.frame(absence_ouis1_test, abc)
names(absence_ouis2_test)<-c("geometry", "Cllthrx")
absence_ouis2_test<-st_as_sf(absence_ouis2_test, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
```


###For the test data subset
```{r setup, include=FALSE}
tampon_absence_ouis2_test <- st_buffer(absence_ouis2_test, dist = 177.29)

raster_absence_ouis2_test <- raster::extract(x = Nom_raster_multlicouche, 
                                y = tampon_absence_ouis2_test , 
                                df = TRUE)


moy_temp_ouis<-tapply(raster_absence_ouis2_test$moy_temp_ouis, raster_absence_ouis2_test$ID, mean)
moy_vap_ouis<-tapply(raster_absence_ouis2_test$moy_vap_ouis, raster_absence_ouis2_test$ID, mean)
moy_ra_ouis<-tapply(raster_absence_ouis2_test$moy_ra_ouis, raster_absence_ouis2_test$ID, mean)
moy_ndvi_ouis<-tapply(raster_absence_ouis2_test$ndvi1, raster_absence_ouis2_test$ID, mean)
moy_pop_ouis<-tapply(raster_absence_ouis2_test$moy_pop_ouis, raster_absence_ouis2_test$ID, mean)
moy_land_ouis<-extract(land, absence_ouis2_test)
moy_vege2_ouis<-extract(vege2, absence_ouis2_test)
moy_vcf_non_vege_ouis<-tapply(raster_absence_ouis2_test$moy_vcf_non_vege_ouis, raster_absence_ouis2_test$ID, mean)
moy_vcf_tree_ouis<-tapply(raster_absence_ouis2_test$moy_vcf_tree_ouis, raster_absence_ouis2_test$ID, mean)
moy_prec_ouis<-tapply(raster_absence_ouis2_test$moy_prec_ouis, raster_absence_ouis2_test$ID, mean)

data_merged_ouis_abs_test <- data.frame(moy_temp_ouis,moy_vap_ouis, moy_ra_ouis,moy_ndvi_ouis, moy_pop_ouis, moy_vcf_non_vege_ouis, moy_vcf_tree_ouis,moy_land_ouis,moy_vege2_ouis, moy_prec_ouis, absence_ouis2_test)
```



###For the main dataset
```{r setup, include=FALSE}
tampons_absence_ouis <- st_buffer(absence_ouis2, dist = 177.29)

temp_abs_ouis <- raster::extract(x = Nom_raster_multlicouche, 
                                y = tampons_absence_ouis , 
                                df = TRUE)


moy_temp_ouis<-tapply(temp_abs_ouis$moy_temp_ouis, temp_abs_ouis$ID, mean)
moy_vap_ouis<-tapply(temp_abs_ouis$moy_vap_ouis, temp_abs_ouis$ID, mean)
moy_ra_ouis<-tapply(temp_abs_ouis$moy_ra_ouis, temp_abs_ouis$ID, mean)
moy_ndvi_ouis<-tapply(temp_abs_ouis$ndvi1, temp_abs_ouis$ID, mean)
moy_pop_ouis<-tapply(temp_abs_ouis$moy_pop_ouis, temp_abs_ouis$ID, mean)
moy_land_ouis<-extract(land, absence_ouis2)
moy_vege2_ouis<-extract(vege2, absence_ouis2)
moy_vcf_non_vege_ouis<-tapply(temp_abs_ouis$moy_vcf_non_vege_ouis, temp_abs_ouis$ID, mean)
moy_vcf_tree_ouis<-tapply(temp_abs_ouis$moy_vcf_tree_ouis, temp_abs_ouis$ID, mean)
moy_prec_ouis<-tapply(temp_abs_ouis$moy_prec_ouis, temp_abs_ouis$ID, mean)

absence_ouis2_sfc<- st_as_sfc(absence_ouis2)
shapefile_sub_sfc<- st_as_sfc(shapefile_sub[6])
distance<-st_distance(absence_ouis2_sfc, shapefile_sub_sfc)
min_distances<- apply(distance, 1, min)


data_merged_ouis_abs_dist <- data.frame(moy_temp_ouis,moy_vap_ouis, moy_ra_ouis,moy_ndvi_ouis, moy_pop_ouis, moy_vcf_non_vege_ouis, moy_vcf_tree_ouis,moy_land_ouis,moy_vege2_ouis, moy_prec_ouis,min_distances, absence_ouis2)
crs <- "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
data_merged_ouis_abs_dist <- st_set_crs(data_merged_ouis_abs_dist, crs)

```

###Presences
```{r setup, include=FALSE}
ouis<-data2$Cllthrx
presence_ouis <- data2[ouis == "1",]

presence_ouis_rest<-st_cast(st_sample(presence_ouis, size=300), "POINT")
st_crs(presence_ouis_rest)<-"+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
presence_ouis<-st_as_sf(presence_ouis, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
presence_ouis_rest<-st_as_sf(presence_ouis_rest, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


presence_ouis1_test <- presence_ouis[-which(presence_ouis$geometry %in% presence_ouis_rest$x),]
presence_ouis1_test <- presence_ouis[!(presence_ouis$geometry %in% presence_ouis_rest$x), ]


ab<-rep(1,300)
presence_ouis2<-data.frame(presence_ouis_rest, ab)
names(presence_ouis2)<-c("geometry", "Cllthrx")
presence_ouis2<-st_as_sf(presence_ouis2, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
presence_ouis1_test<-st_as_sf(presence_ouis1_test, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
```


###For the test data subset
```{r}
tampon_presence_ouis1_test <- st_buffer(presence_ouis1_test , dist = 177.29)

raster_presence_ouis1_test <- raster::extract(x = Nom_raster_multlicouche, 
                                y = tampon_presence_ouis1_test , 
                                df = TRUE)


moy_temp_ouis<-tapply(raster_presence_ouis1_test$moy_temp_ouis, raster_presence_ouis1_test$ID, mean)
moy_vap_ouis<-tapply(raster_presence_ouis1_test$moy_vap_ouis, raster_presence_ouis1_test$ID, mean)
moy_ra_ouis<-tapply(raster_presence_ouis1_test$moy_ra_ouis, raster_presence_ouis1_test$ID, mean)
moy_ndvi_ouis<-tapply(raster_presence_ouis1_test$ndvi1, raster_presence_ouis1_test$ID, mean)
moy_pop_ouis<-tapply(raster_presence_ouis1_test$moy_pop_ouis, raster_presence_ouis1_test$ID, mean)
moy_land_ouis<-extract(land, presence_ouis1_test)
moy_vege2_ouis<-extract(vege2, presence_ouis1_test)
moy_vcf_non_vege_ouis<-tapply(raster_presence_ouis1_test$moy_vcf_non_vege_ouis, raster_presence_ouis1_test$ID, mean)
moy_vcf_tree_ouis<-tapply(raster_presence_ouis1_test$moy_vcf_tree_ouis, raster_presence_ouis1_test$ID, mean)
moy_prec_ouis<-tapply(raster_presence_ouis1_test$moy_prec_ouis, raster_presence_ouis1_test$ID, mean)




data_merged_ouis_pres_test <- data.frame(moy_temp_ouis,moy_vap_ouis, moy_ra_ouis,moy_ndvi_ouis, moy_pop_ouis, moy_vcf_non_vege_ouis, moy_vcf_tree_ouis,moy_land_ouis,moy_vege2_ouis, moy_prec_ouis, presence_ouis1_test)
```


###For the main dataset 
```{r setup, include=FALSE}
tampons_presence_ouis <- st_buffer(presence_ouis2, dist = 177.29)

temp_pres_ouis <- raster::extract(x = Nom_raster_multlicouche, 
                                y = tampons_presence_ouis  , 
                                df = TRUE)

moy_temp_ouis<-tapply(temp_pres_ouis$moy_temp_ouis, temp_pres_ouis$ID, mean)
moy_vap_ouis<-tapply(temp_pres_ouis$moy_vap_ouis, temp_pres_ouis$ID, mean)
moy_ra_ouis<-tapply(temp_pres_ouis$moy_ra_ouis, temp_pres_ouis$ID, mean)
moy_ndvi_ouis<-tapply(temp_pres_ouis$ndvi1, temp_pres_ouis$ID, mean)
moy_pop_ouis<-tapply(temp_pres_ouis$moy_pop_ouis, temp_pres_ouis$ID, mean)
moy_land_ouis<-extract(land, presence_ouis2)
moy_vege2_ouis<-extract(vege2, presence_ouis2)
moy_vcf_non_vege_ouis<-tapply(temp_pres_ouis$moy_vcf_non_vege_ouis, temp_pres_ouis$ID, mean)
moy_vcf_tree_ouis<-tapply(temp_pres_ouis$moy_vcf_tree_ouis, temp_pres_ouis$ID, mean)
moy_prec_ouis<-tapply(temp_pres_ouis$moy_prec_ouis, temp_pres_ouis$ID, mean)


presence_ouis2_sfc<- st_as_sfc(presence_ouis2)
shapefile_sub_sfc<- st_as_sfc(shapefile_sub[6])
distance<-st_distance(presence_ouis2_sfc, shapefile_sub_sfc)
min_distances<- apply(distance, 1, min)
Cllthrxx<-rep(1,300)

data_merged_ouis_pres <- data.frame(moy_temp_ouis_d,moy_vap_ouis_d, moy_ra_ouis_d,moy_pop_ouis_d, moy_ndvi_ouis_d, moy_vcf_non_vege_ouis_d, moy_vcf_tree_ouis_d,moy_land_ouis_d,moy_vege2_ouis_d, moy_prec_ouis_d, presence_ouis_d2)

data_merged_ouis_pres_dist <- data.frame(moy_temp_ouis,moy_vap_ouis, moy_ra_ouis,moy_pop_ouis, moy_ndvi_ouis, moy_vcf_non_vege_ouis, moy_vcf_tree_ouis,moy_land_ouis,moy_vege2_ouis, moy_prec_ouis, min_distances, presence_ouis2)
data_merged_ouis_pres_dist <- st_as_sf(data_merged_ouis_pres_dist, coords = c("coords.x1", "coords.x2"))
crs <- "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
data_merged_ouis_pres_dist <- st_set_crs(data_merged_ouis_pres_dist, crs)
```

###Tables
```{r setup, include=FALSE}
data_complet_merged_ouis_dist <- data.frame()
for(file in c("data_merged_ouis_pres_dist", "data_merged_ouis_abs_dist")){
  temp_file <- get(file)
  data_complet_merged_ouis_dist <- rbind(data_complet_merged_ouis_dist, temp_file[c("geometry","moy_temp_ouis", "moy_ra_ouis", "moy_pop_ouis", "moy_ndvi_ouis", "moy_vap_ouis", "moy_vcf_non_vege_ouis", "moy_vcf_tree_ouis" , "moy_vege2_ouis", "Cllthrx","moy_land_ouis", "moy_prec_ouis","min_distances")])
}


data_ouis_test <- data.frame()
for(file in c("data_merged_ouis_pres_test", "data_merged_ouis_abs_test")){
  temp_file <- get(file)
  data_ouis_test <- rbind(data_ouis_test , temp_file[c("geometry","moy_temp_ouis", "moy_vap_ouis", "moy_ra_ouis", "moy_pop_ouis", "moy_ndvi_ouis", "moy_vcf_non_vege_ouis", "moy_vcf_tree_ouis" , "moy_vege2_ouis", "moy_land_ouis", "moy_prec_ouis", "Cllthrx")])
}

absence_ouistest_sfc<- st_as_sfc(data_ouis_test)
shapefile_sub_sfc<- st_as_sfc(shapefile_sub[6])
distance<-st_distance(absence_ouistest_sfc, shapefile_sub_sfc)
min_distances<- apply(distance, 1, min)
data_ouis_test$min_distances<-min_distances


data_ouis_test_model <- data_ouis_test[-which(data_ouis_test$moy_vege2_ouis == 260),]
data_ouis_test_2.0 <- data.frame()
for(file in c("data_ouis_test_model")){
  temp_file <- get(file)
  data_ouis_test_2.0<- rbind(data_ouis_test_2.0, temp_file[c("geometry", "moy_vcf_tree_ouis", "moy_vcf_non_vege_ouis", "moy_vege2_ouis" , "moy_vap_ouis", "Cllthrx")])
}

donnees_shapefile_sans_na <- na.omit(data_ouis_test_2.0)
levels(donnees_shapefile_sans_na$Cllthrx) <- levels(data_complet_merged_ouis_dist$Cllthrx)
levels(donnees_shapefile_sans_na$moy_vege2_ouis) <- levels(data_complet_merged_ouis_dist$moy_vege2_ouis)
selection_0<-donnees_shapefile_sans_na %>%
  filter(Cllthrx ==0)%>%
  sample_n(60)
selection_1<-donnees_shapefile_sans_na %>%
  filter(Cllthrx ==1)
donnees_shapefile_sans_na<-rbind(selection_1, selection_0)

```

###ACP
```{r}
library(ade4)
ACP_ouis <- dudi.pca(data_complet_merged_ouis_dist[5:7], center = TRUE, scale = TRUE, scannf = FALSE)
summary(ACP_ouis)

s.corcircle(ACP_ouis$co)
dudi.mix
```

#GLM
```{r}
data_complet_merged_ouis_dist$moy_vege2_ouis<-as.factor(data_complet_merged_ouis_dist$moy_vege2_ouis)
data_complet_merged_ouis_dist$moy_land_ouis<-as.factor(data_complet_merged_ouis_dist$moy_land_ouis)
data_complet_merged_ouis_dist$Cllthrx<-as.factor(data_complet_merged_ouis_dist$Cllthrx)
data_complet_merged_ouis_dist$moy_temp_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_temp_ouis)
data_complet_merged_ouis_dist$moy_vap_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_vap_ouis)
data_complet_merged_ouis_dist$moy_ra_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_ra_ouis)
data_complet_merged_ouis_dist$moy_ndvi_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_ndvi_oui)
data_complet_merged_ouis_dist$moy_pop_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_pop_ouis)
data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis)
data_complet_merged_ouis_dist$moy_vcf_tree_ouis<-as.numeric(data_complet_merged_ouis_dist$moy_vcf_tree_ouis)
data_complet_merged_ouis_dist$min_distances<-as.numeric(data_complet_merged_ouis_dist$min_distances)


null.m <- glm(Cllthrx ~ 1, family = binomial(link = "logit"), data =
                data_complet_merged_ouis_dist)

Nom_modele_step_dist <- stepAIC(null.m, Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, family=binomial(link = "logit"), data = data_complet_merged_ouis_dist)


model.final.dist <- glm(Cllthrx ~ moy_vap_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + 
    moy_vcf_tree_ouis, family=binomial(link = "logit"), data=data_complet_merged_ouis_dist)
```

####GLM model predictions and probability projection
```{r}
resolution_reference <- res(temperature1)
extent_reference <- extent(temperature1)
nrow <- nrow(temperature1)
ncol <- ncol(temperature1)
nouveau_raster <- raster(nrows = nrow, ncols = ncol, ext = extent_reference, crs = crs(temperature1))
centroids <- raster::cellFromXY(nouveau_raster, xy = coordinates(nouveau_raster))
points <- SpatialPoints(coordinates(nouveau_raster)[centroids,], proj4string = crs(nouveau_raster))
points_df <- SpatialPointsDataFrame(points, data = data.frame(ID = 1:length(points)), match.ID = FALSE)

points_sf <- st_as_sf(points_df)
shapefile_sub_sf <- st_as_sf(shapefile_sub)
distances <- st_distance(points_sf, shapefile_sub_sf)
min_distances_ras<- apply(distances, 1, min)


points_df$min_distance <- min_distances_ras
library(sf)

points_sf <- st_as_sf(points_df)
min_distance_data <- st_transform(points_sf, crs = "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
points_raster <- rasterize(min_distance_data, nouveau_raster, field = "min_distance")


model.final.dist <- glm(Cllthrx ~ moy_vap_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + 
    moy_vcf_tree_ouis, family=binomial(link = "logit"), data=data_complet_merged_ouis_dist)

prediction <- function(model.final.dist,data_complet_merged_ouis_dist, sefit=FALSE) {
  v <-stats::predict(model.final.dist, data_complet_merged_ouis_dist, se.fit=T, type="response")
  if(sefit==FALSE){p=as.vector(v$fit)} else {p=as.vector(v$se.fit)}
  p
}

points_raster<-crop(x = points_raster, y = precipitation1)
raster_stack <- stack(ndvi2, temperature1, precipitation1, vapeur1, ra1, pop1, land, vege2, vcf_non_vege, vcf_tree, points_raster)
names(raster_stack )<- c("moy_ndvi_ouis", "moy_temp_ouis", "moy_prec_ouis", "moy_vap_ouis", "moy_ra_ouis", "moy_pop_ouis", "moy_land_ouis", "moy_vege2_ouis", "moy_vcf_non_vege_ouis", "moy_vcf_tree_ouis", "min_distances")

carte_prediction<-predict(raster_stack,
                          model=model.final.dist, fun=prediction,
                          index=11)

carte_prediction[is.na(carte_prediction[])] <- 0


carte_prediction_sd<-predict(raster_stack,
model=model.final.dist, fun=prediction, sefit=TRUE,
index=1)
carte_prediction_sd[is.na(carte_prediction_sd[])] <- 0
```

###KAPPA
```{r}
prediction<-predict(model.final.dist, donnees_shapefile_sans_na)
prediction <- ifelse(prediction > 0.5,1,0 )
confusion_matrix <- table(Predicted = predictions, Actual = donnees_shapefile_sans_na$Cllthrx)
a <- confusion_matrix[1,1] 
b <- confusion_matrix[1,2] 
c <- confusion_matrix[2,1] 
d <- confusion_matrix[2,2] 

n <- sum(confusion_matrix)
po <- (a + d) / n 
pe <- ((a + b) * (a + c) + (c + d) * (b + d)) / (n^2) 
kappa <- (po - pe) / (1 - pe)

print(kappa)

```
###Matthews Index  (MCC)
```{r}
MCC <- (a * d - b * c) / sqrt((a + b) * (a + c) * (d + b) * (d + c))#0.748
```

###Model coefficient
```{r setup, include=FALSE}
pseudoR2 <- (model.final.dist$null.deviance - model.final.dist$deviance) / model.final.dist$null.deviance
pseudoR2
exp(confint(model.final.dist))

library(questionr)
odds.ratio(model.final.dist)
```


###Graphics
```{r}
library(GGally)
plot <- ggcoef_model(model.final.dist, exponentiate = TRUE, variable_labels = c(moy_vap_ouis = "Water vapor pressure", moy_vege2_ouis = "Type of vegetation", moy_vcf_non_vege_ouis = "Percentage of non-vegetated surface", moy_vcf_tree_ouis = "Percentage of tree cover"))

plot +
    theme(
        axis.title.x = element_text(size = 14),  
        axis.title.y = element_text(size = 14),
        axis.text.x = element_text(size = 14), 
        axis.text.y = element_text(size = 14)
    )
```

###Water vapor pressure 
```{r}
data_complet_merged_ouis_dist$moy_vap_ouis1<-data_complet_merged_ouis_dist$moy_vap_ouis/100
da<-data.frame(data_complet_merged_ouis_dist$moy_vap_ouis1, data_complet_merged_ouis_dist$Cllthrx)
data_complet_merged_ouis_dist$moy_vap_ouis1<-as.numeric(data_complet_merged_ouis_dist$moy_vap_ouis1)
means_vap_ouis<- c(mean(data_complet_merged_ouis_dist$moy_vap_ouis1[data_complet_merged_ouis_dist$Cllthrx == "0"]), mean(data_complet_merged_ouis_dist$moy_vap_ouis1[data_complet_merged_ouis_dist$Cllthrx == "1"]))
sd_vap_ouis <- c(sd(data_complet_merged_ouis_dist$moy_vap_ouis1[data_complet_merged_ouis_dist$Cllthrx == "0"]), sd(data_complet_merged_ouis_dist$moy_vap_ouis1[data_complet_merged_ouis_dist$Cllthrx == "1"]))


plot_vap <- ggplot(data_complet_merged_ouis_dist, aes(x=moy_vap_ouis1, fill = Cllthrx)) +
  geom_histogram(position="identity", alpha=0.5, color="black")+
  geom_vline(xintercept = means_vap_ouis, color=c("black", "green"), linetype = "dashed",size = 1.25)+
  scale_color_manual(values=c("#999999", "darkgreen"))+
  scale_fill_manual(values=c("#999999", "darkgreen"))+
  labs(title="Distribution of water vapour pressure variable at marmoset presence/absence points",x="Vapour pressure", y = "Frequency", fill="Presence/Absence")+
  theme_classic()
```
###VCF non vegetation
```{r}
means_vcf_non_vege <- c(mean(data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis[data_complet_merged_ouis_dist$Cllthrx == "0"]), mean(data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis[data_complet_merged_ouis_dist$Cllthrx == "1"]))
sd_vcf_non_vege <- c(sd(data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis[data_complet_merged_ouis_dist$Cllthrx == "0"]), sd(data_complet_merged_ouis_dist$moy_vcf_non_vege_ouis[data_complet_merged_ouis_dist$Cllthrx == "1"]))

plot_vcf_non_vege <- ggplot(data_complet_merged_ouis_dist, aes(x=moy_vcf_non_vege_ouis, fill = Cllthrx)) +
  geom_histogram(position="identity", alpha=0.5, color="black")+
  geom_vline(xintercept = means_vcf_non_vege, color=c("black", "green"), linetype = "dashed",size = 1.25)+
  scale_color_manual(values=c("#999999", "darkgreen"))+
  scale_fill_manual(values=c("#999999", "darkgreen"))+
  labs(title="Distribution of percentage of non-vegetated cover variable at marmoset presence/absence points",x="percentage of non-vegetated cover", y = "Frequency", fill="Presence/Absence")+
  theme_classic()
```
###VCF tree
```{r}
means_vcf_tree <- c(mean(data_complet_merged_ouis_dist$moy_vcf_tree_ouis[data_complet_merged_ouis_dist$Cllthrx == "0"]), mean(data_complet_merged_ouis_dist$moy_vcf_tree_ouis[data_complet_merged_ouis_dist$Cllthrx == "1"]))
sd_vcf_tree <- c(sd(data_complet_merged_ouis_dist$moy_vcf_tree_ouis[data_complet_merged_ouis_dist$Cllthrx == "0"]), sd(data_complet_merged_ouis_dist$moy_vcf_tree_ouis[data_complet_merged_ouis_dist$Cllthrx == "1"]))

plot_vcf_tree <- ggplot(data_complet_merged_ouis_dist, aes(x=moy_vcf_tree_ouis, fill = Cllthrx)) +
  geom_histogram(position="identity", alpha=0.5, color="black")+
  geom_vline(xintercept = means_vcf_tree, color=c("black", "green"), linetype = "dashed",size = 1.25)+
  scale_color_manual(values=c("#999999", "darkgreen"))+
  scale_fill_manual(values=c("#999999", "darkgreen"))+
  labs(title="Distribution of percentage of tree cover variable at marmoset presence/absence points",x="Percentage of tree cover", y = "Frequency", fill="Presence/Absence")+
  theme_classic()
```

###Vegetation
```{r}
plot_vege <- ggplot(data_complet_merged_ouis_dist, aes(x = moy_vege2_ouis, fill = Cllthrx)) + 
  geom_bar(position = "dodge") + 
  scale_fill_manual(values = c("#999999", "darkgreen")) +
  labs(x = "Type of vegetation", y = "Frequency", fill = "Presence/Absence") +
  ggtitle("Distribution of type of vegetation at marmoset presence/absence points") +
  theme_classic()
plot(plot_vege)
```
###Statistical test
```{r}
vegetation<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vege", 'vegeRJ')
vegetation<-st_transform(vegetation, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
vegetation$leg_carga<-as.factor(vegetation$leg_carga)
levels(vegetation$leg_carga)[data_complet_merged_ouis_dist$moy_vege2_ouis]
cbind(levels(vegetation$leg_carga),seq(1,length(levels(vegetation$leg_carga))))


vegetation$leg_sup<-as.factor(vegetation$leg_sup)
raster_vegetation_anth <- rasterize(vegetation, precipitation1, "leg_sup")
ouis<-data_complet_merged_ouis_dist$Cllthrx
presence_ouis2<-data_complet_merged_ouis_dist[ouis == "1",]

presence_ouis2<- st_as_sfc(presence_ouis2)
presence_ouis2 <-as(presence_ouis2, "Spatial")
presence<-extract(raster_vegetation_anth, presence_ouis2)
ab<-rep(1,300)
presence_anth<-data.frame(presence, ab)


presence_anth$presence<-as.factor(presence_anth$presence)
presence_anth$ab<-as.factor(presence_anth$ab)

plot_anth <- ggplot(presence_anth, aes(x = presence, fill = presence)) + 
  geom_bar(position = "dodge") + 
  scale_fill_manual(values = c("lightgreen", "darkgreen"),
                    labels = c("Dominant anthropic area", "Dominant natural vegetation")) +
  labs(x = "Type of vegetation", y = "Frequency") +
  ggtitle("Distribution of marmosets according to vegetation type: anthropic or natural") +
  theme_classic() +
  guides(fill = guide_legend(title = NULL)) + 
  theme(legend.position = "right")  



presence_anth14 <- subset(presence_anth, presence %in% c("1", "4"))
contingence14 <- table(presence_anth14$ab, presence_anth14$presence)
j<-chisq.test(contingence14)
j$stdres
```

###ANOVA
```{r}
anova(model.final.dist,test="Chisq")
par_var_expl_lmodel.final <- anova(model.final.dist)[-
1,2]*100 / anova(model.final.dist)[1,4]
print( cbind(rownames(anova(model.final.dist))[-
1],round(par_var_expl_lmodel.final  , digit = 1)))
```

###Probability graph
```{r}
courbe_vcf_non_vege <- ggeffect(model.final.dist, terms = c("moy_vcf_non_vege_ouis"))
a <- ggplot() +
    geom_line(data = courbe_vcf_non_vege, aes(x = x, y = predicted), color = "black", linewidth = 1) +  
    labs(x = "Percent of non-vegetated cover", y = "Predicted probability of occurrence") +
    theme(
        panel.background = element_rect(fill = "white", color = "white"),
        plot.background = element_rect(fill = "white", color = "white"),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)
    )


courbe_vcf_tree <- ggeffect(model.final.dist, terms = c("moy_vcf_tree_ouis"))
b <- ggplot() +
    geom_line(data = courbe_vcf_tree, aes(x = x, y = predicted), color = "black", linewidth = 1) +  
    labs(x = "Pourcent of tree cover", y = "Predicted probability of occurrence") +
    theme(
        panel.background = element_rect(fill = "white", color = "white"),
        plot.background = element_rect(fill = "white", color = "white"),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)
    )

model.final.dist.vap<- glm(Cllthrx ~ moy_vap_ouis1 + moy_vege2_ouis + moy_vcf_non_vege_ouis + 
                          moy_vcf_tree_ouis, family=binomial(link = "logit"), data=data_complet_merged_ouis_dist)

courbe_vap <- ggeffect(model.final.dist.vap, terms = c("moy_vap_ouis1"))
c <- ggplot() +
    geom_line(data = courbe_vap, aes(x = x, y = predicted), color = "black", linewidth = 1) +  
    labs(x = "Water vapor pressure", y = "Predicted probability of occurrence") +
    theme(
        panel.background = element_rect(fill = "white", color = "white"),
        plot.background = element_rect(fill = "white", color = "white"),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)
    )

```

###Area under the curve
```{r}
library(pROC)
library(caret)
pred_test <- predict(model.final.dist, data_complet_merged_ouis_dist, type="response")
test_roc = roc(data_complet_merged_ouis_dist$Cllthrx ~ pred_test, plot = TRUE, print.auc = TRUE)
plot(test_roc,auc.polygon=TRUE)
auc(test_roc)
```

###Cross validation
```{r}
K = 5
nk = floor(nrow(data_complet_merged_ouis_dist)/K)
pred = NULL
Yobs = NULL
rk=list()
set.seed(48)
ii = sample(nrow(data_complet_merged_ouis_dist),nrow(data_complet_merged_ouis_dist)) # permutation aléatoire des indices
for (k in 1:K){
  testi = ii[((k-1)*nk+1):(k*nk)]
  appri = setdiff(1:nrow(data_complet_merged_ouis_dist),testi)
  blm<-glm(Cllthrx ~ moy_vap_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + 
    moy_vcf_tree_ouis, family = binomial(link = "logit"), data = data_complet_merged_ouis_dist)
           
  pr = predict(blm,data_complet_merged_ouis_dist[testi,],type="response")
  pred = c(pred,pr)
  Yobs = c(Yobs,data_complet_merged_ouis_dist$Cllthrx[testi])
  rk[[k]]=roc(data_complet_merged_ouis_dist$Cllthrx[testi],pr)
}
add = FALSE
for (k in 1:K){
  if (k>1) add=TRUE
  plot.roc(rk[[k]],col="gray",add=add,lwd=1)
}
r=roc(Yobs,pred)
plot(r,main=paste("AUROC =",round(r$auc*100)/100),add=TRUE)
grid()
auc(r)
auc(rk[[k]])
```
###Cross validation on the test dataset
```{r}
K = 5
nk = floor(nrow(donnees_shapefile_sans_na)/K)
pred = NULL
Yobs = NULL
rk=list()
set.seed(48)
ii = sample(nrow(donnees_shapefile_sans_na),nrow(donnees_shapefile_sans_na)) # permutation aléatoire des indices
for (k in 1:K){
  testi = ii[((k-1)*nk+1):(k*nk)]
  appri = setdiff(1:nrow(data_complet_merged_ouis_dist),testi)
  blm<-glm(Cllthrx ~ moy_vap_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + 
    moy_vcf_tree_ouis, family = binomial(link = "logit"), data = data_complet_merged_ouis_dist)
           
  pr = predict(blm,donnees_shapefile_sans_na[testi,],type="response")
  pred = c(pred,pr)
  Yobs = c(Yobs,donnees_shapefile_sans_na$Cllthrx[testi])
  rk[[k]]=roc(donnees_shapefile_sans_na$Cllthrx[testi],pr)
}
add = FALSE
for (k in 1:K){
  if (k>1) add=TRUE
  plot.roc(rk[[k]],col="gray",add=add,lwd=1)
}
r=roc(Yobs,pred)
plot(r,main=paste("AUROC =",round(r$auc*100)/100),add=TRUE)
grid()
auc(r)
auc(rk[[k]])

```

###PseudoR²
```{r}
model_summary <- summary(model.final.dist)

r_squared <- model_summary$deviance / model_summary$null.deviance

print(r_squared)
```

###Accuracy
```{r}
donnees_shapefile_sans_na$moy_vege2_ouis<-as.factor(donnees_shapefile_sans_na$moy_vege2_ouis)
predictions_acc <- predict(model.final.dist, donnees_shapefile_sans_na, type = "response")
predictions_acc <- ifelse(predictions_acc > 0.5,1,0 ) 
accuracy <- sum(predictions_acc == donnees_shapefile_sans_na$Cllthrx) / length(donnees_shapefile_sans_na$Cllthrx)
print(paste("Accuracy:", accuracy)) 
```

###Sensitivity, specificity, F1-score validation data
```{r}
new_predictions <- predictions_acc
new_confusion_matrix <- table(donnees_shapefile_sans_na$Cllthrx, new_predictions)
new_sensitivity <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[2, ])
print(paste("Sensitivity (Recall) on new data:", new_sensitivity))
new_specificity <- new_confusion_matrix[1, 1] / sum(new_confusion_matrix[1, ])
print(paste("Specificity on new data:", new_specificity))
new_precision <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[, 2])
new_recall <- new_sensitivity
new_f1_score <- 2 * (new_precision * new_recall) / (new_precision + new_recall)
print(paste("F1-score on new data:", new_f1_score))
```

###Resistance
```{r}
carte_prediction<-(carte_prediction)*100
carte_prediction[is.na(carte_prediction[])] <- 0
carte_prediction<-floor(carte_prediction)

isBecomes <- cbind(c((0:10), (11:20), (21:30), (31:40), (41:50), (51:60), (61:70), (71:80), (81:90), (91:100)), c(rep(100,11), rep(90,10), rep(80,10),rep(70,10), rep(60,10), rep(50,10), rep(40,10),rep(30, 10), rep(20,10) , rep(10,10)))
carte_resistance<- reclassify(carte_prediction, rcl = isBecomes)
carte_resistance<-as.factor(carte_resistance)

r_8bits <- writeRaster(carte_resistance, filename="C:/Users/bambi/OneDrive/Bureau/article/carte_resistance1.tif", datatype="INT1U")

tampons_fusion<-st_read("C:/Users/bambi/OneDrive/Bureau/article", 'tampons_merged')
habitat<-st_read("C:/Users/bambi/OneDrive/Bureau/article/Projet2", 'patches')
raster_resistance <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_resistance1.tif")

route<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin/vecteur_route_pb.shp")

routes_raster <- rasterize(route, raster_resistance, field = 1)
pixels_routes <- raster_resistance * routes_raster

modified_raster_route <- reclassify(pixels_routes, cbind(from = c(9:100), to =105))
mon_modified_masked_raster <- reclassify(modified_raster_route, cbind(from = c(0:80), to = NA))
rc <- cbind(NA, 0)
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)
rc <- cbind(105, 100)
final <- reclassify(r_zero,rc, updatevalue=TRUE)

okok<-final
okok[final  <= 90] <- raster_resistance[final == 0]

mon_masked_raster <- mask(okok, tampons_fusion)
mon_habitat_raster<-mask(okok, habitat)

mon_modified_masked_raster <- reclassify(mon_masked_raster, cbind(from = c(2:100), to = c(101, NA)))
rc <- cbind(c(NA, 1), c(0, 101))
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)

mon_modified_masked_raster_habitat <- reclassify(mon_habitat_raster, cbind(from = c(2:100), to = c(102)))
rc <- cbind(c(NA, NA), c(0, 0))
r_zero_habitat <- reclassify(mon_modified_masked_raster_habitat, rc, updatevalue=TRUE)

final<-r_zero_habitat
final[r_zero_habitat  <= 100] <- okok[r_zero_habitat == 0]

resultat <- r_zero
resultat[r_zero   <= 100] <- final[r_zero  == 0]
```


###Area calculation
```{r}
raster_resistance9<-raster("C:/Users/bambi/OneDrive/Bureau/article/carte_resistance9.tif")
final_presence <- rasterToPolygons(raster_resistance9, dissolve = TRUE)
final_presence1<-st_as_sf(final_presence, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs" )
niveaux <- "102"
shapefile_102 <- subset(final_presence1, carte_resistance9 =="102")
st_area(shapefile_102)
niveaux101 <- "101"
shapefile_101 <- subset(final_presence1, carte_resistance9 == niveaux101)
st_area(shapefile_101)
```

#Random Forest model and adjustment 
```{r}
library(randomForest)
library(caret)
data_complet_merged_ouis_dist$Cllthrx<-as.factor(data_complet_merged_ouis_dist$Cllthrx)
model <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, data = data_complet_merged_ouis_dist, ntree = 500, mtry = sqrt(9))
print(model)

plot(model$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")
set.seed(123)
validationError <- predict(model_sansvariable, data_ouis_test_2.0, type="response")data_complet_merged_ouis_dist$Cllthrx<-as.factor(data_complet_merged_ouis_dist$Cllthrx)
fit <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, data = data_complet_merged_ouis_dist, ntree = 5000, mtry = sqrt(9))
plot(fit$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")

set.seed(123)
fit <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, data = data_complet_merged_ouis_dist, ntree = 3000, mtry = 4)
plot(fit$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")

model <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, data = data_complet_merged_ouis_dist, ntree = 3000, mtry = 4)
print(model)


train_errors <- vector()
test_errors <- vector()
train_sizes <- seq(0.01, 1, by = 0.1)
donnees_shapefile_sans_na$Cllthrx <- as.factor(donnees_shapefile_sans_na$Cllthrx)
data_complet_merged_ouis_dist$Cllthrx <- as.factor(data_complet_merged_ouis_dist$Cllthrx)

for (size in train_sizes) {
    subset_indexes <- sample(1:nrow(donnees_shapefile_sans_na), size = floor(size * nrow(donnees_shapefile_sans_na)))
    subset_train_data <- donnees_shapefile_sans_na[subset_indexes, ]
    
    model <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis + moy_pop_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + moy_vcf_tree_ouis + min_distances, data = data_complet_merged_ouis_dist, ntree = 3000, mtry = 4)
    
    train_pred <- predict(model, data_complet_merged_ouis_dist)
    train_errors <- c(train_errors, sum(train_pred == data_complet_merged_ouis_dist$Cllthrx) / length(train_pred))
    
    test_pred <- predict(model, donnees_shapefile_sans_na)
    test_errors <- c(test_errors, sum(test_pred == donnees_shapefile_sans_na$Cllthrx) / length(test_pred))
}

plot(train_sizes, train_errors, type = "b", pch = 16, col = "blue", ylim = c(0, 1), xlab = "Taille de l'ensemble d'apprentissage", ylab = "Exactitude")
points(train_sizes, test_errors, type = "b", pch = 16, col = "red")
legend("bottomright", c("Ensemble d'apprentissage", "Ensemble de test"), col = c("blue", "red"), pch = 16, lty = 1)
```

###Importance
```{r}
variable_importance <- importance(model)
print(variable_importance)
model$importance[order(model$importance[,1], decreasing=TRUE),]
```

###Probability projection 
```{r}
data_complet_merged_ouis_dist$Cllthrx<-as.factor(data_complet_merged_ouis_dist$Cllthrx)
prediction <- function(model, data_complet_merged_ouis_dist, sefit = FALSE) {
  v <- predict(model, data_complet_merged_ouis_dist, se.fit = sefit, type = "response")
  if (sefit == FALSE) {
    p <- as.vector(v)
  } else {
    p <- as.vector(v$se.fit)
  }
  p
}

points_raster<-crop(x = points_raster, y = precipitation1)
raster_stack <- stack(ndvi2, temperature1, precipitation1, vapeur1, ra1, pop1, land, vege2, vcf_non_vege, vcf_tree, points_raster)
names(raster_stack )<- c("moy_ndvi_ouis", "moy_temp_ouis", "moy_prec_ouis", "moy_vap_ouis", "moy_ra_ouis", "moy_pop_ouis", "moy_land_ouis", "moy_vege2_ouis", "moy_vcf_non_vege_ouis", "moy_vcf_tree_ouis", "min_distances")

carte_prediction_RF <- 1-raster::predict(model=model, object=raster_stack, type="prob")

pred <- raster::predict(model, object = raster_stack, type = "prob", index=2)
ecarts_types <- sqrt(pred * (1 - pred))
```
###Calculation of species presence probabilities based on forest fragments
```{r}
carte_prediction_RF<-raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1.tif")

extraction<-raster::extract(x = carte_prediction_RF, y = aireetude, df = TRUE)
extraction <- na.omit(extraction)
extraction_ouis<-tapply(extraction $carte_prediction_RF, extraction$ID, mean)

values <- c(0.74722221, 0.28908974, 0.50192592, 0.10288889, 0.17636394, 0.30682157, 0.40066082, 0.99183333, 0.34045513, 0.23441333, 0.36757895, 0.10011765, 0.06796078)
data_extraction <- data.frame(ID = 1:13, Value = values)
# Assigner les noms de localités à chaque ID
localites <- c("Imbaú III", " Aldeia I", "Aldeia II", "Cabo Frio", "Dourada", "Gaviões", "Imbaú I", "Imbaú II", "Poço das Antas", "Pirineus", "Rio Vermelho", "União I", "União II")
data_extraction$Localites <- factor(localites, levels = localites)

data_extraction
Localities<-data_extraction$Localites
Mean<-data_extraction$Value
ggplot(data_extraction, aes(x = Localities, y = Mean)) +
  geom_bar(stat = "identity", fill = "darkgreen")+
  theme(panel.background = element_rect(fill = "white"))

```
###KAPPA
```{r}
prediction<-predict(model, donnees_shapefile_sans_na)
prediction <- ifelse(prediction > 0.5,1,0 )

confusion_matrix <- table(Predicted = predictions, Actual = donnees_shapefile_sans_na$Cllthrx)

a <- confusion_matrix[1,1] 
b <- confusion_matrix[1,2] 
c <- confusion_matrix[2,1] 
d <- confusion_matrix[2,2] 


n <- sum(confusion_matrix) 
po <- (a + d) / n 
pe <- ((a + b) * (a + c) + (c + d) * (b + d)) / (n^2) 
kappa <- (po - pe) / (1 - pe)

print(kappa)
```
###Matthews Index (MCC)
```{r}
MCC <- (a * d - b * c) / sqrt((a + b) * (a + c) * (d + b) * (d + c))
```

###Accuracy
```{r}
donnees_shapefile_sans_na$moy_vege2_ouis<-as.factor(donnees_shapefile_sans_na$moy_vege2_ouis)
predictions_acc <- predict(model.final.dist, donnees_shapefile_sans_na, type = "response")
predictions_acc <- ifelse(predictions_acc > 0.5,1,0 ) 
accuracy <- sum(predictions_acc == donnees_shapefile_sans_na$Cllthrx) / length(donnees_shapefile_sans_na$Cllthrx)
print(paste("Accuracy:", accuracy)) 
```

###Cross validation training data
```{r}
library(randomForest)
library(pROC)
data_complet_merged_ouis_dist$Cllthrx<-as.factor(data_complet_merged_ouis_dist$Cllthrx)
K <- 5
nk <- floor(nrow(data_complet_merged_ouis_dist) / K)
pred <- NULL
Yobs <- NULL
rk <- list()
set.seed(48)
ii <- sample(nrow(data_complet_merged_ouis_dist), nrow(data_complet_merged_ouis_dist)) # Permutation aléatoire des indices

for (k in 1:K) {
  testi <- ii[((k-1) * nk + 1):(k * nk)]
  appri <- setdiff(1:nrow(data_complet_merged_ouis_dist), testi)
  model <- randomForest(Cllthrx ~ moy_vap_ouis + moy_ra_ouis + moy_temp_ouis + moy_ndvi_ouis+ moy_pop_ouis+ moy_vege2_ouis+ moy_vcf_non_vege_ouis +moy_vcf_tree_ouis +min_distances, data = data_complet_merged_ouis_dist, ntree = 500, mtry = sqrt(9))
  pr <- predict(model, data_complet_merged_ouis_dist[testi,], type = "prob")[, 2]
  pred <- c(pred, pr)
  Yobs <- c(Yobs, data_complet_merged_ouis_dist$Cllthrx[testi])

  rk[[k]] <- roc(data_complet_merged_ouis_dist$Cllthrx[testi], pr)
}

add <- FALSE
for (k in 1:K) {
  if (k > 1) add <- TRUE
  plot.roc(rk[[k]], col = "gray", add = add, lwd = 1)
}
r <- roc(Yobs, pred)
plot(r, main = paste("AUROC =", round(r$auc * 100) / 100), add = TRUE)
grid()
auc(r)
auc(rk[[K]])
```

###Cross validation sample test
```{r}
K <- 5
nk <- floor(nrow(donnees_shapefile_sans_na) / K)
pred <- NULL
Yobs <- NULL
rk <- list()
set.seed(48)
ii <- sample(nrow(donnees_shapefile_sans_na), nrow(donnees_shapefile_sans_na)) # Permutation aléatoire des indices
for (k in 1:K) {
    testi <- ii[((k-1) * nk + 1):(k * nk)]
    appri <- setdiff(1:nrow(donnees_shapefile_sans_na), testi)
    model <-randomForest(Cllthrx ~ moy_vap_ouis + moy_vege2_ouis + moy_vcf_non_vege_ouis + moy_vcf_tree_ouis,data = data_complet_merged_ouis_dist, ntree = 3000, mtry = 4)
    pr <- predict(model, donnees_shapefile_sans_na[testi,], type = "prob")[, 2]
    pred <- c(pred, pr)
    Yobs <- c(Yobs,donnees_shapefile_sans_na$Cllthrx[testi])
    rk[[k]] <- roc(donnees_shapefile_sans_na$Cllthrx[testi], pr)
}

add <- FALSE
for (k in 1:K) {
    if (k > 1) add <- TRUE
    plot.roc(rk[[k]], col = "gray", add = add, lwd = 1)
}

r <- roc(Yobs, pred)
plot(r, main = paste("AUROC =", round(r$auc * 100) / 100), add = TRUE)
grid()
auc(r)
auc(rk[[K]])
```


###Accuracy
```{r}
predictions_acc <- predict(model, donnees_shapefile_sans_na, type = "response")
predictions_acc <- ifelse(predictions_acc > 0.5,1,0 )
accuracy <- sum(predictions_acc== donnees_shapefile_sans_na$Cllthrx) / length(donnees_shapefile_sans_na$Cllthrx)
print(paste("Accuracy:", accuracy))

```

###Sensitivity, specificity, F1-score train data
```{r}
train_predictions <- predict(model, data_complet_merged_ouis_dist)
train_confusion_matrix <- table(data_complet_merged_ouis_dist$Cllthrx, train_predictions)
train_sensitivity <- train_confusion_matrix[2, 2] / sum(train_confusion_matrix[2, ])
print(paste("Sensitivity (Recall) on training data:", train_sensitivity))
train_specificity <- train_confusion_matrix[1, 1] / sum(train_confusion_matrix[1, ])
print(paste("Specificity on training data:", train_specificity))
train_precision <- train_confusion_matrix[2, 2] / sum(train_confusion_matrix[, 2])
train_recall <- train_sensitivity
train_f1_score <- 2 * (train_precision * train_recall) / (train_precision + train_recall)
print(paste("F1-score on training data:", train_f1_score))
```


###Sensitivity, specificity, F1-score validation data
```{r}
new_predictions <- predict(model, donnees_shapefile_sans_na) 
new_predictions <- ifelse(new_predictions > 0.5,1,0 ) 
new_confusion_matrix <- table(donnees_shapefile_sans_na$Cllthrx, new_predictions)
new_sensitivity <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[2, ])
print(paste("Sensitivity (Recall) on new data:", new_sensitivity))
new_specificity <- new_confusion_matrix[1, 1] / sum(new_confusion_matrix[1, ])
print(paste("Specificity on new data:", new_specificity))
new_precision <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[, 2])
new_recall <- new_sensitivity
new_f1_score <- 2 * (new_precision * new_recall) / (new_precision + new_recall)
print(paste("F1-score on new data:", new_f1_score))
```

###Correlation between models  
```{r}
pred_RF <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1.tif")
pred_GLM<- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction1i.tif")
pred_RF<-crop(x=pred_RF, y=precipitation1)
pred_GLM<-crop(x=pred_GLM, y=precipitation1)

val_pred_RF <- getValues(pred_RF)
val_pred_GLM <- getValues(pred_GLM)

corr <- cor(val_pred_RF, val_pred_GLM)
cor.test(val_pred_RF, val_pred_GLM)


library(vegan)
mat_pred_RF<- as.matrix(pred_RF)
mat_pred_GLM <- as.matrix(pred_GLM)
matrice1 <- dist(mat_pred_RF)  
matrice2 <- dist(mat_pred_GLM)  

resultat_mantel <- mantel(matrice1, matrice2, method = "spearman", permutations = 999)
```

###Resistance matrix
```{r}
raster_prediction <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1.tif")
carte_prediction<-(raster_prediction)*100
carte_prediction[is.na(carte_prediction[])] <- 0
carte_prediction<-floor(carte_prediction)

isBecomes <- cbind(c((0:10), (11:20), (21:30), (31:40), (41:50), (51:60), (61:70), (71:80), (81:90), (91:100)), c(rep(100,11), rep(90,10), rep(80,10),rep(70,10), rep(60,10), rep(50,10), rep(40,10),rep(30, 10), rep(20,10) , rep(10,10)))
carte_resistance<- reclassify(carte_prediction, rcl = isBecomes)
carte_resistance<-as.factor(carte_resistance)

r_8bits <- writeRaster(carte_resistance,filename="C:/Users/bambi/OneDrive/Bureau/article/carte_resistance_RF4.tif", datatype="INT1U")

tampons_fusion<-st_read("C:/Users/bambi/OneDrive/Bureau/article", 'tampons_merged')
habitat<-st_read("C:/Users/bambi/OneDrive/Bureau/article/ProjetRF4_50", 'patches')
raster_resistance <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_resistance_RF4.tif")

route<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin/vecteur_route_pb.shp")

routes_raster <- rasterize(route, raster_resistance, field = 1)
pixels_routes <- raster_resistance * routes_raster
plot(pixels_routes)


modified_raster_route <- reclassify(pixels_routes, cbind(from = c(9:100), to =105))
mon_modified_masked_raster <- reclassify(modified_raster_route, cbind(from = c(0:80), to = NA))
rc <- cbind(NA, 0)
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)
rc <- cbind(105, 100)
final <- reclassify(r_zero,rc, updatevalue=TRUE)


okok<-final
okok[final  <= 90] <- raster_resistance[final == 0]

mon_masked_raster <- mask(okok, tampons_fusion)
mon_habitat_raster<-mask(okok, habitat)

mon_modified_masked_raster <- reclassify(mon_masked_raster, cbind(from = c(2:100), to = c(101, NA)))
rc <- cbind(c(NA, 1), c(0, 101))
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)

mon_modified_masked_raster_habitat <- reclassify(mon_habitat_raster, cbind(from = c(2:100), to = c(102)))
rc <- cbind(c(NA, NA), c(0, 0))
r_zero_habitat <- reclassify(mon_modified_masked_raster_habitat, rc, updatevalue=TRUE)


final<-r_zero_habitat
final[r_zero_habitat  <= 100] <- okok[r_zero_habitat == 0]
resultat <- r_zero
resultat[r_zero   <= 100] <- final[r_zero  == 0]
```

###Area calculation
```{r}
carte_prediction_RF10_50_2<-raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF10_50_2.tif")
final_presence <- rasterToPolygons(carte_prediction_RF10_50_2, dissolve = TRUE)
final_presence1<-st_as_sf(final_presence, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs" )
niveaux <- "102"
shapefile_102 <- subset(final_presence1, carte_prediction_RF10_50_2=="102")
st_area(shapefile_102)
niveaux101 <- "101"
shapefile_101 <- subset(final_presence1, carte_prediction_RF10_50_2 == niveaux101)
st_area(shapefile_101)
```

### Id vegetation
```{r}
vegetation$leg_carga<-as.factor(vegetation$leg_carga)

levels(vegetation$leg_carga)[data_complet_merged_ouis_dist$moy_vege2_ouis]
View(vegetation)
plot(data_complet_merged_ouis_dist$moy_vege2_ouis)
cbind(levels(vegetation$leg_carga),seq(1,length(levels(vegetation$leg_carga))))
cbind(levels(vegetation$leg_carga),seq(1,length(levels(vegetation$leg_sup))))
```
#Sapajus
```{r}
data_article<-st_read("C:/Users/bambi/OneDrive/Bureau/article", 'master_pa')
data_article<-st_transform(data_article, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
data_article <- st_intersection(data_article, h)

data1<-data_article
sap<-data1$Sapajus
presence_sap <- data1[sap == "1",]
absence_sap <- data1[sap == "0",]
tampons_presence_sap <- st_buffer(presence_sap, dist = 912.61)
tampons_presence_sap_pourmap <- st_buffer(presence_sap, dist = 564.19)
absence_sap<-st_as_sf(absence_sap, coords=data1$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
presence_sap<-st_as_sf(presence_sap, coords=data1$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

points_non_contenus <- data1[6]
for(i in 1:nrow(tampons_presence_sap)) {
  points_non_contenus <- points_non_contenus[!st_intersects(points_non_contenus, tampons_presence_sap[i,], sparse = FALSE),]
}

data2 <- data.frame()
for(file in c("points_non_contenus", "presence_sap")){
  temp_file <- get(file)
  data2 <- rbind(data2, temp_file[c("geometry", "Sapajus")])
}
data2<- na.omit(data2)


sap<-data2$Sapajus
Sapajus <- data2[sap == "1",]

data2_0<- st_as_sfc(points_non_contenus)
empty_geoms <- st_is_empty(data2_0)
data2_0 <- st_drop_geometry(data2_0[!empty_geoms])
data_sp_0 <-as(data2_0, "Spatial")
cells <- cellFromXY(temperature1, data_sp_0)
dups <- duplicated(cells)
data_pixel_0 <- data_sp_0[!dups, ]
cat(nrow(data_sp_0) - nrow(data_pixel_0), "records are removed")


data_pixel_sf_obj <- as(data_pixel_0, "sf")
data_pixel_0<-st_as_sf(data_pixel_sf_obj, coords=data_pixel_sf_obj$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


data2_1<- st_as_sfc(Sapajus)
empty_geoms <- st_is_empty(data2_1) 
data2_1 <- st_drop_geometry(data2_1[!empty_geoms])
data_sp_1 <-as(data2_1, "Spatial")


data_pixel_sf_obj1 <- as(data_sp_1, "sf")
data_sp_1<-st_as_sf(data_pixel_sf_obj1, coords=data_pixel_sf_obj1$geometry, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


absent<-rep(0,872)
data_pixel_0<-data.frame(data_pixel_0, absent)
names(data_pixel_0)<-c("geometry", "Sapajus")
data_pixel_0<-st_as_sf(data_pixel_0, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

pre<-rep(0,78)
data_pixel_1<-data.frame(data_sp_1, pre)
names(data_pixel_1)<-c("geometry", "Sapajus")
data_pixel_1<-st_as_sf(data_pixel_1, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


data2<- data.frame()
for(file in c("data_pixel_0", "Sapajus")){
  temp_file <- get(file)
  data2 <- rbind(data2, temp_file[c("geometry", "Sapajus")])
}
```

###Absences
```{r}
absence_sap <- data2[sap == "0",]
absence_sap<-st_as_sf(absence_sap,"+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs" )


sap<-data2$Sapajus
Sapajus <- data2[sap == "1",]
absence_sap<-data2[sap == "0",]

nc_nonempty_absence_sap <- absence_sap[!st_is_empty(absence_sap), ]
absence_sap1<-st_as_sfc(absence_sap)

abd<-rep(0,78)
absence_sap1<-st_cast(st_sample(nc_nonempty_absence_sap, size=78), "POINT")
absence_sap1<-data.frame(absence_sap1, abd)
names(absence_sap1)<-c("geometry","Sapajus")
absence_sap1<-st_as_sf(absence_sap1, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")


tampons_absence_sap <- st_buffer(absence_sap1, dist = 912.61)
tampons_absence_sap<-data.frame(tampons_absence_sap)
names(tampons_absence_sap)<-c("geometry", "Sapajus")
tampons_absence_sap<-st_as_sf(tampons_absence_sap, crs="+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

tampons_absence_sap <- as_Spatial(tampons_absence_sap)

temp_abs_sap <-raster::extract(x = Nom_raster_multlicouche, y = tampons_absence_sap, df = TRUE)

head(temp_abs_sap)
moy_prec_sap<-tapply(temp_abs_sap $moy_prec_ouis, temp_abs_sap $ID, mean)
moy_temp_sap<-tapply(temp_abs_sap $moy_temp_ouis, temp_abs_sap $ID, mean)
moy_vap_sap<-tapply(temp_abs_sap $moy_vap_ouis, temp_abs_sap $ID, mean)
moy_ra_sap<-tapply(temp_abs_sap $moy_ra_ouis, temp_abs_sap $ID, mean)
moy_ndvi_sap<-tapply(temp_abs_sap $ndvi1, temp_abs_sap $ID, mean)
moy_pop_sap<-tapply(temp_abs_sap $moy_pop_ouis, temp_abs_sap $ID, mean)
moy_land_sap<-extract(land, absence_sap1)
moy_vege2_sap<-extract(vege2, absence_sap1)
moy_vcf_non_vege_sap<-tapply(temp_abs_sap $moy_vcf_non_vege_ouis, temp_abs_sap $ID, mean)
moy_vcf_tree_sap<-tapply(temp_abs_sap $moy_vcf_tree_ouis, temp_abs_sap $ID, mean)
moy_prec_sap<-tapply(temp_abs_sap $moy_prec_ouis, temp_abs_sap $ID, mean)

absence_sap2_sfc<- st_as_sfc(absence_sap1)
shapefile_sub<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin/vecteur_route_pb.shp")
shapefile_sub_sfc<- st_as_sfc(shapefile_sub[6])
distance<-st_distance(absence_sap2_sfc, shapefile_sub_sfc)
min_distances<- apply(distance, 1, min)

data_merged_sap_abs <- data.frame(moy_temp_sap,moy_vap_sap, moy_ra_sap,moy_ndvi_sap,moy_pop_sap, moy_vcf_non_vege_sap, moy_vcf_tree_sap,moy_land_sap,moy_vege2_sap, moy_prec_sap,min_distances, absence_sap1)

```

###Presences
```{r}
sap<-data2$Sapajus
presence_sap <- data2[sap == "1",]
presence_sap<-st_as_sf(presence_sap,"+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs" )

tampons_presence_sap <- st_buffer(presence_sap, dist = 912.61)
tampons_presence_sap<-st_write(tampons_presence_sap,"C:/Users/bambi/OneDrive/Bureau/article/tampons_sap.shp")

temp_pres_sap <- raster::extract(x = Nom_raster_multlicouche, 
                                y = tampons_presence_sap  , 
                                df = TRUE)

head(temp_pres_ouis)
moy_temp_sap<-tapply(temp_pres_sap$moy_temp_oui, temp_pres_sap$ID, mean)
moy_vap_sap<-tapply(temp_pres_sap$moy_vap_oui, temp_pres_sap$ID, mean)
moy_ra_sap<-tapply(temp_pres_sap$moy_ra_oui, temp_pres_sap$ID, mean)
moy_ndvi_sap<-tapply(temp_pres_sap$ndvi1, temp_pres_sap$ID, mean)
moy_pop_sap<-tapply(temp_pres_sap$moy_pop_ouis, temp_pres_sap$ID, mean)
moy_land_sap<-extract(land, presence_sap)
moy_vege2_sap<-extract(vege2, presence_sap)
moy_vcf_non_vege_sap<-tapply(temp_pres_sap$moy_vcf_non_vege_ouis, temp_pres_sap$ID, mean)
moy_vcf_tree_sap<-tapply(temp_pres_sap$moy_vcf_tree_ouis, temp_pres_sap$ID, mean)
moy_prec_sap<-tapply(temp_pres_sap$moy_prec_ouis, temp_pres_sap$ID, mean)

presence_sap2_sfc<- st_as_sfc(presence_sap)
shapefile_sub<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin/vecteur_route_pb.shp")
shapefile_sub_sfc<- st_as_sfc(shapefile_sub[6])
distance<-st_distance(presence_sap2_sfc, shapefile_sub_sfc)
min_distances<- apply(distance, 1, min)


data_merged_sap_pres <- data.frame(moy_temp_sap,moy_vap_sap, moy_ra_sap,moy_pop_sap, moy_ndvi_sap, moy_vcf_non_vege_sap, moy_vcf_tree_sap,moy_land_sap,moy_vege2_sap, moy_prec_sap, min_distances, presence_sap)
```


###Table
```{r}
data_complet_merged_sap <- data.frame()
for(file in c("data_merged_sap_pres", "data_merged_sap_abs")){
  temp_file <- get(file)
  data_complet_merged_sap <- rbind(data_complet_merged_sap, temp_file[c("geometry","moy_vege2_sap", "moy_pop_sap","moy_vcf_non_vege_sap", "moy_temp_sap", "moy_vap_sap", "moy_ra_sap",  "moy_ndvi_sap", "moy_vcf_tree_sap" , "moy_land_sap", "moy_prec_sap", "min_distances", "Sapajus")])
}
```

###ACP
```{r}
library(ade4)
ACP_sap<-dudi.pca(data_complet_merged_sap[,c(15:21, 24)], center =TRUE, scale = TRUE, scannf =FALSE)
summary(ACP_sap)
s.corcirle(ACP_sap$co)
dudi.mix

library(FactoMineR)
PCALarge<-PCA(data_complet_merged_sap [, c(15:21, 24)], graph =TRUE)
corrplot(cor(data_complet_merged_sap[, c(15:21, 24)], use="complete.obs"))

library(corrplot)
corr_matrix<-cor(data_complet_merged_sap [, c(15:21, 24)], use = "complete.obs")
corrplot(cor_matrix, method ="number", tl.cex = 0.7)

```


#GLM Sapajus 
```{r}
data_complet_merged_sap$moy_vege2_sap<-as.factor(data_complet_merged_sap$moy_vege2_sap)
data_complet_merged_sap$moy_land_sap<-as.factor(data_complet_merged_sap$moy_land_sap)
data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)
data_complet_merged_sap$moy_temp_sap<-as.numeric(data_complet_merged_sap$moy_temp_sap)
data_complet_merged_sap$moy_vap_sap<-as.numeric(data_complet_merged_sap$moy_vap_sap)
data_complet_merged_sap$moy_ra_sap<-as.numeric(data_complet_merged_sap$moy_ra_sap)
data_complet_merged_sap$moy_ndvi_sap<-as.numeric(data_complet_merged_sap$moy_ndvi_sap)
data_complet_merged_sap$moy_pop_sap<-as.numeric(data_complet_merged_sap$moy_pop_sap)
data_complet_merged_sap$moy_vcf_non_vege_sap<-as.numeric(data_complet_merged_sap$moy_vcf_non_vege_sap)
data_complet_merged_sap$moy_vcf_tree_sap<-as.numeric(data_complet_merged_sap$moy_vcf_tree_sap)
data_complet_merged_sap$min_distances<-as.numeric(data_complet_merged_sap$min_distances)


null.m <- glm(Sapajus ~ 1, family = binomial(link = "logit"), data =
                data_complet_merged_sap)


Nom_modele_step <- stepAIC(null.m, Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap+ min_distances, family=binomial(link = "logit"), data = data_complet_merged_sap)


model.final.sap <- glm(Sapajus ~ min_distances, family=binomial(link = "logit"), data=data_complet_merged_sap)
```


###GLM model predictions and probability projection
```{r}
prediction <- function(model.final.sap, data_complet_merged_sap, sefit=FALSE) {
  v <-stats::predict(model.final.sap, data_complet_merged_sap, se.fit=T, type="response")
  if(sefit==FALSE){p=as.vector(v$fit)} else {p=as.vector(v$se.fit)}
  p
}

raster_stack_sap <- stack(ndvi2, temperature1, precipitation1, vapeur1, ra1, pop1, land, vege2, vcf_non_vege, vcf_tree, points_raster)
names(raster_stack_sap )<- c("moy_ndvi_sap", "moy_temp_sap", "moy_prec_sap", "moy_vap_sap", "moy_ra_sap", "moy_pop_sap", "moy_land_sap", "moy_vege2_sap", "moy_vcf_non_vege_sap", "moy_vcf_tree_sap", "min_distances")

carte_prediction<-predict(raster_stack_sap,
                          model=model.final.sap, fun=prediction,
                          index=11)

carte_prediction[is.na(carte_prediction[])] <- 0
plot(carte_prediction)


carte_prediction_sd<-predict(raster_stack_sap,
model=model.final.sap, fun=prediction, sefit=TRUE,
index=1)
carte_prediction_sd[is.na(carte_prediction_sd[])] <- 0
writeRaster(carte_prediction_sd, "C:/Users/bambi/OneDrive/Bureau/article/GLM carte pred/carte_prediction_GLM_sap_sd.tif", overwrite=TRUE)
```

###Model coefficient
```{r}
coef<-exp(coefficients(model.final.sap))
exp(confint(model.final.sap))
```

###PseudoR²
```{r}
model_summary <- summary(model.final.sap)
r_squared <- model_summary$deviance / model_summary$null.deviance
print(r_squared)
```

###ANOVA 
```{r}
anova(model.final.sap,test="Chisq")
par_var_expl_lmodel.final <- anova(model.final.sap)[-
1,2]*100 / anova(model.final.sap)[1,4]
print( cbind(rownames(anova(model.final.sap))[-
1],round(par_var_expl_lmodel.final  , digit = 1)))
```


###Graphics and types of vegetation
```{r}
library(GGally)
library(ggplot2)
library(jtools)

plot <- ggcoef_model(model.final.sap, exponentiate =TRUE, variable_labels = c(min_distances = "Distance to roads"))

plot +
    theme(
        axis.title.x = element_text(size = 14),  
        axis.title.y = element_text(size = 14),
        axis.text.x = element_text(size = 14), 
        axis.text.y = element_text(size = 14)
    )


data_complet_merged_sap$moy_vap_sap1<-data_complet_merged_sap$moy_vap_sap/100
da<-data.frame(data_complet_merged_sap$moy_vap_sap1, data_complet_merged_sap$Sapajus)
data_complet_merged_sap$moy_vap_sap1<-as.numeric(data_complet_merged_sap$moy_vap_sap1)
data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)



courbe_min_distances_sap <- ggeffect(model.final.sap, terms = c("min_distances"))
plot(courbe_min_distances_sap )
a <- ggplot() +
    geom_line(data = courbe_min_distances_sap, aes(x = x, y = predicted), color = "black", linewidth = 1) +  
    labs(x = "Distance to roads", y = "Predicted probability of presence") +
    theme(
        panel.background = element_rect(fill = "white", color = "white"),
        plot.background = element_rect(fill = "white", color = "white"),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14)
    )



vegetation<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/vege", 'vegeRJ')
vegetation<-st_transform(vegetation, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")#en utm23CRS
vegetation$leg_carga<-as.factor(vegetation$leg_carga)
levels(vegetation$leg_carga)[data_complet_merged_sap$moy_vege2_sap]
cbind(levels(vegetation$leg_carga),seq(1,length(levels(vegetation$leg_carga))))

vegetation$leg_sup<-as.factor(vegetation$leg_sup)
raster_vegetation_anth <- rasterize(vegetation, precipitation1, "leg_sup")

sap<-data_complet_merged_sap$Sapajus
Sapajus <- data_complet_merged_sap[sap == "1",]
absence_sap<-data_complet_merged_sap[sap == "0",]
absence_sap1<-st_as_sfc(absence_sap)

absence_sap1 <-as(absence_sap1, "Spatial")
absence_anth<-extract(raster_vegetation_anth, absence_sap1)
ab<-rep(0,78)
absence_anth<-data.frame(absence_anth, ab)


sap<-data_complet_merged_sap$Sapajus
presence_sap<-data_complet_merged_sap[sap == "1",]
presence_sap<- st_as_sfc(presence_sap)
presence_sap <-as(presence_sap, "Spatial")
absence_anth<-extract(raster_vegetation_anth, presence_sap)
ab<-rep(1,78)
presence_anth<-data.frame(absence_anth, ab)

data_anth <- data.frame()
for(file in c("presence_anth", "absence_anth")){
  temp_file <- get(file)
  data_anth <- rbind(data_anth, temp_file[c("absence_anth", "ab")])
}
presence_anth$absence_anth<-as.factor(presence_anth$absence_anth)
presence_anth$ab<-as.factor(presence_anth$ab)

plot_anth <- ggplot(presence_anth, aes(x = presence, fill = presence)) + 
  geom_bar(position = "dodge") + 
  scale_fill_manual(values = c("lightgreen", "darkgreen"),
                    labels = c("Dominant anthropic area", "Dominant natural vegetation")) +
  labs(x = "Type of vegetation", y = "Frequency") +
  ggtitle("Distribution of capuchins according to vegetation type: anthropic or natural") +
  theme_classic() +
  guides(fill = guide_legend(title = NULL)) +  
  theme(legend.position = "right")


presence_anth14 <- subset(presence_anth, absence_anth %in% c("1", "4"))
contingence14 <- table(presence_anth14$ab, presence_anth14$absence_anth)
j<-chisq.test(contingence14)
j$stdres



sap<-data_complet_merged_sap$Sapajus
presence_sap2<-data_complet_merged_sap[sap == "1",]
presence_sap2<- st_as_sfc(presence_sap2)
presence_sap2 <-as(presence_sap2, "Spatial")

vegetation$leg_sup<-as.factor(vegetation$leg_sup)
raster_vegetation_anth <- rasterize(vegetation, precipitation1, "leg_sup")

presence<-extract(raster_vegetation_anth, presence_sap2)
ab<-rep(1,78)
presence_anth<-data.frame(presence, ab)


presence_anth$presence<-as.factor(presence_anth$presence)
presence_anth$ab<-as.factor(presence_anth$ab)

plot_anth <- ggplot(presence_anth, aes(x = presence, fill = presence)) + 
  geom_bar(position = "dodge") + 
  scale_fill_manual(values = c("darkgreen", "darkgreen")) +
  labs(x = "Type of vegetation", y = "Frequency") +
  ggtitle("Distribution of capuchins according to vegetation type: anthropic or natural") +
  theme_classic() +
  guides(fill=guide_legend(title=NULL)) # Remove legend title
plot_anth <- plot_anth + theme(legend.position = "bottom") 
plot_anth
```



###Cross validation
```{r}
K = 5
nk = floor(nrow(data_complet_merged_sap)/K)
pred = NULL
Yobs = NULL
rk=list()
set.seed(48)
ii = sample(nrow(data_complet_merged_sap),nrow(data_complet_merged_sap)) # permutation aléatoire des indices
for (k in 1:K){
  testi = ii[((k-1)*nk+1):(k*nk)]
  appri = setdiff(1:nrow(data_complet_merged_sap),testi)
  blm<-glm(Sapajus ~ min_distances + moy_ndvi_sap, family = binomial(link = "logit"), data = data_complet_merged_sap)
           
  pr = predict(blm,data_complet_merged_sap[testi,],type="response")
  pred = c(pred,pr)
  Yobs = c(Yobs,data_complet_merged_sap$Sapajus[testi])
  rk[[k]]=roc(data_complet_merged_sap$Sapajus[testi],pr)
}
add = FALSE
for (k in 1:K){
  if (k>1) add=TRUE
  plot.roc(rk[[k]],col="gray",add=add,lwd=1)
}
r=roc(Yobs,pred)
plot(r,main=paste("AUROC =",round(r$auc*100)/100),add=TRUE)
grid()
auc(r)
auc(rk[[k]])



pred_test <- predict(model.final.sap, data_complet_merged_sap, type="response")
test_roc = roc(data_complet_merged_sap$Sapajus ~ pred_test, plot = TRUE, print.auc = TRUE)
plot(test_roc,auc.polygon=TRUE)
auc(test_roc)
```

###KAPPA
```{r}
prediction_sap<-predict(model.final.sap, data_complet_merged_sap, type='response')
prediction_sap <- ifelse(prediction_sap > 0.5,1,0 )

confusion_matrix <- table(Predicted = prediction_sap, Actual = data_complet_merged_sap$Sapajus)

a <- confusion_matrix[1,1] 
b <- confusion_matrix[1,2] 
c <- confusion_matrix[2,1] 
d <- confusion_matrix[2,2] 

n <- sum(confusion_matrix) 
po <- (a + d) / n 
pe <- ((a + b) * (a + c) + (c + d) * (b + d)) / (n^2) 
kappa <- (po - pe) / (1 - pe)
print(kappa)
```


###Matthews index (MCC)
```{r}
MCC <- (a * d - b * c) / sqrt((a + b) * (a + c) * (d + b) * (d + c))
```


###Sensitivity, specificity, F1-score validation data
```{r}
predictions_acc <- predict(model.final.sap, data_complet_merged_sap, type = "response")
predictions_acc <- ifelse(predictions_acc > 0.5,1,0 )
new_predictions <- predictions_acc
new_confusion_matrix <- table(data_complet_merged_sap$Sapajus, new_predictions)

new_sensitivity <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[2, ])
print(paste("Sensitivity (Recall) on new data:", new_sensitivity))

new_specificity <- new_confusion_matrix[1, 1] / sum(new_confusion_matrix[1, ])
print(paste("Specificity on new data:", new_specificity))

new_precision <- new_confusion_matrix[2, 2] / sum(new_confusion_matrix[, 2])
new_recall <- new_sensitivity
new_f1_score <- 2 * (new_precision * new_recall) / (new_precision + new_recall)
print(paste("F1-score on new data:", new_f1_score))
```

###Accuracy
```{r}
predictions_acc <- predict(model.final.sap, data_complet_merged_sap, type = "response")
predictions_acc <- ifelse(predictions_acc > 0.5,1,0 )
accuracy <- sum(predictions_acc== data_complet_merged_sap$Sapajus) / length(data_complet_merged_sap$Sapajus)
print(paste("Accuracy:", accuracy))
```

#Random Forest model
```{r}
library(randomForest)
library(caret)
data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)
model <- randomForest(Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap + min_distances, data = data_complet_merged_sap, ntree = 500, mtry = sqrt(9))
print(model)

plot(model$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")
set.seed(123)
data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)
fit <- randomForest(Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap +min_distances, data = data_complet_merged_sap, ntree = 8100, mtry = sqrt(9))
plot(fit$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")

 
set.seed(123)
fit <- randomForest(Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap +min_distances, data = data_complet_merged_sap, ntree = 8100, mtry = 5)
plot(fit$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")


model_sap <- randomForest(Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap +min_distances, data = data_complet_merged_sap, ntree = 8100, mtry = 5)
print(model_sap)
```

###KAPPA
```{r}
prediction_sap<-predict(model_sap, data_complet_merged_sap)

confusion_matrix <- table(Predicted = prediction_sap, Actual = data_complet_merged_sap$Sapajus)
a <- confusion_matrix[1,1] 
b <- confusion_matrix[1,2] 
c <- confusion_matrix[2,1]
d <- confusion_matrix[2,2] 

n <- sum(confusion_matrix) 
po <- (a + d) / n 
pe <- ((a + b) * (a + c) + (c + d) * (b + d)) / (n^2) 
kappa <- (po - pe) / (1 - pe)

print(kappa)
```

###Matthews Indice (MCC)
```{r}
MCC <- (a * d - b * c) / sqrt((a + b) * (a + c) * (d + b) * (d + c))
```

###Importance
```{r}
model_sap$importance[order(model_sap$importance[,1], decreasing=TRUE),]
```

### ANOVA
```{r}
anova(model_sap,test="Chisq")
par_var_expl_lmodel.final <- anova(model_sap)[-
1,2]*100 / anova(model_sap)[1,4]
print( cbind(rownames(anova(model_sap))[-
1],round(par_var_expl_lmodel.final  , digit = 1)))
```

###AUC cross validation
```{r}
data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)
K <- 5
nk <- floor(nrow(data_complet_merged_sap) / K)

pred <- NULL
Yobs <- NULL

rk <- list()

set.seed(48)
ii <- sample(nrow(data_complet_merged_sap), nrow(data_complet_merged_sap))

for (k in 1:K) {
  testi <- ii[((k-1) * nk + 1):(k * nk)]
  appri <- setdiff(1:nrow(data_complet_merged_sap), testi)
  
  model <- randomForest(Sapajus ~ moy_vap_sap + moy_ra_sap + moy_temp_sap + moy_ndvi_sap+ moy_pop_sap+ moy_vege2_sap+ moy_vcf_non_vege_sap +moy_vcf_tree_sap +min_distances, data = data_complet_merged_sap, ntree = 8100, mtry = 5)
  
  pr <- predict(model, data_complet_merged_sap[testi,], type = "prob")[, 2]
  
  pred <- c(pred, pr)
  Yobs <- c(Yobs, data_complet_merged_sap$Sapajus[testi])
  
  rk[[k]] <- roc(data_complet_merged_sap$Sapajus[testi], pr)
}

add <- FALSE
for (k in 1:K) {
  if (k > 1) add <- TRUE
  plot.roc(rk[[k]], col = "gray", add = add, lwd = 1)
}

r <- roc(Yobs, pred)

plot(r, main = paste("AUROC =", round(r$auc * 100) / 100), add = TRUE)
grid()

auc(r)
```

###Accuracy
```{r}
predictions_acc <- predict(model_sap, data_complet_merged_sap, type = "response")

accuracy <- sum(predictions_acc== data_complet_merged_sap$Sapajus) / length(data_complet_merged_sap$Sapajus)
print(paste("Accuracy:", accuracy))
```


###Sensibilité, spécificité, F1-score train data
```{r}
train_predictions <- predict(model_sap, data_complet_merged_sap)
train_confusion_matrix <- table(data_complet_merged_sap$Sapajus, train_predictions)
train_sensitivity <- train_confusion_matrix[2, 2] / sum(train_confusion_matrix[2, ])
print(paste("Sensitivity (Recall) on training data:", train_sensitivity))

train_specificity <- train_confusion_matrix[1, 1] / sum(train_confusion_matrix[1, ])
print(paste("Specificity on training data:", train_specificity)) 

train_precision <- train_confusion_matrix[2, 2] / sum(train_confusion_matrix[, 2])
train_recall <- train_sensitivity
train_f1_score <- 2 * (train_precision * train_recall) / (train_precision + train_recall)
print(paste("F1-score on training data:", train_f1_score))
```

###Predictions and probability projection
```{r}
data_complet_merged_sap$predictions <- predictions
data_complet_merged_sap_longlat <- st_transform(data_complet_merged_sap, crs = "+proj=longlat +datum=WGS84")
library(leaflet)
leaflet(data_complet_merged_sap_longlat) %>%
  addTiles() %>%
  addCircles(weight = 1, radius = 100, color = "blue", fillOpacity = 0.5, 
             fillColor = ~colorQuantile("Blues", predictions)(predictions),
             popup = ~paste("Prediction:", predictions)) %>%
  addLegend(position = "bottomright", pal = colorQuantile("Blues", predictions), 
            values = predictions, title = "Predictions") %>%
  addScaleBar(position = "bottomleft")



data_complet_merged_sap<-st_transform(data_complet_merged_sap, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")

data_complet_merged_sap$Sapajus<-as.factor(data_complet_merged_sap$Sapajus)
prediction <- function(model, data_complet_merged_sap, sefit = FALSE) {
  v <- predict(model, data_complet_merged_sap, se.fit = sefit, type = "response")
  if (sefit == FALSE) {
    p <- as.vector(v)
  } else {
    p <- as.vector(v$se.fit)
  }
  p
}

points_raster<-crop(x = points_raster, y = precipitation1)
raster_stack <- stack(ndvi2, temperature1, precipitation1, vapeur1, ra1, pop1, land, vege2, vcf_non_vege, vcf_tree, points_raster)
names(raster_stack )<- c("moy_ndvi_sap", "moy_temp_sap", "moy_prec_sap", "moy_vap_sap", "moy_ra_sap", "moy_pop_sap", "moy_land_sap", "moy_vege2_sap", "moy_vcf_non_vege_sap", "moy_vcf_tree_sap", "min_distances")

carte_prediction_RF <- 1-raster::predict(model=model, object=raster_stack, type="prob")
plot(carte_prediction_RF)
plot(aireetude[0],add=T)

pred <- raster::predict(model, object = raster_stack, type = "prob", index=2)
ecarts_types <- sqrt(pred * (1 - pred))
```


###Resistance
```{r}
raster_prediction <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1_sap.tif")
raster_prediction[is.na(raster_prediction[])] <- 0
plot(raster_prediction)
writeRaster(raster_prediction,"C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1_sapsansNA.tif" )

raster_prediction_sd_sap <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1sd_sap.tif")
raster_prediction_sd_sap[is.na(raster_prediction_sd_sap[])] <- 0.5
writeRaster(raster_prediction_sd_sap,"C:/Users/bambi/OneDrive/Bureau/article/raster_prediction_sd_sapsansNA.tif" )

carte_prediction<-(raster_prediction)*100
carte_prediction[is.na(carte_prediction[])] <- 0
carte_prediction<-floor(carte_prediction)

isBecomes <- cbind(c((0:10), (11:20), (21:30), (31:40), (41:50), (51:60), (61:70), (71:80), (81:90), (91:100)), c(rep(100,11), rep(90,10), rep(80,10),rep(70,10), rep(60,10), rep(50,10), rep(40,10),rep(30, 10), rep(20,10) , rep(10,10)))
carte_resistance<- reclassify(carte_prediction, rcl = isBecomes)
carte_resistance<-as.factor(carte_resistance)

r_8bits <- writeRaster(carte_resistance,filename="C:/Users/bambi/OneDrive/Bureau/article/carte_resistance_RF_sap1.tif", datatype="INT1U")

tampons_fusion<-st_read("C:/Users/bambi/OneDrive/Bureau/article/tampons_sap.shp")
habitat<-st_read("C:/Users/bambi/OneDrive/Bureau/article/RF2_sap", 'patches')
raster_resistance <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_resistance_RF_sap1.tif")

route<-st_read("C:/Users/bambi/OneDrive/Bureau/cours M2/stage/qgis test/script données stage M2/bassin/vecteur_route_pb.shp")

routes_raster <- rasterize(route, raster_resistance, field = 1)
pixels_routes <- raster_resistance * routes_raster
plot(pixels_routes)


modified_raster_route <- reclassify(pixels_routes, cbind(from = c(9:100), to =105))
mon_modified_masked_raster <- reclassify(modified_raster_route, cbind(from = c(0:80), to = NA))
rc <- cbind(NA, 0)
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)
rc <- cbind(105, 100)
final <- reclassify(r_zero,rc, updatevalue=TRUE)


okok<-final
okok[final  <= 90] <- raster_resistance[final == 0]

mon_masked_raster <- mask(okok, tampons_fusion)
mon_habitat_raster<-mask(okok, habitat)

mon_modified_masked_raster <- reclassify(mon_masked_raster, cbind(from = c(2:100), to = c(101, NA)))
rc <- cbind(c(NA, 1), c(0, 101))
r_zero <- reclassify(mon_modified_masked_raster , rc, updatevalue=TRUE)

mon_modified_masked_raster_habitat <- reclassify(mon_habitat_raster, cbind(from = c(2:100), to = c(102)))
rc <- cbind(c(NA, NA), c(0, 0))
r_zero_habitat <- reclassify(mon_modified_masked_raster_habitat, rc, updatevalue=TRUE)

final<-r_zero_habitat
final[r_zero_habitat  <= 100] <- okok[r_zero_habitat == 0]

resultat <- r_zero
resultat[r_zero   <= 100] <- final[r_zero  == 0]
```

###Surface area
```{r}
raster_resistance9<-raster("C:/Users/bambi/OneDrive/Bureau/article/carte_resistance9.tif")
final_presence <- rasterToPolygons(raster_resistance9, dissolve = TRUE)
final_presence1<-st_as_sf(final_presence, "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs" )
niveaux <- "102"
shapefile_102 <- subset(final_presence1, raster_resistance9 =="102")

st_area(shapefile_102)
niveaux101 <- "101"
shapefile_101 <- subset(final_presence1, carte_resistance9 == niveaux101)
st_area(shapefile_101)
```

### Id vegetation
```{r}
vegetation$leg_carga<-as.factor(vegetation$leg_carga)

levels(vegetation$leg_carga)[data_complet_merged_sap$moy_vege2_sap]
View(vegetation)
plot(data_complet_merged_sap$moy_vege2_sap)
cbind(levels(vegetation$leg_carga),seq(1,length(levels(vegetation$leg_carga))))
```


###Calculation of the probability of species presence based on forest fragments
```{r}
carte_prediction_RF<-raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1_sap.tif")
carte_prediction_RF[is.na(carte_prediction_RF[])] <- 0

extraction<-raster::extract(x = carte_prediction_RF, y = aireetude, df = TRUE)
extraction <- na.omit(extraction)

extraction_sap<-tapply(extraction $carte_prediction_RF, extraction$ID, mean)
values <- c(0.2553909, 0.3369264, 0.2693827, 0.000000, 0.4676749, 0.4261511, 0.3341780, 0.1527778, 0.6071397, 0.4765852, 0.1979918, 0.4636287, 0.4357531)
data_extraction <- data.frame(ID = 1:13, Value = values)

localites <- c("Imbaú III", " Aldeia I", "Aldeia II", "Cabo Frio", "Dourada", "Gaviões", "Imbaú I", "Imbaú II", "Poço das Antas", "Pirineus", "Rio Vermelho", "União I", "União II")
data_extraction$Localites <- factor(localites, levels = localites)

data_extraction
Localities<-data_extraction$Localites
Mean<-data_extraction$Value
ggplot(data_extraction, aes(x = Localities, y = Mean)) +
  geom_bar(stat = "identity", fill = "darkgreen")+
  theme(panel.background = element_rect(fill = "white"))

```

#Correlation of predictions for the two species
```{r}
sapajus_predi <- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1_sap.tif")
sapajus_predi[is.na(sapajus_predi[])] <- 0

ouis_predi<- raster("C:/Users/bambi/OneDrive/Bureau/article/carte_prediction_RF1.tif")

ouis_predi<-crop(x=ouis_predi, y=precipitation1)
ouis_predi[is.na(ouis_predi[])] <- 0
sapajus_predi<-crop(x=sapajus_predi, y=precipitation1)


val_sap <- getValues(sapajus_predi)
val_ouis <- getValues(ouis_predi)


cor.test(val_sap, val_ouis)

corr_matrix <- corLocal(sapajus_predi, ouis_predi, test=TRUE)

plot(val_sap, val_ouis, 
     xlab = "Probability of presence of capuchins ", 
     ylab = "Probability of presence of marmosets ",
     pch = 19,
     col = "darkgreen",
     cex = 0.4)

abline(lm(val_ouis ~ val_sap), col = "red")

legend("topright", legend = "Linear regression", col = "red", lwd = 1)
```