SCRIPT R - STAGE PGM 2019

#### INITIALISATION ####
setwd("...:/ESPACE DE TRAVAIL/...")
rm(list=ls())
library(adehabitatHR)

## Loading required package: sp

## Loading required package: deldir

## deldir 0.1-16

## Loading required package: ade4

## Loading required package: adehabitatMA

## Loading required package: adehabitatLT

## Loading required package: CircStats

## Loading required package: MASS

## Loading required package: boot

library(rgdal)

## rgdal: version: 1.3-9, (SVN revision 794)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/R/packages/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/R/packages/rgdal/proj
##  Linking to sp version: 1.3-1

library(maptools)

## Checking rgeos availability: TRUE


#### REMPLIR PAR UTILISATEUR ####
for (i in 1:6) {
id<-i
SD<-4
reel<-1 #1 pour réel ou 2 pour tirages aléatoires
ecrire<-2 #1 pour oui (écrit les domaines vitaux dans shp + IMAGES) / 2 pour non

#### OBJETS ####
l<-NULL

if (id==6) {
data<-read.table("data6.csv", h=T, sep=";")
title="Émetteur 6 "
pos<-"topright"
l$x <- 877000
l$y <- 6282400
ep<-0.02
}
if (id==5) {
data<-read.table("data5.csv", h=T, sep=";")
title="Émetteur 5"
pos<-"topleft"
l$x <- 879200
l$y <- 6281160
ep<-0.01
}
if (id==4) {
data<-read.table("data4.csv", h=T, sep=";")
title="Émetteur 4"
pos<-"topleft"
l$x <- 878600
l$y <- 6281200
ep<-0.04
}
if (id==3) {
data<-read.table("data3.csv", h=T, sep=";")
title="Émetteur 3"
pos<-"topleft"
l$x <- 877000
l$y <- 6281750
ep<-0.05
}
if (id==2) {
data<-read.table("data2.csv", h=T, sep=";")
title="Émetteur 2"
pos<-"bottomleft"
l$x <- 879800
l$y <- 6280700
ep<-0.02
}
if (id==1) {
data<-read.table("data1.csv", h=T, sep=";")
title="Émetteur 1"
pos<-"bottomleft"
l$x <- 876700
l$y <- 6281300
ep<-0.02
}

LIGNES<-length(data$DATE)

AZ1<-data$AZ1
AZ2<-data$AZ2
X1<-data$X1
Y1<-data$Y1
X2<-data$X2
Y2<-data$Y2

################ CALCUL COORDONNEES REELLEES ################
if (reel==1) {
REL<-data.frame(matrix(NA,ncol=4,nrow=LIGNES))
names(REL)<-c("BETA_AZ1","BETA_AZ2","X","Y")

for (i in (1:LIGNES)) { #calcul des angles beta
  REL[i,1]<-(90-data$AZ1[i])*(pi/180)
  REL[i,2]<-(90-data$AZ2[i])*(pi/180)
}

for (i in (1:LIGNES)) { #calcul des coordonnées
  REL[i,3]<-(X1[i]*tan(REL$BETA_AZ1[i])-X2[i]*tan(REL$BETA_AZ2[i])+Y2[i]-Y1[i])/(tan(REL$BETA_AZ1[i])-tan(REL$BETA_AZ2[i]))
  REL[i,4]<-abs((((X2[i]-X1[i])*tan(REL$BETA_AZ1[i])*tan(REL$BETA_AZ2[i]))-Y2[i]*tan(REL$BETA_AZ1[i])+Y1[i]*tan(REL$BETA_AZ2[i]))/(tan(REL$BETA_AZ1[i])-tan(REL$BETA_AZ2[i])))
}
}

################ TIRAGE ALEATOIRE ################
#Tire des azimuths selon loi normale et arrondis à l'unité
if (reel==2) {
FIN<-data.frame(matrix(NA,ncol=6,nrow=LIGNES))
names(FIN)<-c("AZ1","AZ2","BETA_AZ1","BETA_AZ2","X","Y")

for (j in (1:500)){
for (i in (1:LIGNES)) {
  FIN[i,1]<-abs(round(rnorm(1,AZ1[i],SD),0))
  FIN[i,2]<-abs(round(rnorm(1,AZ2[i],SD),0))
  
}
}

for (i in (1:LIGNES)) { #calcul des angles beta
FIN[i,3]<-(90-FIN$AZ1[i])*(pi/180)
FIN[i,4]<-(90-FIN$AZ2[i])*(pi/180)
}

for (i in (1:LIGNES)) { #calcul des coordonnées
  FIN[i,5]<-(X1[i]*tan(FIN$BETA_AZ1[i])-X2[i]*tan(FIN$BETA_AZ2[i])+Y2[i]-Y1[i])/(tan(FIN$BETA_AZ1[i])-tan(FIN$BETA_AZ2[i]))
  FIN[i,6]<-abs((((X2[i]-X1[i])*tan(FIN$BETA_AZ1[i])*tan(FIN$BETA_AZ2[i]))-Y2[i]*tan(FIN$BETA_AZ1[i])+Y1[i]*tan(FIN$BETA_AZ2[i]))/(tan(FIN$BETA_AZ1[i])-tan(FIN$BETA_AZ2[i])))

}
}

#write.table(FIN,"FIN.csv",sep=";") #possibilité d'enregistrer la table créée


################ CONVERSION EN LAMBERT 93 ################
if (reel==1) {
  x<-REL$X
  y<-REL$Y
  title1<-"Positions des points réels"
}
if (reel==2) {
  x<-FIN$X
  y<-FIN$Y
  title1<-"Positions des tirages aléatoires"
}

xy<-data.frame(x,y)
coordinates(xy) <- xy[, c('x', 'y')]
plot(xy, axes=TRUE, main=title1, cex.axis=.95) #position des points

proj4string(xy)<-CRS("+init=epsg:4326") # WGS 84
CRS.new<-CRS("+init=epsg:2154") #L93
data.newcrs<-spTransform(xy, CRS.new)
#plot(data.newcrs, axes=TRUE, main="Conversion", cex.axis=.95) #vérification de la conversion
lambert<-data.frame(data.newcrs@coords,id)
lambert_plot<-data.frame(data.newcrs@coords,id)
#write.table(lambert,"lambert.csv", sep=";")
coordinates(lambert_plot) <- lambert_plot[, c('x', 'y')]

################ CALCUL DU DOMAINE VITAL - MCP ################
mcp_100 <- mcp(lambert_plot[, 3], percent = 100)
mcp_95 <- mcp(lambert_plot[, 3], percent = 95)
mcp_90 <- mcp(lambert_plot[, 3], percent = 90)
mcp_80 <- mcp(lambert_plot[, 3], percent = 80)
mcp_70 <- mcp(lambert_plot[, 3], percent = 70)
mcp_50 <- mcp(lambert_plot[, 3], percent = 50)
mcp_15 <- mcp(lambert_plot[, 3], percent = 15)

#### CENTROIDE + PLOT ####
#x11()
centroid <- apply(lambert[, 1:2], 2, mean)
plot(lambert$x, lambert$y, asp = 1, xlab="X", ylab="Y", main=paste("Position des points de l'émetteur",id))
points(centroid[1], centroid[2], pch = 19, col = 'red')
if (reel==1 && ecrire==1){
dev.copy(png,paste("./IMAGES/centroide-emet",id,".png"))
dev.off()
}
#on sélectionne les couleurs
col100<-"grey"
col95<-'#777777'
col90<-'#777777'
col80<-'#444444'
col70<-'#222222'
col50<-'#444444'
col15<-'blue'

#x11()
plot(mcp_100, axes=F,col=col100)#, main=paste("Probabilité de voir l'",title,"à un instant t"))
#col = '#999999'
axis(1)
axis(2)
plot(mcp_95, add = TRUE, col = col95)
plot(mcp_50, add = TRUE, col = col50)
#plot(mcp_70, add = TRUE, col = col70)
#plot(mcp_15, add = TRUE, col = col15)
plot(lambert_plot, cex = .75, pch = 19, col = 'red', add = TRUE)
points(centroid[1], centroid[2],cex=1, pch = 7, col = 'green')
legend(pos,title=title, pch=c(15,15,15,7,19), 
       c(paste("100% : ",round(mcp_100$area,2)," ha"),
         paste("95% : ",round(mcp_95$area,2)," ha"),
         paste("50% : ",round(mcp_50$area,2)," ha"),
         #paste("70% : ",round(mcp_70$area,2)," Ha"),
         "Centroïde","Point de position"), 
       col=c(col100,col95,col50,"green","red"))
#affichage de l’échelle
SpatialPolygonsRescale(layout.scale.bar(ep),offset=c(l$x,l$y),scale=100, fill=c("black"),plot.grid=F)
text(l$x+100/2,l$y,paste("100m","\n\n",sep=""))
if (reel==1 && ecrire==1){
  dev.copy(png,paste("./IMAGES/MCP/MCP-emet",id,".png"))
  dev.off()
}

plot(mcp_95) #écrit dans un fichier shape plus bas
mcp_95

#x11()
mcparea <- mcp.area(lambert_plot[, 3],percent=seq(50, 100, by = 5))
if (reel==1 && ecrire==1){
  dev.copy(png,paste("./IMAGES/MCP/MCP-mcparea-emet",id,".png"))
  dev.off()
}
mcparea

################ SORTIES ################
#CREATION FICHIER QGIS :
if (reel==1 && ecrire==1){
  writePolyShape(mcp_95, paste("./SHAPES/homerange-MCP-emet",id))
  #writePolyShape(homerange, paste("./SHAPES/homerange-kernel-emet",id))
}

if (reel==2 && ecrire==1){
  writePolyShape(mcp_95, paste("./SHAPES/SD-homerange-MCP-emet",id))
}



}
#### IMAGES ####

SCRIPT-STAGE PGM 2019.R

Vendredi 14 juin 2019 - 12:00:00 - created by Adrien CHARBONNEAU