#########################################################################
# Cree le 14/12/2016 par Anne Mouget
# Derniere mise a jour le 09/05/2017 par Anne Mouget

# Objectif : reunir dans un meme script toutes les etapes de traitement des donnees
# Ce script propose un script pour le traitement des donnees issues de l'echososndeur RoxAnn

# Pour l'utilisation du script, se referrer au protocole
# "Protocole d’utilisation du système RoxAnn© pour la classification des fonds lacustres" Mouget et al., 2017

# Etape 0 : Recuperation des donnees utiles, dialogue avec l'operateur
# Etape 1 : Chargement et installation des packages 
# Etape 2 : Chargement et nettoyage des donnees
# Etape 3 : Preparation du fond de carte
# Etape 4 : Affichage (et enregistrement) de la carte de bathymetrie
# Etape 5 : Affichage de la carte des substrats
#	# Etape 5.1 : validation du modele et application a l'ensemble des donnees
#	# Etape 5.2 : interpolation des donnees de probabilite de presence
#	# Etape 5.3 : compilation des informations de presence
#	# Etape 5.4 : preparation des donnees a l'affichage
#	# Etape 5.5 : affichage (et enregistrement) de la carte 
# Etape 6 : aide a l'edition des carres RoxAnn
#	# Etape 6.1 : rectangles simples
#	# Etape 6.2 : croisements entre carres
#	# Etape 6.3 : nettoyage
#	# Etape 6.4 : reorganisation et edition des coordonnes
#	# Etape 6.5 : affichage
# Etape 7 : Creation et affichage d'un fichier recapitulatif

# Modele choisit : Random Forest
# version 4.1

#########################################################################

# Nettoyage de l'espace de travail R avant de commencer

rm(list=ls()) 
graphics.off() 

#########################################################################################
#########################################################################################

# Etape 0 : Dialogue avec l'operateur pour la selection de donnees

#########################################################################
#########################################################################
# La discussion avec l'operateur necessite l'installation (automatique) d'un package
chooseCRANmirror(ind=1)
if (is.element('tcltk', installed.packages()[,1]) == FALSE) {install.packages('tcltk') }
library(tcltk)

# Deux fonctions sont crees pour simplifier le dialogue avec l'operateur
attente <- function(texte) {
    tt <- tktoplevel()
    tkpack( tkbutton(tt, text='Continuer', command=function()tkdestroy(tt)),
        side='bottom')
    tkbind(tt,'<Key>', function()tkdestroy(tt) )
	fr1=tkframe(tt)
	tkpack(fr1)
	t1=tklabel(fr1, text=texte)
	tkpack(t1)
	tkwait.window(tt)
}

reponse <- function(texte) {
	base=tktoplevel()
	fr1=tkframe(base)
	tkpack(fr1)
	t1=tklabel(fr1, text=texte)
	tkpack(t1)
	fr3=tkframe(base)
	tkpack(fr3)
	v_prop_n=tclVar("")
	e1=tkentry(fr3, textvariable=v_prop_n)
	tkpack(e1,side="right") 
	ok=function(){
		chiffre <<- tclvalue(v_prop_n)
		tkdestroy(base)
	}
	okbutton=tkbutton(fr3, text="OK", command=ok)
	tkpack(okbutton) 
	tkwait.window(base)
	return(chiffre)
}


# Dans un meme repertoire doit se trouver :
# les 2 bases de donnees correctement formatees (calibration et mesures)
# la base de donnee PlanEau (de BD Carthage) dezipee et avec toutes ses dependances

# Selection du fichier de calibration
attente("pour lancer la sélection du fichier de CALIBRATION, appuyer sur continuer")
nomFichierCalib <- file.choose()
nomFichierCalib <- gsub("\\\\","/",nomFichierCalib)
sep="	"
dec="."
tableCalib <- read.table(nomFichierCalib, sep=sep, head=F, dec=dec)

# Le chemin d'acces au repertoire de travail est recupere a partir du chemin du fichier de calibration
chemin <- dirname(nomFichierCalib)
setwd(dir=chemin)

# Pour vérifier que l'importation des donnees s'est correctement deroulee, il est possible d'executer un summary
#summary(tableCalib) 

# Selection des donnees de mesures non calibrees brutes
attente("pour lancer la sélection du fichier de MESURES, appuyer sur continuer")
nomFichierMesures <- file.choose()
nomFichierMesures <- gsub("\\\\","/",nomFichierMesures)
sep="	"
dec="."
tableCarto <- read.table(nomFichierMesures, sep=sep, head=F, dec=dec)

# Pour vérifier que l'importation des donnees s'est correctement deroulee, il est possible d'executer un summary
#summary(tableCarto) 

# Entrer ici la profondeur du transducteur (en metre)
profondeurTransducteur <- reponse("Entrer ici la profondeur d'immersion du transducteur (en metre) ex : 0.40")
profondeurTransducteur <- as.numeric(profondeurTransducteur)

# Entrer ici la profondeur maximale du plan d'eau (en metres)
# Cette donnee est utilisee pour supprimer les valeurs aberrantes
profondeurMax <- reponse("Entrer ici la profondeur maximale du lac en metre (ex : 40)")
profondeurMax <- as.numeric(profondeurMax)
profondeurMax <- -profondeurMax

# chargement des donnees de tous les lacs et plans d'eau francais
# Base de donnees carthage, disponible gratuitement au telechargement a l'adresse suivante :
# http://services.sandre.eaufrance.fr/telechargement/geo/ETH/PlanEau/FXX/2014/PlanEau_FXX-shp.zip
# Mettre la BDD chargee (dezipee et avec ses dependances (plusieurs fichiers aux formats varies) dans le repertoire
# Dans nomfichierCarte, indiquer le nom du fichier .shp

nomFichierCarte <- "PlanEau.shp"


# Entrer ici le nom ou partie du nom du lac d'etude 
# /!\ les accents passent mal sur R. Ne pas mettre de nom avec accent
# ex : étang d'entressen, il y a un accent sur etang donc je ne rentre que "entressen"
# si plusieurs lacs repondent a ce nom, unea recherche est automatiquement faite, basee sur les coordonnees

nomLac <- reponse("Entrer le nom du lac sans accent ni majuscule pour la recherche de carte") 

#Entrer ici le nom du lac pour le nommer dans les fichiers et les titres
nomLac2 <- reponse("Entrer le nom du lac qui sera affiche dans les titres")


#########################################################################
#########################################################################
# Etape 1

#########################################################################
# Installation si necessaire des packages

chooseCRANmirror(ind=1)
if (is.element('caret', installed.packages()[,1]) == FALSE) {install.packages('caret') }
if (is.element('maptools', installed.packages()[,1]) == FALSE) {install.packages('maptools') }
if (is.element('akima', installed.packages()[,1]) == FALSE) {install.packages('akima') }
if (is.element('stringr', installed.packages()[,1]) == FALSE) {install.packages('stringr') }
if (is.element('rgl', installed.packages()[,1]) == FALSE) {install.packages('rgl') }
if (is.element('maps', installed.packages()[,1]) == FALSE) {install.packages('maps') }
if (is.element('rgdal', installed.packages()[,1]) == FALSE) {install.packages('rgdal') }
if (is.element('e1071', installed.packages()[,1]) == FALSE) {install.packages('e1071') }
if (is.element('randomForest', installed.packages()[,1]) == FALSE) {install.packages('randomForest') }
if (is.element('spatstat', installed.packages()[,1]) == FALSE) {install.packages('spatstat') }
if (is.element('geosphere', installed.packages()[,1]) == FALSE) {install.packages('geosphere') }
if (is.element('lubridate', installed.packages()[,1]) == FALSE) {install.packages('lubridate') }


# Chargement des packages

library(caret)	# RandomForest (choix du meilleur modele)
library(e1071)	# RandomForest (construction du modele)
library(randomForest)	# RandomForest
library(maptools)	# Pour lire les format shapefile
library(akima)	# Pour l'interpolation lineaire
library(stringr)	# Pour travailler avec des chaines de caractere
library(rgl)	# Pour l'affichage 3D
library(maps)	# Pour afficher l'échelle des cartes
library(rgdal)	# Pour transformer le referentiel des cartes
library(spatstat)	# Pour savoir si les points sont bien situés dans le lac
library(geosphere)	# Pour la conversion des distances de degre a metre
library(lubridate)	# Pour les calculs de duree

#########################################################################
#########################################################################
# Etape 2

#########################################################################

# Nettoyage des bases de donnees
# Creation de tables de travail a partir des fichiers charges

#########################################################################
#########################################################################
# Preparation du fichier de calibration

# on nomme les différentes colonnes des donnees brutes
colnames(tableCalib) <- c("codeRox","Lat","Long","Prof","E1_2","E2_2","Heure","Date","Substrat")

# E1 et E2 etant eleves au carre dans le donnees brutes, on en recupere les racines
# que l'on reinjecte dans la base de donnees
E1 <- sqrt(tableCalib$E1_2)
E2 <- sqrt(tableCalib$E2_2)
tableCalib <- cbind(tableCalib, E1, E2)
tableCalib <- na.omit(tableCalib)

#########################################################################
# Preparation du fichier de cartographie (mesures sans calibration)

# on nomme les différentes colonnes des donnees brutes
colnames(tableCarto) <- c("codeRox","Lat","Long","Prof","E1_2","E2_2","Heure","Date")

# E1 et E2 etant eleves au carre dans le donnees brutes, on en recupere les racines
# que l'on reinjecte dans la base de donnees
E1 <- sqrt(tableCarto$E1_2)
E2 <- sqrt(tableCarto$E2_2)
tableCarto <- cbind(tableCarto, E1, E2)
tableCarto <- na.omit(tableCarto)

#########################################################################
# Nettoyage des jeux de donnees (suppression des donnees aberrantes)
# Suppression des lignes dont la profondeur est inférieure à la profondeur maximale
# Suppression des points avec une pente trop importante par rapport à leurs voisin

# ecartProfMax est la difference de profondeur maximale (en m)
# acceptee entre 2 points de memes coordonnees
ecartProfMax <- 10

# distMax est la distance maximale entre 2 points a comparer (en m)
distMax <- 10

# penteMax est la pente maximum autorisee
# penteMax = ecart de profondeur/distance entre les points
penteMax <- 5

#########################################################################
# Nettoyage de la table de calibration
tableCalib <- subset(tableCalib, tableCalib$Prof > profondeurMax)

ligne <- 2
while (ligne <= nrow(tableCalib)) {
	distLat <- abs(tableCalib$Lat[ligne]-tableCalib$Lat[ligne-1])
	distLong <- abs(tableCalib$Long[ligne]-tableCalib$Long[ligne-1])
	dist <- (sqrt(distLat^2 + distLong^2))*111000 	#dist la distance en m
	if (dist < distMax) {
		ecartProf <- abs(tableCalib$Prof[ligne]-tableCalib$Prof[ligne-1])
		if (dist != 0) {
			pente <- ecartProf/dist
			if (pente > penteMax) {
			tableCalib <- tableCalib[-ligne,]
			}
			else {ligne <- ligne + 1}
		}
		if (dist == 0) {
			if (ecartProf > ecartProfMax) {
				if (tableCalib$Prof[ligne] < tableCalib$Prof[ligne-1]) {
					tableCalib <- tableCalib[-ligne,]
				}
				else {
					tableCalib <- tableCalib[-(ligne-1)]
					ligne <- ligne + 1
				}
			}
			else {ligne <- ligne + 1}
		}
	}
	else {ligne <- ligne+1}	
} 

#########################################################################
# Nettoyage de la table de mesures
tableCarto <- subset(tableCarto, tableCarto$Prof > profondeurMax)

ligne <- 2
while (ligne <= nrow(tableCarto)) {
	distLat <- abs(tableCarto$Lat[ligne]-tableCarto$Lat[ligne-1])
	distLong <- abs(tableCarto$Long[ligne]-tableCarto$Long[ligne-1])
	dist <- (sqrt(distLat^2 + distLong^2))*111000 	#dist la distance en m
	if (dist < distMax) {
		ecartProf <- abs(tableCarto$Prof[ligne]-tableCarto$Prof[ligne-1])
		if (dist != 0) {
			pente <- ecartProf/dist
			if (pente > penteMax) {
			tableCarto <- tableCarto[-ligne,]
			}
			else {ligne <- ligne + 1}
		}
		if (dist == 0) {
			if (ecartProf > ecartProfMax) {
				if (tableCarto$Prof[ligne] < tableCarto$Prof[ligne-1]) {
					tableCarto <- tableCarto[-ligne,]
				}
				else {
					tableCarto <- tableCarto[-(ligne-1)]
					ligne <- ligne + 1
				}
			}
			else {ligne <- ligne + 1}
		}

	}
	else {ligne <- ligne+1}	
} 


#########################################################################
#########################################################################
# Selection des variables d'interet

# Pour la table de calibration
tableCalibRed <- tableCalib[,c("Lat","Long","E1","E2","Prof","Substrat")] 

# Pour la table de cartographie
tableCartoRed <- tableCarto[,c("Lat","Long","E1","E2","Prof")]

# On ajoute le profondeur du d'immersion du transducteur sur l'ensemble des profondeurs
tableCalibRed$Prof <- tableCalibRed$Prof - profondeurTransducteur
tableCartoRed$Prof <- tableCartoRed$Prof - profondeurTransducteur

#########################################################################
#########################################################################
# Etape 3
# Preparation du fond de carte du lac etudie

#########################################################################

# On importe les donnees en precisant qu'elles sont en lambert 93
fondCarte <- readShapePoly(nomFichierCarte, proj4string = CRS("+init=epsg:2154"))

# on trandorme les coordonnees du lambert 93 au wgs84
fondCarte <- spTransform (fondCarte, CRS ("+init=epsg:4326") )
fondCarte0 <- fondCarte

# suppression des données dont le nom (NomEntiteH) est inconnu (NA)
pos <- which(is.na(fondCarte@data [, "NomEntiteH"]))
fondCarte <- fondCarte[-pos, ]

# recherche de la (ou des) ligne(s) comprenant le nom du lac
pos <- data.frame()
for (i in 1:length(fondCarte$NomEntiteH)) {
	if (str_detect(fondCarte$NomEntiteH[i], nomLac) == T) {
		pos <- rbind(pos,i)
	}
}

# si plusieurs lacs correspondent a la recherche, on compare les resultats aux coordonnees GPS
if (dim(pos)[1] > 1) {
	for (i in 1:dim(pos)[1]) {
		lacTemp <- fondCarte[pos[i,1],]
		if ((tableCartoRed[1,1] <= lacTemp@bbox[2,2]) &
			(tableCartoRed[1,1] >= lacTemp@bbox[2,1]) &
			(tableCartoRed[1,2] <= lacTemp@bbox[1,2]) &
			(tableCartoRed[1,2] >= lacTemp@bbox[1,1])) {
				posF <- pos[i,1]
		}
	}
pos <- posF
}
if (is.null(dim(pos)[1]) != T) {pos <- pos[1,1]}

# Si aucun lac ne correspond à la recherche (nom different), les coordonnees sont comparees a l'ensemble des lacs
if (is.null(pos)) {
	for (i in 1:nrow(fondCarte0)) {
		lacTemp <- fondCarte0[i,]
		if ((tableCartoRed[1,1] <= lacTemp@bbox[2,2]) &
			(tableCartoRed[1,1] >= lacTemp@bbox[2,1]) &
			(tableCartoRed[1,2] <= lacTemp@bbox[1,2]) &
			(tableCartoRed[1,2] >= lacTemp@bbox[1,1])) {
				posF <- pos[i,1]
		}
	}
pos <- posF
}


lac <- fondCarte[pos,]

# a afficher pour s'assurer qu'il s'agit du bon lac
#windows()
#plot(lac,axes=T, main=paste("fond de carte du plan d'eau de",nomLac2)) 


# on cree la grille d'interpolation. Elle forme un carre qui englobe la totalite du lac
# t1 et t2 determinent la taille de chaque maille
# une grille fine sera precise mais demande un temps de traitement plus long
c1 <- 10^-5 #la taille des cellules en longitude (en degre)
c2 <- 10^-5 #la taille des cellules en latitude (en degre)

# l1 et l2 donnent le nombre de mailles
l1 <- (lac@bbox[1,2] - lac@bbox[1,1])/c1
l2 <- (lac@bbox[2,2] - lac@bbox[2,1])/c2

# Par defaut, les limites de la grille d'interpolation sont les limites du lac
x0 <- seq(lac@bbox[1,1], lac@bbox[1,2], length = l1)
y0 <- seq(lac@bbox[2,1], lac@bbox[2,2], length = l2)

# Mais il est possible de les changer pour les limites de la zone mesuree
#x0 <- seq(min(tableCarto$Long), max(tableCarto$Long), length = l1)
#y0 <- seq(min(tableCarto$Lat), max(tableCarto$Lat), length = l2)

#########################################################################
#########################################################################
# Etape 4 : affichage de la carte des profondeurs

#########################################################################
# Interpolation lineaire avec la package akima 

profMin <- round(min(tableCarto$Prof))

xProf <- tableCarto$Long
yProf <- tableCarto$Lat
zProf <- tableCarto$Prof
carteProf <- interp(xProf,yProf,zProf, x0, y0, linear=T, extrap=F, duplicate="mean")

# Lorsqu'une profondeur a été estimée en dehors du lac, il est effacé
for (longitude in 1:dim(carteProf$z)[1]) {
	for (latitude in 1:dim(carteProf$z)[2]) {
		if (is.na(carteProf$z[longitude, latitude]) != T) {
			if (inside.owin(carteProf$x[longitude], carteProf$y[latitude], lac) == F) {
				carteProf$z[longitude, latitude] <- NA
			}
		}
	}
}

# Affichage de la carte des profondeurs
couleurProfMin <- "cyan"
couleurProfMax <- "black"
couleurProf <- colorRampPalette(c(couleurProfMax, couleurProfMin))(-profMin+1)

# L'emplacement des legendes est regie par ces variables
xLegende <- max(x0)
yLegendeCouleur <- max(y0)
yLegendeTrajectoire <- min(y0)
yLegendeEchelle <- (9*min(y0)+max(y0))/10

# Taille de la fenetre adaptee au lac
xWindows <- (max(x0)-min(x0))*10000
yWindows <- (max(y0)-min(y0))*10000

# Affichage a proprement parler de la carte de bathymetrie
windows(xWindows, yWindows)
par(xpd=TRUE, mai=c(1,1,1,2))
image(carteProf$x, carteProf$y, carteProf$z, 
	xlim=c(min(x0), max(x0)), ylim=c(min(y0), max(y0)),
	main=paste("Carte de bathymetrie de",nomLac2),
	col=couleurProf,
	xlab="Longitude", ylab="Latitude", asp=1)
plot(lac, add=TRUE, border="black", lwd=3)
points(tableCartoRed$Lat ~ tableCartoRed$Long, cex=0.1, col="red")
legend(xLegende,yLegendeCouleur, 
	legend=seq(profMin, 0, by=-round(profMin/10)), 
	fill=colorRampPalette(c("black", "cyan"))(length(seq(profMin, 1, by=-round(profMin/10)))), 
	title="Profondeur en metres")
legend(xLegende, yLegendeTrajectoire,
	legend="trajectoire", col="red", lty=1)
map.scale(x=xLegende, y=yLegendeEchelle, relwidth = 0.1, metric = TRUE, ratio=F)

# Enregistrement de la carte aux formats pdf et png

savePlot(paste0("carteProfondeur",nomLac2), "pdf")
savePlot(paste0("carteProfondeur",nomLac2), "jpeg")

#########################################################################
#########################################################################
# Etape 5 : affichage de la carte des substrats

#	# Etape 5.1 : calibration et validation du modele

#########################################################################

# Calibration du modele sur un jeu de donnees d'apprentissage "app"
# Puis validation par verification des predictions sur un jeu de validation "val"
# Utilisation des variables explicatives E1, E2 et Prof
# Le modele choisit par defaut est random forest

# separation du jeu de calibration en 2 : apprentissage et validation
# moitie pour apprentissage, moitie pour la validation 

pourcentageApp <- 0.5
appindex = sample(1:nrow(tableCalibRed), round(nrow(tableCalibRed)*pourcentageApp), replace = FALSE)
tableCalibRed$Substrat <- factor(tableCalibRed$Substrat,exclude=NULL)
app = tableCalibRed[appindex,]
val = tableCalibRed[-appindex,]

modelTest <- train(Substrat ~ E1 + E2 + Prof, data = app, method = "rf") 

# Pour plus d'informations sur la qualite du modele choisit, il est possible d'executer les lignes ci-dessous
#print(modelTest$finalModel)

# Ces lignes donnent le taux de substrat correctement dans le jeu d'apprentissage (sur lequel s'est calibre le modele)
appBC <- mean(modelTest$finalModel$predicted==app$Substrat)
paste("le taux de substrats correctement classes dans le jeu d'apprentissage est",
	round(appBC,digits = 3))

# Application du modele au jeu de validation
pred <- predict(modelTest, val)

# Ces lignes donnent le taux de substrat correctement dans le jeu de validation
valBC <- mean(pred==val$Substrat)
paste("le taux de substrats correctement classes dans le jeu de validation est",
	round(valBC,digits = 3))

#########################################################################
# Une fois le modele valide (accepte), calibration du modele sur l'ensemble du jeu de données (tableCalibRed)

model <- train(Substrat ~ E1 + E2 + Prof, data = tableCalibRed, method = "rf") 
modelFinal <- model$finalModel

# Pour connaitre le taux de substrats du jeu de calibration correctement classes,
# il est possible d'executer les lignes ci-dessous
calibBC <- mean(model$finalModel$predicted==tableCalibRed$Substrat)
paste("le taux de substrats correctement classes dans le jeu de calibration est",
	round(calibBC,digits = 3))

#########################################################################
# Application du modele sur l'ensemble des donnees sans calibration

substratPred <- predict(modelFinal, tableCartoRed)
tableSub <- cbind(tableCartoRed, substratPred)
#summary(tableSub)

# on cree dans le repertoire de travail un fichier .csv avec les substrats predit le long de la trajectoire
# cette etape n'est pas essentiel mais permet de conserver un historique des donnees le long de la trajectoire

write.csv2(tableSub,file="substratsModele.csv", row.names=F)

#########################################################################
########################################################################## 
# Etape 5 : affichage de la carte des substrats

#	# Etape 5.2 : Interpolation des probabilites de presence

#########################################################################

tableP0 <- tableSub[,c("Lat","Long","Prof","substratPred")]
tableP0 <- na.omit(tableP0)

# On cree une matrice avec une colonne par substrat et dans chaque ligne presence ou absence de ce substrat (0/1)
# La matrice a des colonnes predefinies pour eviter les categories melange.
# Si substrat multiple (ex: sable_gravier), presence dans les deux substrats (sable=1, gravier=1)

# nSub le nombre de substrats predefinis
# tableP (table pourcentage) a pour nom de colonnes les noms des substrats predefnis + latitude et longitude

nSub <- 8
tableP <- matrix(0, nrow(tableP0),nSub+2) 
names <- c("Lat","Long","vegetaux","vase","sable","gravier","pierre","galet","bloc","IND")
colnames(tableP) <- names

for (ligne in 1:nrow(tableP)) {
	tableP[ligne,1] <- tableP0[ligne,1]
	tableP[ligne,2] <- tableP0[ligne,2]	
	for (cherch in 3:(nSub+2)) {
		if (str_detect(tableP0[ligne,4],colnames(tableP)[cherch])==TRUE) {
			tableP[ligne,cherch] <- 1
		}
		if (rowSums(tableP[,3:9])[ligne]==0) {
			tableP[ligne,"IND"] <- 1
		}
	}
}
tableP <- as.data.frame(tableP)
#summary(tableP)

#on supprime les colonnes dont la somme vaut 0

colonne <- 3
while (colonne <= ncol(tableP)) {
	if (sum(tableP[,colonne])==0) {
		tableP <- tableP[,-colonne]
	}
	else {colonne <- colonne + 1}
}
#summary(tableP)

#########################################################################
# On cree un tableau et une carte par substrat presentant le probabilite de
# presence de ce substrat a un endroit donne
# Cette probabilite est calculee par interpolation lineaire
# Les cartes sont enregistrees au format jpeg dans le repertoire de travail

couleurAbsence <- "white"
couleurPresence <- "red"
couleurPresence <- colorRampPalette(c(couleurAbsence, couleurPresence))(6)

xSub <- tableP$Long
ySub <- tableP$Lat

#Pour chaque substrat, interpolation lineaire et affichage de la carte
for (sub in 3:length(names(tableP))) {
	zSub <- tableP[,sub]
	nomSub <- names(tableP)[sub]
	interpSub <- interp(xSub, ySub, zSub, x0, y0, linear=T, duplicate="mean")

	for (longitude in 1:dim(interpSub$z)[1]) {
		for (latitude in 1:dim(interpSub$z)[2]) {
			if (is.na(interpSub$z[longitude, latitude]) != T) {
				if (interpSub$z[longitude, latitude] > 1) {interpSub$z[longitude,latitude] <- 1}
				if (interpSub$z[longitude, latitude] < 0) {interpSub$z[longitude,latitude] <- 0}
			}
		}
	}
	eval(parse(text=paste0("interp.",nomSub,"=interpSub")))
	derTab <- eval(parse(text=paste0("interp.",nomSub,"=interpSub")))


# L'emplacement des legendes est regie par ces variables
xLegende <- max(x0)
yLegendeCouleur <- max(y0)
yLegendeTrajectoire <- min(y0)
yLegendeEchelle <- (9*min(y0)+max(y0))/10

xWindows <- (max(x0)-min(x0))*10000
yWindows <- (max(y0)-min(y0))*10000

# Affichage a proprement parler
windows(xWindows,yWindows)
par(xpd=TRUE, mai=c(1,1,1,2))
  
# Affichage de la carte de probabilite de presence du substrat sub
image(interpSub$x, interpSub$y, interpSub$z, 
	main=paste("Probabilité de presence du substrat",nomSub), col=couleurPresence,
	xlab="Longitude", ylab="Latitude", asp=1)

# Ajout du fond de carte
plot(lac, add=TRUE, border="black", lwd=3)

# Ajout de la legende
legend(xLegende, yLegendeCouleur, legend=seq(0,1,0.2), fill=couleurPresence, bg="white")
legend (xLegende, yLegendeTrajectoire, legend="trajectoire", lty=1, col="black")

# Ajout des points de calibration
points(tableCarto$Lat ~ tableCarto$Long, col="black", type="l")
map.scale(x=xLegende, y=yLegendeEchelle, relwidth = 0.15, metric = TRUE, ratio=F)
savePlot(paste0("probaPresence",nomSub), "pdf")
savePlot(paste0("probaPresence",nomSub), "jpeg")
dev.off()
}


#########################################################################
#########################################################################
# Etape 5 : affichage de la carte des substrats

#	# Etape 5.3 : Compilation des informations de presence

#########################################################################

# CarteSub est une matrice definissant le type de substrat attribue a chaque point
# la regle est :
# si 1 substrat > seuil -> c'est ce substrat
# si 2 substrats > seuil -> juxtaposition de ces substrats (melange) note substrat1_substrat2
# si plus de 2 substrats > seuil -> pas de substrat mais "melange"
# si aucun substrat > seuil -> IND
# seuil de 30% utilise par Francois Guiton (2012)
# seuil propose ici : 50% (a 30%, beaucoup de zones de melanges)

# carteSub est un tableau où les noms de colonne sont les longitudes et les noms de ligne les latitudes

seuil <- 0.5
carteSub <- derTab
carteSub$z[,] <- "IND"

for (longitude in 1:dim(carteSub$z)[1]) {
	for (latitude in 1:dim(carteSub$z)[2]) {
		for (sub in 3:(length(names(tableP))-1)) {
			nom <- names(tableP)[sub]
			titre <- get(paste0("interp.",nom))
			if (is.na(titre$z[longitude,latitude]) != T) {
				if ((titre$z[longitude,latitude]>=seuil) & (carteSub$z[longitude,latitude]=="IND")) { 
					carteSub$z[longitude,latitude] <- names(tableP)[sub]
				}
				if ((titre$z[longitude,latitude]>=seuil) & (carteSub$z[longitude,latitude]!="IND")) {
					if (str_detect(carteSub$z[longitude,latitude], names(tableP)[sub]) == F) {
						sub0 <- carteSub$z[longitude,latitude]
						carteSub$z[longitude,latitude] <- paste0(sub0,"_",nom)
					}
				}
			}
		}
	}
}


# On considere que lorsqu'il y a 3 substrats differents ou plus, il ne s'agit plus de plusieurs sediments distincts 
# mais d'un mélange
seuilMelange <- 2

for (longitude in 1:dim(carteSub$z)[1]) {
	for (latitude in 1:dim(carteSub$z)[2]) {
		if (length( unlist( strsplit(carteSub$z[longitude,latitude],"_") ) ) > seuilMelange) {
			carteSub$z[longitude,latitude] <- "melange"
		}
	}
}

#summary(carteSub$z)

# Lorsqu'un substrat a été estimé en dehors du lac, il est effacé
for (longitude in 1:dim(carteSub$z)[1]) {
	for (latitude in 1:dim(carteSub$z)[2]) {
		if (carteSub$z[longitude, latitude] != "IND") {
			if (inside.owin(carteSub$x[longitude], carteSub$y[latitude], lac) == F) {
				carteSub$z[longitude, latitude] <- "IND"
			}
		}
	}
}
#########################################################################
#########################################################################
# Etape 5 : affichage de la carte des substrats

#	# Etape 5.4 : Preparation des donnees pour l'affichage

#########################################################################

# on cree un tableau du meme format que tableSub mais au lieu d'un nom de substrat, on y implemente un chiffre
# on cree egalement une table de correspondance substrat_chiffre
# tableCorr est la table de correspondance

tableCorr <- matrix(0,0,1)
tableCorr <- rbind(tableCorr, carteSub$z[1,1])
for (longitude in 1:dim(carteSub$z)[1]) {
	for (latitude in 1:dim(carteSub$z)[2]) {
		if (is.na(carteSub$z[longitude,latitude]) != TRUE) {
			nom <- carteSub$z[longitude,latitude]
			ligne <- 1
			while (ligne<=nrow(tableCorr)) {
				if (tableCorr[ligne]==nom) {ligne <- 10^10}
				else {ligne <- ligne+1}
			}
			if (ligne != 10^10) {
				tableCorr <- rbind(tableCorr,nom)
			}
		}
	}
}

tableCorr2 <- as.character(unique(levels(tableSub$substratPred)))
tableCorr2 <- unique(c(tableCorr,tableCorr2))


# SubNum est une table semblable a carteSub mais avec des chiffres au lieu des noms
# la modification se fait en utilisant la table de correspondance tableCorr

subNum <- carteSub
subNum$z <- matrix(0,dim(carteSub$z)[1],dim(carteSub$z)[2])
for (longitude in 1:dim(carteSub$z)[1]) {
	for (latitude in 1:dim(carteSub$z)[2]) {
		nom <- carteSub$z[longitude,latitude]
		for (ligne in 1:length(tableCorr2)) {
			if (nom==tableCorr2[ligne]) {
				subNum$z[longitude,latitude] <- ligne
			}
		}
	}
}

# SubNum est une table semblable a carteSub mais avec des chiffres au lieu des noms
# la modification se fait en utilisant la table de correspondance tableCorr



#########################################################################
# Preparation des couleurs pour l'affichage
# Creation d'une matrice de correspondance substrat-couleur afin que toutes les
# cartes aient la meme legende

tableCoul <- matrix(0,30,3)
colnames(tableCoul) <- c("sub1","sub2","couleur")
tableCoul[1:9,1] <- c("IND","melange",names[3:9])
tableCoul[1:9,2] <- tableCoul[1:9,1]

a <- 3
b <- 4
for (ligne in 10:30) {
	tableCoul[ligne,1] <- paste0(tableCoul[a,1],"_",tableCoul[b,1])
	tableCoul[ligne,2] <- paste0(tableCoul[b,1],"_",tableCoul[a,1])
	if (b<9) { b <- b+1 }
	else { 
		a <- a + 1
		b <- a + 1 
	}
} 
tableCoul[,3] <- c("white","black",rainbow(7),
	colorRampPalette(c("dark green", "pink","dark blue"))(10), 
	colorRampPalette(c("yellow3", "aquamarine", "light green"))(11))


# Extraction des substrats présents sur le lac et de leurs couleurs associees
couleurSub <- matrix(0,length(tableCorr),2)
for (ligne in 1:length(tableCorr)) {
	for (ligneC in 1:nrow(tableCoul)) {
		if ((tableCorr[ligne]==tableCoul[ligneC,1]) | (tableCorr[ligne]==tableCoul[ligneC,2])) {
			couleurSub[ligne,1] <- tableCorr[ligne]
			couleurSub[ligne,2] <- tableCoul[ligneC,3]
		}
	}
}

#Extraction des substrats presents dans le jeu de calibration
couleurSub2 <- matrix("white",nrow(tableSub),1)
for (ligne in 1:nrow(tableSub)) {
	for (ligneC in 1:nrow(tableCoul)) {
		if ((tableSub$substratPred[ligne] == tableCoul[ligneC,1]) | (tableSub$substratPred[ligne] == tableCoul[ligneC,2])) {
			couleurSub2[ligne,1] <- tableCoul[ligneC,3]
		}
	}
}

# tableCorr2 donne les substrats indiques dans la legende
# On supprime donc les doublons de type substrat1_substrat2 et substrat2_substrat1
ligne1 <- 1
while (ligne1 <= length(tableCorr2)) {
	for (ligne2 in 10:nrow(tableCoul)) {
		if (tableCorr2[ligne1] == tableCoul[ligne2,1]) {
			ligne3 <- 1
			while (ligne3 <= length(tableCorr2)) {
				if (tableCorr2[ligne3] == tableCoul[ligne2,2]) {
					tableCorr2 <- tableCorr2[-ligne3]
					ligne1 <- 1
				}
				else {ligne3 <- ligne3 + 1}
			}
		}
	}
	ligne1 <- ligne1 + 1
}

# On cree une table de couleur pour la legende
coulLeg <- matrix(0,length(tableCorr2), 1) 
for (ligne in 1:nrow(coulLeg)) {
	for (ligneC in 1:nrow(tableCoul)) {
		if ((tableCorr2[ligne] == tableCoul[ligneC,1]) | (tableCorr2[ligne] == tableCoul[ligneC,2])) {
			coulLeg[ligne,1] <- tableCoul[ligneC,3]
		}
	}
}

#########################################################################
#########################################################################
# Etape 5 : affichage de la carte des substrats

#	# Etape 5.5 : Affichage et enregistrement de la carte des substrats

#########################################################################
# L'emplacement des legendes est regie par ces variables
xLegende <- max(x0)
yLegendeCouleur <- max(y0)
yLegendeTrajectoire <- min(y0)
yLegendeEchelle <- (9*min(y0)+max(y0))/10

xWindows <- (max(x0)-min(x0))*10000
yWindows <- (max(y0)-min(y0))*10000

# Affichage a proprement parler
windows(xWindows,yWindows)
par(xpd=TRUE, mai=c(1,1,1,3))

# Affichage du resultat de l'interpolation
a <- image(subNum$x,subNum$y,subNum$z,col=couleurSub[,2],
	main=paste("Carte des substrats de",nomLac2), 
	xlab="Longitude", ylab="Latitude", asp=1)

# On ajoute les legendes du graphique
legend(xLegende, yLegendeCouleur, 
	legend=tableCorr2, fill=coulLeg, xpd=NA, title="Substrats")
legend (xLegende, yLegendeTrajectoire, 
	legend=c("points de calibration", "trajectoire"),
	pch=c(4,NA), lty=c(NA,1), col="black")

# On ajoute l'echelle du graphique
map.scale(x=xLegende, y=yLegendeEchelle, relwidth = 0.1, metric = TRUE, ratio=F)

# Sur le graphique sont ajoutes les points de mesure (meme legende que la carte)
points(tableSub$Lat ~ tableSub$Long, lwd=2, pch=1, col=couleurSub2)

# Pour visualiser ces points qui se confondent car meme couleur, on trace (finement) la trajectoire du bateau
points(tableSub$Lat ~ tableSub$Long, lwd=0.5, pch=".", col="black")

# On ajoute les points de calibration
points(tableCalib$Lat ~ tableCalib$Long, col="black", pch=4)

# On ajoute le fond de carte du lac
plot(lac, add=TRUE, border="black", lwd=3)

# Enregistrement de la carte aux formats pdf et jpeg dans le repertoire de travail
savePlot(paste0("carteSubstrat",nomLac2), "pdf")
savePlot(paste0("carteSubstrat",nomLac2), "jpeg")

#########################################################################
# Affichage de la carte de repartition des substrats en 3 D (latitude, longitude et profondeur)

# On cree une matrice avec 4 colonnes : Lat, Long, Prof, Sub et codeSub
# a partir du subNum et carteSub

carte3D <- matrix(0,nrow(carteSub$z)*ncol(carteSub$z),6)
colnames(carte3D) <- c("Long", "Lat", "Prof", "Substrat", "CodeSubstrat", "Couleur")
for (i in 1:ncol(carteSub$z)) {
	ligneDebut <- 1 + (i-1)*nrow(carteSub$z)
	ligneFin <- i*nrow(carteSub$z)
	carte3D[ligneDebut:ligneFin, 1] <- carteSub$x
	carte3D[ligneDebut:ligneFin, 2] <- rep(carteSub$y[i],nrow(carteSub$z))
	carte3D[ligneDebut:ligneFin, 3] <- carteProf$z[,i]
	carte3D[ligneDebut:ligneFin, 4] <- carteSub$z[,i]
	carte3D[ligneDebut:ligneFin, 5] <- subNum$z[,i]
}
for (i in 1:nrow(carte3D)) {
	for (j in 1:nrow(tableCoul)) {
		if ((carte3D[i,4] == tableCoul[j,1]) | (carte3D[i,4] == tableCoul[j,2])) {
			carte3D[i,6] <- tableCoul[j,3]
		}
	}
}

substratsInterpoles <- carte3D[,c("Long","Lat","Prof","Substrat")]
substratsInterpoles <- na.omit(substratsInterpoles)
write.csv2(substratsInterpoles, file=paste0("substratsInterpole",nomLac,".csv"), row.names=F)
	
open3d()
par3d(windowRect = c(20, 30, 800, 800))
plot3d(as.numeric(carte3D[,1]), as.numeric(carte3D[,2]), as.numeric(carte3D[,3]),
	col=carte3D[,6], xlab="Longitude", ylab="Latitude", zlab="Profondeur",
	aspect=c(5,5,1))
legend3d("topright", legend=tableCorr2, fill=coulLeg,
	title="carte 3D des substrats")

# Enregistrement au format PNG mais perte de la possibilité de tourner le graphique
rgl.snapshot(paste0("carte3D",nomLac,".png"), fmt = "png", top = TRUE )

#########################################################################
# Affichage du substrat dans un graphique 3D (E1, E2 et profondeur)
# il est possible de faire tourner le graphique avec le clic gauche de la souris

open3d()
par3d(windowRect = c(20, 30, 800, 800))
plot3d(tableSub$Prof ~ tableSub$E1 * tableSub$E2, col=couleurSub2)
legend3d("topright", legend=unique(tableSub$substratPred), fill=unique(couleurSub2),
	title="repartition des substrats selon E1, E2 et profondeur")

# Enregistrement au format PNG mais perte de la possibilité de tourner le graphique
rgl.snapshot("graph3D.png", fmt = "png", top = TRUE )

#########################################################################
#########################################################################
# Etape 6 : aide a l'edition des carres RoxAnn

#	# Etape 6.1 : Un rectangle simple englobant tous les points d'un meme substrat

#########################################################################
# Creation d'une table de travail
# Seuls E1 et E2 sont utilises pour estimer le substrat

table1 <- tableSub[,c("E1","E2","substratPred")]
colnames(table1) <- c("E1", "E2", "Substrat")
table2 <- tableCalibRed[,c("E1","E2","Substrat")]
table <- rbind(table1, table2)

#########################################################################
# Creation des rectangles

substrat <- matrix(0, length(levels(table$Substrat)),1)
taille <- matrix(0, length(levels(table$Substrat)),4)
colnames(taille) <- c("E1min","E1max","E2min","E2max")
for (sub in 1:length(levels(table$Substrat))) {
	substrat[sub,1] <- levels(table$Substrat)[sub]
	tableExtrait <- subset(table, table$Substrat==substrat[sub,1])
	taille[sub,1] <- min(tableExtrait$E1)
	taille[sub,2] <- max(tableExtrait$E1)
	taille[sub,3] <- min(tableExtrait$E2)
	taille[sub,4] <- max(tableExtrait$E2)
}
substrat0 <- substrat
taille0 <- taille

#########################################################################
########################################################################## Etape 6 : Aide a l'edition des carres RoxAnn

#	# Etape 6.2 : Reperage des croisements entre carres

#########################################################################

for (i in 1:(nrow(taille)-1)) {
	for (j in (i+1):nrow(taille)) {
		for (colonnei in 1:4) {
			compteur <- 0
			if ((colonnei==1) | (colonnei==2)) {colonnej <- 1}
			if ((colonnei==3) | (colonnei==4)) {colonnej <- 3}
			if ((taille[i,colonnei]>=taille[j,colonnej]) & (taille[i,colonnei]<=taille[j,(colonnej+1)])) {
				newSub <- paste0(substrat[i],"_",substrat[j])
				E1min <- max(taille[i,1],taille[j,1])
				E1max <- min(taille[i,2],taille[j,2])
				E2min <- max(taille[i,3],taille[j,3])
				E2max <- min(taille[i,4],taille[j,4])
				newTaille <- c(E1min, E1max, E2min, E2max)
								
				
				for (ligne in 1:nrow(taille)) {
					if ((taille[ligne,1]==newTaille[1]) &
						(taille[ligne,2]==newTaille[2]) &
						(taille[ligne,3]==newTaille[3]) &
						(taille[ligne,4]==newTaille[4])) {
						 compteur <- compteur + 1
					}
				}
				if (compteur ==0) {
					substrat <- rbind(substrat,newSub)
					taille <- rbind(taille, newTaille)
				}
			}
		}
	}
}

#########################################################################
#########################################################################
# Etape 6 : aide a l'edition des carres RoxAnn

#	# Etape 6.3 : Nettoyage du jeu de donnes ainsi obtenu

#########################################################################
# enlever "IND" lorsqu'il y a d'autres substrats
# ex : sable_IND devient sable

for (ligne in 1:nrow(substrat)) {
	if (str_detect(substrat[ligne],"_")==T) {
		if (str_detect(substrat[ligne],"IND")==T) {
			substrat[ligne] <- gsub("IND_","",substrat[ligne])
			substrat[ligne] <- gsub("_IND","",substrat[ligne])
		}
	}
}

#########################################################################
# lorsqu'il y a deux fois le meme substrat dans un melange, on ne garde qu'une occurrence
# ex : sable_sable_gravier devient sable_gravier

# dans substrat0, on ne conserve que les substrats unitaires (sable, vase...)
a <- data.frame()
for (ligne in 1:nrow(substrat0)) {
	if (str_detect(substrat0[ligne],"_")==T) {
		a <-rbind(a,ligne)		
	}
}
substrat0 <- substrat0[-a[[1]],]

# Pour chaque substrat de "substrat", on regarde s'ils n'ont pas plusieurs fois le 
# meme substrat unitaire

for (nom in 1:length(substrat0)) {
	for (ligne in 1:nrow(substrat)) {
		if (str_detect(substrat[ligne],"_")==T) {
			a <- str_locate_all(substrat[ligne],substrat0[nom])
			if (dim(a[[1]])[1]>1) {
				substrat[ligne] <- sub(paste0(substrat0[nom],"_"),"",substrat[ligne])
			}
		}
	}
}

#########################################################################
# Lorsque le nombre de sediments possibles vaut 3 ou plus, 
# on considere qu'il s'agit d'un melange

seuilMelange <- 2
for (ligne in 1:nrow(substrat)) {
	if (str_detect(substrat[ligne],"_")==T) {
		a <- str_locate_all(substrat[ligne],"_")
		if (dim(a[[1]])[1]> (seuilMelange-1)) {
			substrat[ligne] <- "melange"
		}
	}
}

#########################################################################
# Supprimer les rectangles qui n'apparaissent pas car completement recouverts

# Auparavant, on limite le nombre de décimales a 2 chiffres apres la virgule
# (car RoxAnn a une precision de cet ordre)

taille <- round(taille,2)

compteur <- data.frame()

for (ligne in 1:(nrow(taille)-1)) {
	for (ligne2 in (ligne+1):nrow(taille)) {
		if ((taille[ligne,1] >= taille[ligne2,1]) &
			(taille[ligne,2] <= taille[ligne2,2]) &
			(taille[ligne,3] >= taille[ligne2,3]) &
			(taille[ligne,4] <= taille[ligne2,4])) {
			compteur <- rbind(compteur,ligne)
		}
	}
}


substrat2 <- as.character(data.frame())
taille2 <- data.frame()

if (is.null(compteur[1,1]) != T) {
	for (ligne in 1:nrow(taille)) {
		test <- 0
		for (c in 1:nrow(compteur)) {
			if (ligne==compteur[c,1]) { test <- 1 }
		}
		if (test==0) {
			taille2 <- rbind(taille2,taille[ligne,])
			substrat2 <- rbind(substrat2,substrat[ligne,])
		}
	}
}
if (is.null(compteur[1,1]) == T) {
	substrat2 <- as.character(substrat)
	taille2 <- taille
}

colnames(taille2) <- colnames(taille)

#########################################################################
#########################################################################
# Etape 6 : aide a l'edition des carres RoxAnn

#	# Etape 6.4 : Reorganisation de l'ordre des rectangles
#	# Et edition de la liste des coordonnees des carres RoxAnn

#########################################################################
# Ranger les lignes dans l'ordre suivant :
# (1) substrat unique
# (2) substrat double
# (3) melange


#on cree un tableau unique avec les noms et les valeurs des rectangles
carre <- cbind(substrat2,taille2)	

#on separe les lignes selon que leur substrat est unique double, ou un melange
subUnique <- matrix(0,0,dim(carre)[2])
subDouble <- matrix(0,0,dim(carre)[2])
subMelange <- matrix(0,0,dim(carre)[2])

for (ligne in 1:dim(carre)[1]) {
	if (str_detect(carre[ligne,1],"_")==T) {
		subDouble <- rbind(subDouble,carre[ligne,])
	}
	else if (carre[ligne,1]=="melange") {
		subMelange <- rbind(subMelange, carre[ligne,])
	}
	else {
		subUnique <- rbind(subUnique,carre[ligne,])
	}
}

subUnique <- subUnique[order(subUnique[,1],decreasing=F), ]
subDouble <- subDouble[order(subDouble[,1],decreasing=F), ]


liste <- rbind(subUnique, subDouble, subMelange)


#ajout d'une ligne en bas avec les mini et maxi
E1m <- min(liste[,2])
E1M <- max(liste[,3])
E2m <- min(liste[,4])
E2M <- max(liste[,5])
extremes <- c("extremes",E1m,E1M,E2m,E2M)
levels(liste[,1]) <- c(levels(liste[,1]),"extremes")
liste <- rbind(liste,extremes)

write.csv(liste,"coordRoxBox.csv")

#########################################################################
#########################################################################
# Etape 6 : aide a l'edition des carres RoxAnn

#	# Etape 6.5 : Affichage des rectangles finaux

#########################################################################


# Preparation des couleurs pour avoir des couleurs communes entre les points et les rectangles
 
coulBox1 <- unique(liste[,1])[-length(unique(liste[,1]))]
coulBox2 <- levels(table[,3])
coulBox <- c(as.character(coulBox1),as.character(coulBox2))
coulBox <- unique(coulBox)	
codeCouleur <- rainbow(length(coulBox))

tableCoulBox <- matrix(0,nrow(table),1)
for (ligne in 1:nrow(table)) {
	for (col in 1:length(coulBox)){
		if (table$Substrat[ligne]==coulBox[col]) {
			tableCoulBox[ligne] <- col
		}
	}
}
table <- cbind(table, tableCoulBox)

tableCoulBox2 <- matrix(0,nrow(liste),1)
for (ligne in 1:nrow(liste)) {
	for (col in 1:length(coulBox)){
		if (liste[ligne,1]==coulBox[col]) {
			tableCoulBox2[ligne] <- col
		}
	}
}
liste <- cbind(liste, tableCoulBox2)

# Affichage a proprement parler

windows(150,100)
par(xpd=TRUE, mai=c(1,1,1,3))
plot.new()
palette(rainbow(length(coulBox))) 

plot(E1 ~ E2, data=table, 
	main="Visualisation des carres RoxAnn et des points de calibration")

for (sub in 1:(nrow(liste)-1)) {
	x <- c(as.numeric(liste[sub,4]), as.numeric(liste[sub,5]), as.numeric(liste[sub,5]), as.numeric(liste[sub,4]))
	y <- c(as.numeric(liste[sub,2]), as.numeric(liste[sub,2]), as.numeric(liste[sub,3]), as.numeric(liste[sub,3]))
	points(x,y, col=as.numeric(tableCoulBox2[sub]))
	polygon(x,y, col=as.numeric(tableCoulBox2[sub]))
}

points(E1 ~ E2, data=table, col=tableCoulBox)
legend(1.01*max(table$E2),max(table$E1),legend=coulBox, fill=1:length(coulBox), xpd=NA)

savePlot(paste0("graphiqueRoxBox",nomLac), "pdf")
savePlot(paste0("graphiqueRoxBox",nomLac), "jpeg")

#########################################################################
#########################################################################
# Etape 7 : Creation et affichage d'un fichier recapitulatif (format texte)

#########################################################################

# Calcul de la vitesse le long de la trajectoire (en km)
distance <- matrix(0,nrow(tableCartoRed)-1,1)
for (ligne in 1:nrow(distance)) {
	p1 <- c(tableCartoRed$Long[ligne],tableCartoRed$Lat[ligne])
	p2 <- c(tableCartoRed$Long[ligne+1],tableCartoRed$Lat[ligne+1])
	distance[ligne,1] <- distCosine(p1,p2)/1000
}

# Calcul du temps entre chaque point de mesure (en heure)
temps <- matrix(0,nrow(tableCartoRed)-1,1)
for (ligne in 1:nrow(distance)) {
	t1 <- strptime(tableCarto$Heure[ligne], format="%H:%M:%S")
	t2 <- strptime(tableCarto$Heure[ligne+1], format="%H:%M:%S")
	temps[ligne,1] <- time_length(interval(t1, t2), "hours")
}

# Calcul de la vitesse en chaque point (en km/h)
vitesse <- distance
for (ligne in 1:nrow(vitesse)) {
	vitesse[ligne, 1] <- distance[ligne, 1]/temps[ligne, 1]
}

# Vitesse exprimee en km/h
vitesseMoyenne <- round(mean(vitesse),2)
vitesseMax <- round(max(vitesse),2)
vitesseEcartType <- round(sd(vitesse),2)

# On rappelle l'estimation du taux d'erreur du modele
erreurModele <- 1-round(calibBC,2)

# Informations concernant les substrats
# Nombre de substrat
nombreSubstrats <- length(tableCorr2) 

# tableau recapitulatif des substrats
tableauSubstrat <- matrix(0, nombreSubstrats, 3)
rownames(tableauSubstrat) <- tableCorr2
colnames(tableauSubstrat) <- c("nb cellules", "% de surface", "profondeur moyenne")

for (longitude in 1:nrow(carteSub$z)) {
	for (latitude in 1:ncol(carteSub$z)) {
		for (ligne in 1:nrow(tableauSubstrat)) {
			if (carteSub$z[longitude, latitude] == rownames(tableauSubstrat)[ligne]) {
				tableauSubstrat[ligne, 1] <- tableauSubstrat[ligne, 1]+1
			}
		}
	}
}
nSub <- sum(tableauSubstrat[-1,])
for (ligne in 2:nrow(tableauSubstrat)) {
	tableauSubstrat[ligne, 2] <- round((tableauSubstrat[ligne,1]/nSub)*100,2)
}

# Estimation de la profondeur moyenne de chaque substrat
# La moyenne n'est calculee qu'a partir des donnees le long de la trajectoire

for (ligne in 2:nrow(tableauSubstrat)) {
	profondeur <- subset(tableSub, tableSub$substratPred==rownames(tableauSubstrat)[ligne])
	tableauSubstrat[ligne, 3] <- mean(profondeur$Prof)
}

tableauSubstrat <- tableauSubstrat[,-1]
tableauSubstrat <- subset(tableauSubstrat, rownames(tableauSubstrat) != "IND")
tableauSubstrat <- tableauSubstrat[order(tableauSubstrat[,1],decreasing=TRUE),] 

sink("recapitulatif.txt")
cat("\n","**Donnees concernant la vitesse d'acquisition**", "\n")
cat(paste("La vitesse moyenne est de", vitesseMoyenne, "km/h"), "\n")
cat(paste("La vitesse maximale atteinte est de", vitesseMax, "km/h"), "\n")
cat(paste("Et l'ecart type vaut", vitesseEcartType),"\n")
cat("\n")
cat("**Donnee concernant l'efficacite du modele de prediction du substrat le long de la trajectoire**","\n")
cat(paste("Le taux d'erreur du modele est estime a", erreurModele, ", soit", erreurModele*100, "% de substrats non correctement classes"),"\n")
cat("\n")
cat("**Donnees concernant les substrats estimes**","\n")
cat("Attention, les profondeurs moyennes sont calculees uniquement a partir des mesures le long de la trajectoire","\n")
cat(paste(length(tableCorr2)-1, "substrats ont ete definis comme presents sur la zone d'etude"),"\n")
cat("\n")
print(tableauSubstrat)
sink() 

#########################################################################
#########################################################################

#Fin du script

##################################################################################################################################################