Resource selection function

#Importer base de données poisson
load("Fish72.Rdata")
fish<-Fish72; rm(Fish72)

#Importer le raster du lac pour représenter graphiquement les localisations du poisson sur a carte
load("Bathymetry_Ledoux.Rdata")

#Regarder la structure du tableau de données poisson
head(fish)

#Représenter le raster
image(ledoux_depth,col=palette,ylim=c(780,1150),main="The 596 locations of the Fish recorded in 4h")

#Ajouter les localisations de poisson sur la carte
points(fish$X,fish$Y,pch=16,cex=0.7)

#Représenter le raster
image(ledoux_depth,col=palette,ylim=c(780,1150),main="Step comparison")

#Ajouter les pas observés sur la carte (longueur plus angle)
lines(fish$X[c(1:300,400,500,590)],fish$Y[c(1:300,400,500,590)],lty=4)
lines(fish$X[c(300,400)],fish$Y[c(300,400)],lwd=2)
lines(fish$X[c(400,500)],fish$Y[c(400,500)],lwd=2)


#Ajouter les pas simulés sur la carte
lines(c(fish$X[300],800),c(fish$Y[300],1050),col="red",lwd=1.5)
lines(c(fish$X[300],840),c(fish$Y[300],970),col="red",lwd=1.5)
lines(c(fish$X[300],950),c(fish$Y[300],1025),col="red",lwd=1.5)

lines(c(fish$X[400],750),c(fish$Y[400],1050),col="forestgreen",lwd=1.5)
lines(c(fish$X[400],700),c(fish$Y[400],1070),col="forestgreen",lwd=1.5)
lines(c(fish$X[400],790),c(fish$Y[400],1055),col="forestgreen",lwd=1.5)

# La fonction as.ltraj permet de transformer une matrice avec le temps et les lieux en une succession de pas, avec la longueur du pas, l'angle entre un pas et le précédent, la durée du pas etc.

#NB : lorsque l'on a plusieurs poissons ou périodes échantillonnés, cette fonction permet de le traiter en précisant le nom de chaque individu et le burst (c'est à dire une succession de relocalisation pour un animal précis à une date précise). Il en résultera une liste contenant toutes les trajectoires calculées.

track_fish=as.ltraj(fish[,c("X","Y")],fish$posix_est,id=fish$Name,
                     burst=fish$burst,slsp="missing",proj4string = crs("+proj=aeqd +lat_0=46.802927 +lon_0=-73.277817 +x_0=1000 +y_0=1000 +datum=WGS84 +units=m +no_defs"))

# Voyons ce que nous avons. Ici, on a un poisson suivi pendant quatre heures, on a juste un burst donc un élément dans la liste.
head(track_fish[[1]])

#La fonction redisltraj nous permet d'homogénéiser chaque pas à des pas d'une minute en spécifiant la durée du pas en secondes. Il recalculera simplement les emplacements de chaque pas par interpolation.

fixed_track <- redisltraj(track_fish, 60, type="time")

head(fixed_track[[1]])

#Représenter la carte bathymétrique en fond, avec un grossissement sur la zone utilisée
image(ledoux_depth,col=palette,ylim=c(800,1100),xlim=c(700,1200),main="Raw location vs. Fixed track")

#Ajouter les localisations du poisson
points(fish$X,fish$Y,pch=16,cex=0.7)

#Ajouter la trajectoire calculée du poisson suivi
lines(fixed_track[[1]]$x,fixed_track[[1]]$y,col="salmon",lwd=1.5)

# Définir la profondeur minimale où le poisson est capable d'aller
lim_depth=min(fish$Depth,na.rm=T)
print(lim_depth)

# Créer un raster avec seulement les zones du lac avece ce minimum de profondeur

# Pour ce faire, j'ai inélégamment transformer le raster de bathymétrie en un tableau.
map_temp=as.data.frame(ledoux_depth, row.names=NULL, optional=FALSE, xy=TRUE,
                     na.rm=FALSE, long=FALSE)

# Puis j'ai sélectionner seulement les endroits du lac suffisamment profonds pour accueillir le poisson observé.

map_temp=map_temp[which(map_temp$layer>=lim_depth),]

# Je crée un raster vide des même dimensions que les raster original du lac
 e <- extent(c(min(map_temp$x),max(map_temp$x),min(map_temp$y),max(map_temp$y)))
  ncellx=(e[2]-e[1])
  ncelly=(e[4]-e[3])
  empty_r<- raster(e, ncol=ncellx, nrow=ncelly)

# Finalement j'intègre les endroits sélectionnés dans le raster vide
r <- rasterize(map_temp[, 1:2], empty_r, map_temp[,3], fun=mean)

# J'ai construit une fonction de contrainte qui doit avoir un objet de type SpatialPixelsDataFrame donc nous convertissons le raster
map <- as(r, "SpatialPixelsDataFrame")
crs(map)<-NA

# Ici, la zone du lac où nous allons autoriser la simulation de trajectoires, dépendamment de la distribution en profondeur du poisson observé
image(ledoux_depth,col="black",ylim=c(780,1150),xlim=c(300,1600))
par(new=T)
image(r,ylim=c(780,1150),col=palette,xlim=c(300,1600))

# fonction de contrainte qui supprime l'emplacement si les coordonnées d'une trajectoire simulée (x) tombaient en dehors de la zone contrainte (par). Cette fonction de contrainte est l'un des arguments intégrés à la fonction simulant les trajectoires.

consfun <- function(x, par){
  coordinates(x) <- x[,1:2]
  ov <- over(x, geometry(par))
  return(all(!is.na(ov)))
}

  mo <- NMs.randomCRW(na.omit(fixed_track), rangles=TRUE, rdist=TRUE,
                      constraint.func=consfun,
                      constraint.par = map, nrep=9,fixedStart=TRUE)

 Simulated_tracks=testNM(mo)

# Transformer la liste des trajectoires en data frame, en ajoutant une colonne "no" pour l'id de la trajectoire ; 1 étant celui observé, et 2 à 10, les 9 simulations
no=1
RSF_database=cbind(fixed_track[[1]][,c(1,2)],no)

for (i in 1:9){
  no=i+1
  RSF_database=rbind(RSF_database,cbind(Simulated_tracks[[1]][[i]][,c(1,2)],no))
}

# Un vecteur de couleur pour avoir la trajectoire observée en couleur et les simulées en gris.
couleurs=c("salmon",gray.colors(9, start = 0.3, end = 0.9, gamma = 2.2, rev = FALSE))

# Représenter la carte
image(ledoux_depth,col=palette,ylim=c(780,1150),main="Simulated and observed trajectories")

#Ajouter toutes les lignes de trajectoire
for (i in 10:1){
  temp=subset(RSF_database,no==i)
  lines(temp$x,temp$y,col=couleurs[i])
}

distri_depth=fish$Depth
hist(distri_depth,main="",xlab="Depth (m)")

##Passons rapidement sur le recalcul de la profondeur réelle pour la trajectoire observée

fixed_track[[1]]$depth=NA

#Boucle for pour assigner une profondeur à chaque localisation dans la trajectoire observée
for (i in 1:nrow(fixed_track[[1]])){

  # le temps à la localisation i
  time_stamp=as.numeric(substr(fixed_track[[1]]$date[i],12,13))*3600+
                  as.numeric(substr(fixed_track[[1]]$date[i],15,16))*60+
                  as.numeric(substr(fixed_track[[1]]$date[i],18,19))

  # Regrouper les temps de toutes les mesures de profondeur
  time_table=as.numeric(substr(fish[,"posix_est"],12,13))*3600+
                  as.numeric(substr(fish[,"posix_est"],15,16))*60+
                  as.numeric(substr(fish[,"posix_est"],18,19))

  # Selectionner la profondeur où le temps de la mesure est la plus proche du temps de la localisation i             
  ind_timing1=which(abs(time_stamp-time_table)==min(abs(time_stamp-time_table),na.rm=T))
          print(c(i,length(ind_timing1)))

  # Si les deux (précédent et suivant) mesures de profondeur sont disantes de manière égale dans le temps, on moyenne les profondeurs entre les deux pour avoir une profondeur à la localisation i
   if(length(ind_timing1)==2){
                fixed_track[[1]]$depth[i]=mean(fish[,"Depth"][ind_timing1])
   }

  # Si on a qu'une mesure de profondeur la plus proche
   if(length(ind_timing1)==1){

  # Si la mesure de profondeur la plus proche est après le temps de localistion i, on prend cette mesure et celle avant pour ensuite faire une interpolation linéaire pour retrouver la profondeur à la localisation i
    if(time_stamp-time_table[ind_timing1]<0){
      ind_timing2=ind_timing1-1
      time=time_table[ind_timing1]-time_table[ind_timing2]
      depth_var=fish[ind_timing1,"Depth"]-fish[ind_timing2,"Depth"]
      depth_inter=(time_stamp-time_table[ind_timing2])*(depth_var/time)+fish[ind_timing2,"Depth"]
      fixed_track[[1]]$depth[i]=depth_inter
    }

  # Si la mesure de profondeur la plus proche est avant le temps de localistion i, on prend cette mesure et celle après pour ensuite faire une interpolation linéaire pour retrouver la profondeur à la localisation i
     if(time_stamp-time_table[ind_timing1]>0){
      ind_timing2=ind_timing1+1
      time=time_table[ind_timing2]-time_table[ind_timing1]
      depth_var=fish[ind_timing2,"Depth"]-fish[ind_timing1,"Depth"]
      depth_inter=(time_stamp-time_table[ind_timing1])*(depth_var/time)+fish[ind_timing1,"Depth"]
      fixed_track[[1]]$depth[i]=depth_inter

    }
  # Si la mesure de profondeur la plus proche est au temps exact de la localisation i, on prend directement cette mesure pour assigner une profondeur à la localisation i
    if(time_stamp-time_table[ind_timing1]==0){
      fixed_track[[1]]$depth[i]=fish[ind_timing1,"Depth"]
    }
  }
}

  # Créer une colonne pour recevoir la profondeur de la colonne d'eau pour chaque localisation simulée
RSF_database$max_depth=NA

  # Utilise la fonction extract du package raster pour récupérer la profondeur de la colonne d'eau pour chaque localisation simulée
RSF_database$max_depth=extract(ledoux_depth,RSF_database[,c(1:2)])

  # Vecteur vide
depth_random=NA
  #Boucle for où on colle chaque échantillonnage aléatoire de profondeur dans le vecteur
  for (i in which(RSF_database$no!=1)){
    depth_random=c(depth_random,sample(distri_depth[which(distri_depth<=RSF_database$max_depth[i])],1))
  }
  #Ensuite, dans la base de données regroupant toutes les trajectoires observée et simulées, on ajoute la colonne profondeur

RSF_database$depth=c(fixed_track[[1]]$depth,depth_random[-1])

#Voyons comment se comporte notre distribution en profondeur

ggplot(RSF_database,aes(x=depth))+
  geom_histogram(position="dodge",bins=15)+
  facet_wrap(~no)

  # Créer un vecteur d'intervalle de profondeur
depth_vector=seq(-0.5,16.5,1)

  # Charger la base de données zooplancton
load("large_Z.Rdata")

  # Créer la colonne avec les futures concentrations en zooplancton
RSF_database$LZ=NA

  # Créer un raster vide pour créer un raster zooplancton qui sera extrait pour chaque couche de profondeur
e <- extent(ledoux_depth)
  #Résolution horizontale souhaitée : 5m par 5m
ncellx=(e[2]-e[1])/5
ncelly=(e[4]-e[3])/5
r<- raster(e, ncol=ncellx, nrow=ncelly)

  # Boucle for pour assigner les concentrations en zooplanction à la fois aux localisations simulées mais aussi observées, pour chaque couche de profondeur
for (prof in 1:(length(depth_vector)-1)){

  # Sous-échantillon de la base de données trajectoires avec seulement les localisations dans la bonne couche de profondeur
  temp=RSF_database[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1]),]

  # Si on a des localisations dans cette couche
    if(nrow(temp)>0){

  # Selection des mesures de zooplancton pour la couche correspondante
      tempZ=Large_zooplankton[which(Large_zooplankton$Depth==(depth_vector[prof]+0.5)),c("X","Y","Predicted")]

  # Transformer en un raster
      temp_rast_LZ <- rasterize(tempZ[, 1:2], r, tempZ[,3], fun=mean)

  # Extraire les concentrations de zooplancton pour toutes les localisations dans cette couche de profondeur
      RSF_database$LZ[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1])]=extract(temp_rast_LZ,RSF_database[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1]),c(1,2)])
    }
}

head(RSF_database)

RSF_database$weight=ifelse(RSF_database$no==1,1,5000)

  #Pour les régression conditionnelle, on a besoin d'une variable VRAI/FAUX. Ici, seule la trajectoire no°1 est VRAIE puisque c'est l'observée
RSF_database$case=RSF_database$no==1

  # Mettre à l'échelle les variables
RSF_database$LZ=scale(RSF_database$LZ)
RSF_database$max_depth=scale(RSF_database$max_depth)

  # Ecrire le modèle de régression
RSF_model= glm(case ~ max_depth + LZ,data = RSF_database,weight=weight, family = binomial(link = "logit"))

summary (RSF_model)

RSS_zoo=exp(RSF_model$coefficients[3])
paste0("RSS = ",round(RSS_zoo,3))

  # Ici une base de données contenant des estimations de régression calculées 10 fois avec différents nombres de trajectoires simulées pour nos poissons
N_Simul=read.csv("S4_base.csv")
N_Simul_Z=N_Simul[which(N_Simul$param=="Large zooplankton"),]


ggplot(N_Simul_Z, aes(y=coeff,group=N_random_tracks))+
  geom_boxplot(outlier.shape = NA)+
  theme_bw()+
  xlab("Number of simulated trajectories")+
  scale_x_continuous(breaks=c(seq(-0.4,0.4,0.8/6)),
                   labels=as.character(seq(0,300,50)))+
  geom_vline(xintercept=0,lty=2)+
  theme(axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("RSF estimate"))

  # Impoter la base de données où nous avons calculé 500 fois la régression basée sur une trajectoire simulée et 150 autres trajectoires simulées

H0_500=read.csv("HO_500.csv")

    # Impoter la base de données où nous avons calculé un nombre de fois variable la régression basée sur une trajectoire simulée et 150 autres trajectoires simulées

H0_data= read.csv("S5DB.csv",row.names=1)

ggplot(H0_500, aes(x=RSF_estimate))+
  geom_histogram(fill="gray80",colour="gray50",bins=18,aes(y=(..count../sum(..count..)) * 100))+
  theme_bw()+
  geom_vline(xintercept = quantile(H0_500$RSF_estimate,c(0.025,0.975)),lty=2)+
  xlab("RSF estimate")+
  theme(axis.title.y = element_text(size=14),axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("Percentage of occurrence (%)"))

  ggplot(H0_data, aes(y=inf,x=N_H0_RSF))+
  geom_smooth(col="gray50")+


  geom_smooth(aes(y=sup,x=N_H0_RSF),col="gray50")+
  geom_point()+
  geom_point(aes(y=sup,x=N_H0_RSF))+
  theme_bw()+
  geom_vline(xintercept = 500,lty=2)+
  xlab("Number of RSF iterations")+
  theme(axis.title.y = element_text(size=14),axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("Boundaries of the 5% confidence interval of RSF estimates"))