coopR: Fonctions R pour la Coop

#' @title Chargement des fichiers COOP depuis fichier xls stocke dans le T:R&D (uniquement pour utilisateur ONF !)
#'
#' @description Fonction (normalement temporaire) pour lire les donnees Coop a partir des fichier xls stocke dans T:R&D
#'
#' Ce programme est temporaire puisqu a terme, tout le monde aura acces a la BDD. Il est, de plus, sensible au changement de nom de dossier. Donc attention!
#' 
#' @details Cette fonction est un deuxieme jet. Les fonctions d'importation depuis Excel dependent de la
#' \code{method} choisie : \code{Java} pour la fonction \code{\link[XLConnect]{readWorksheet}} du package \code{XLConnect}
#'
#' @param initiales : vecteur contenant les initiales du ou des dispositifs que l'on veut charger (par defaut, tous !)
#' @param format : plusieurs possibilite de format 
#' - "xls" : list de data.frame correspondant exactement au fichier xls de la Coop ou au chargement de la Coop par la fonction load_bdd_to_xls_ch 
#' - "R" : format completement different, identique au chargement par la fonction load_bdd_to_R_ch et qui sert pour la suite des calculs
#' @param cl_circ data.frame a 2 colonnes : circonference puis effectif
#' @param SurfPla surface en m2
#' @param N nombre de tige par ha, 100 par defaut.
#' @return le diametre dominant
#' Une liste de data.frame:
#' - pour (format = "xls") : Placette, Mesure, Arbre, Ss_unite_surface, Inventaire_liste, Intervention, Historique_traitement, Numero_arbre identiques aux tables issues de la fonction load_bdd_to_xls_ch
#' - pour (format = "R") : unite, mesure, mesure2, arbre, arbre_echantillon, inventaire_somme, surface_somme identiques aux tables issues de la fonction
#' @note Fonction reprise du package interne 'onfR'
#' @note Merci de faire remonter les critiques et ameliorations de la fonction
#' @author Quentin Girard         
#' @references Protocole Coop chene...
#'
#' @seealso  dico.ch, load_bdd_to_R_ch, load_bdd_to_xls_ch, load_bdd_to_R_ch
#' 
#' @keywords function
#'
#' @include coopR-package.R
#' @family coopR
#' @export
#' @importFrom ggplot2 ggplot aes geom_line geom_point theme_bw 
#' @importFrom XLConnect readWorksheet getSheets loadWorkbook
#' @examples
#' ## Utilisons l'exemple avec le format = "R"
#'  require(XLConnect)
#'  require(doBy)
#'  require(plyr)
#'  ex <- load_onf_ch(initiales = c("MBM", "TR", "MR"), format = "R")
#'  str(ex)
#' ## graphe de verification de mesure2
#'  g <- ggplot(ex$mesure2) +
#'   aes(x = no_mesure, y = id_unite, col = suivi) +
#'   geom_line() +
#'   geom_point(aes(shape = as.factor(av_ap_id))) + scale_shape_manual(values = c("0" = 1, "1" = 3, "2" = 4))  +
#'   theme_bw()
#'  print(g)
#' 

################################################################################
# Cette fonction (normalement temporaire), charge les donnees Coop 
# a partir des fichier xls stocke dans T:R&D
#
# Auteur : Q. Girard, ONF RD, Pole de Boigny
# Date : octobre 2015
################################################################################


load_onf_ch <- function(initiales = NULL, format = "xls") {

require(XLConnect)
require(doBy)
require(plyr)
# verification des initiales 
  dispo_existant <- c("AGR","BV","CX","DACHP","DACHS","GB","GR","IZ","LFCHP","LFCHS","MBM","MR","OR","PA","RV","SA","SC","SP","TR","TRA","VE")
  if (is.null(initiales)) {
    initiales <- dispo_existant
  } else if (length(setdiff(initiales, dispo_existant)) != 0) {
    stop(paste("le(s) dispositif(s)", list(setdiff(initiales, dispo_existant)), "n'existe(nt) pas, verifiez l'orthographe"))
  }  
  
# list de sortie                                     
  if (format == "xls") {
    list.tp <- list("Placette"              = data.frame(), "Historique_traitement" = data.frame(),
                    "Arbre"                 = data.frame(), "Mesure"                = data.frame(),
                    "Numero_arbre"          = data.frame(), "Ss_unite_surface"      = data.frame(),
                    "Inventaire_liste"      = data.frame(), "Intervention"          = data.frame()) 
  } else if (format == "R") {                        
    list.tp <- list("Placette"              = data.frame(), "Arbre"                 = data.frame(), 
                    "Mesure"                = data.frame(), "Ss_unite_surface"      = data.frame(),
                    "Inventaire_liste"      = data.frame()) 
  } else {
    stop("format doit etre egal a 'xls' (chargement complet du format xls) ou 'bdd' (chargement des tables pour les calculs)")
  }
  
# itineraire vers xls                                                                               
  iti <- "//arbre.foret.onf.fr/Reseaux/R&D/03-Groupes de travail/06- Coop de donnees/Groupe coop chenes/Coop"
# Verification de l'existance d'un seul fichier a chaque fois (fait au debut pour ne pas tourner inutilement
  for (val in initiales) {                                
    if (length(dir(path = paste(iti, val, "/z_data_compl/", sep=""), full.names = T, pattern = ".xls")) >= 2) {
      stop("Il y a deux fichiers excel dans le dossier Coop", val, "/z_data_compl. Transferez le plus ancien dans le dossier Coop", val, "/z_data_compl/ancien_perime")
    }
  }
  
# pas de doubles fichiers excel ==> c'est parti  
# listes des intitules de colonnes qui pose probleme
  var.jet1 <- c()                       
  var.jet2 <- c() 
  ini.vide <- c()                                 
# Pour tous les sites etudies...
  for (val in initiales) {    
#... recherche de l'adresse...
    fichier <- dir(path = paste(iti, val, "/z_data_compl/", sep=""), full.names = T, pattern = ".xls")                            
#... verification qu'il y a le fichier...
    if (length(fichier) == 0) {
      ini.vide <- c(ini.vide, val)
    } else {                             
# ... ouverture du lien vers le fichier
      wb <- loadWorkbook(fichier)      
# pour toutes les feuilles de donnees du fichier...
      for(feuil in names(list.tp)) {             
# ... chargement de la feuille  
        data.tp <- readWorksheet(wb, feuil)   
        # verification des intitules (que ce qui sont repertories dans dico.ch)  =====> TEMPORAIRE ! Besoin d'une feuille dans dico_var a ce sujet
        var.jet1 <- c(var.jet1, setdiff(names(data.tp), dico.ch$dico.var$variable))
        data.tp <- data.tp[, intersect(names(data.tp), dico.ch$dico.var$variable)]
# pas encore de site integre...        
        if (nrow(list.tp[[feuil]]) == 0) {   
          list.tp[[feuil]] <- data.tp                                            
        } else {
# ... ou bien si des colonnes ne sont pas communes aux differents sites                                        
          var.jet2 <- c(var.jet2, 
                        setdiff(names(data.tp), names(list.tp[[feuil]])),   # colonne de la nouvelle table absente des precedentes
                        setdiff(names(list.tp[[feuil]]), names(data.tp)))   # colonne des tables precedentes absentes de la nouvelle
          list.tp[[feuil]] <- rbind.fill(list.tp[[feuil]], data.tp)
          #list.tp[[feuil]] <- rbind(list.tp[[feuil]][,intersect(names(data.tp), names(list.tp[[feuil]]))], 
          #                          data.tp[,intersect(names(data.tp), names(list.tp[[feuil]]))])               
        }                  
      }
    }
  }                                   
# Pour information
  if (length(ini.vide) != 0) {                                                                                           
    warning(paste("Il n'y pas de fichier excel dans le dossier z_data_compl pour les dispositifs", list(ini.vide)))
  }            
  if (length(var.jet1) != 0) {
    warning(paste("Les variables", paste(list(unique(var.jet1))), "ne sont pas repertoriees dans dico.ch : elles ne sont pas conservees."))
  }                                 
  if (length(var.jet2) != 0) {
    warning(paste("Les variables ", paste(list(unique(var.jet2))), "ne sont pas presentes dans toutes les tables"))# : elles ne sont pas conservee."))
  }
# quelques complement
  list.tp[["Placette"]]$an_naissance <- with(list.tp[["Placette"]], ifelse(is.na(an_cernes030), 
                                                                           calc_veg(annee_plant, mois_plant_id) - age_plant,
                                                                           calc_veg(an_cernes030, mois_cernes030_id) - nb_cernes030))   
  list.tp[["Arbre"]]$c130 <- with(list.tp[["Arbre"]], ifelse(is.na(circonference130), diam_130 * pi, circonference130))
# Attention, correction temporaire !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                                             
  #list.tp[["Ss_unite_surface"]] <- subset(list.tp[["Ss_unite_surface"]], (id_placette == 20001 & av_ap_id == 2 & no_mesure == 2))  
  #list.tp[["Inventaire_liste"]] <- subset(list.tp[["Inventaire_liste"]], !(id_placette == 20001 & av_ap_id == 2 & no_mesure == 2))
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
# Si on veut un nouveau format, celui utilise pour les calculs par la suite
  if (format == "R") {
# table unite
    unite_pl <- data.frame(list.tp[["Placette"]][c("lib_disp", "id_placette","co_traitement",
                                                   "essence_principale1_id","essence_principale2_id","surface","an_naissance")], 
                           pl_sp = "pl") 
      unite_pl$id_unite <- unite_pl$id_placette
      #names(unite_pl)[which(names(unite_pl) == "id_placette")] <- "id_unite"    
    unite_sp <- data.frame(subset(list.tp[["Placette"]], existence_ss_unite_id == 1, 
                                  select = c("lib_disp", "id_placette","co_traitement",
                                             "essence_principale1_id","essence_principale2_id","an_naissance")),
                           pl_sp = "sp")
      unite_sp <- merge(unite_sp, subset(list.tp[["Ss_unite_surface"]], no_mesure ==-2 & substr(no_ss_unite,1,2) == "sp", 
                        select = c("id_placette","no_ss_unite","surface"), all.x = T))
      unite_sp$id_unite <- paste(unite_sp$id_placette, unite_sp$no_ss_unite, sep = "_")
    unite <- rbind(unite_pl, unite_sp[,intersect(names(unite_pl), names(unite_sp))])  

# table mesure
    mesure_pl <- subset(list.tp[["Mesure"]], type_suivi_id != 4, select = c("id_placette", "no_mesure", "an_vegetation", "type_suivi_id"))
      names(mesure_pl)[names(mesure_pl) == "id_placette"] <- "id_unite"      
    mesure_sp <- subset(list.tp[["Mesure"]], type_suivi_id %in% c(4:7), select = c("id_placette", "no_mesure", "an_vegetation", "type_suivi_id"))
      mesure_sp <- merge(mesure_sp, unite_sp[, c("id_placette", "id_unite")])
    mesure <- rbind(mesure_pl, mesure_sp[,intersect(names(mesure_pl), names(mesure_sp))])
    mesure <- merge(mesure, unite[,c("id_unite", "an_naissance")])
    mesure$age <- mesure$an_vegetation - mesure$an_naissance
    mesure <- mesure[, -which(names(mesure) == "an_naissance")]
    
# table arbre
    list.tp[["Arbre"]]$pl_sp <- with(list.tp[["Arbre"]], ifelse(is.na(id_mesure_echantillon_sp), "pl", "sp"))  
    list.tp[["Arbre"]] <- merge(list.tp[["Arbre"]], unite[, c("id_placette", "pl_sp", "id_unite")])
    names(list.tp[["Arbre"]])[which(names(list.tp[["Arbre"]]) == "no_arb_ech")] <- "no_arbre"
    arbre <- subset(list.tp[["Arbre"]], is.na(id_mesure_echantillon_bd), 
                    select = c("id_unite","no_mesure","no_arbre","essence_id","co_ech","etat_san_id","statut_av_id","c130","eclairci_realise_id"))
    # suppression des morts?

# table arbre_echantillon                
    arbre_echantillon <- merge(list.tp[["Arbre"]], mesure[,c("id_unite", "no_mesure", "type_suivi_id")], by = c("id_unite", "no_mesure"))
    arbre_echantillon <- subset(arbre_echantillon, co_ech %in% c(1, 2) | (co_ech== 3 & type_suivi_id %in% c(3, 7)), 
                                select = c("id_unite", "no_mesure", "no_arbre", "essence_id", "co_ech", "etat_san_id", "statut_av_id", "statut_ap_id",
                                           "c130", "hauteur_ech", "ht_pb_s", "ht_pb_g", "ht_pb_i", "ht_pb_v"))
                                           
# table inventaire_somme
    list.tp[["Inventaire_liste"]]$pl_sp <- with(list.tp[["Inventaire_liste"]], ifelse(substr(lib_bandelette,1,2) == "sp", "sp", "pl"))
    list.tp[["Inventaire_liste"]] <- merge(list.tp[["Inventaire_liste"]], unite[,c("id_placette", "pl_sp", "id_unite")])
    inventaire_somme <- summaryBy(effectif ~ id_unite + no_mesure + av_ap_id + no_inventaire + circ_diam_id + essence_id + borne_inf + borne_sup,
                                  list.tp[["Inventaire_liste"]], FUN = sum, na.rm = T, keep.names = T)                   
    inventaire_somme$c130 <- with(inventaire_somme, (borne_inf + borne_sup)/2 * ifelse(circ_diam_id == 1, pi, 1))
    inventaire_somme <- subset(inventaire_somme, select=-c(circ_diam_id, borne_inf, borne_sup))
    
# table surface_somme          
    list.tp[["Ss_unite_surface"]]$pl_sp <- with(list.tp[["Ss_unite_surface"]], ifelse(substr(no_ss_unite,1,2) == "sp", "sp", "pl"))
    list.tp[["Ss_unite_surface"]] <- merge(list.tp[["Ss_unite_surface"]], unite[,c("id_placette", "pl_sp", "id_unite")])
    surface_somme <- summaryBy(surface ~ id_unite + no_mesure + av_ap_id + no_inventaire,
                               subset(list.tp[["Ss_unite_surface"]], no_mesure != -2), FUN = sum, na.rm = T, keep.names = T)
    names(surface_somme)[which(names(surface_somme) == "surface")] <- "surf_inv"
       
# table mesure2 : table avec, pour chaque placette et campagne, une ligne par inventaire ou "etat" (avant-apres intervention) du peuplement. 
  # Le principe est d'avoir une ligne par inventaire (s'il y a plusieurs inventaires lors d'une meme campagne) et par stade du peuplement (avant et apres intervention pour une meme mesure, ou bien une ligne s'il n'y en a pas eu):
      # Pour le calcul a partir des bandelettes, c'est facile (on prend les indices de la table surface_somme)
    m2_bd <- merge(surface_somme, mesure)
      m2_bd$suivi    <- "classe"
      m2_bd <- m2_bd[,c("id_unite", "no_mesure", "an_vegetation", "av_ap_id", "no_inventaire", "surf_inv", "suivi", "age")]   
          
      # Pour les transitions bandelette- placette,
    m2_bd_pl <- subset(merge(mesure, unite[,c("id_unite", "pl_sp", "surface")]), type_suivi_id %in% c(3,7) & pl_sp == "pl")
      m2_bd_pl <- merge(m2_bd_pl, summaryBy(no_inventaire ~ id_unite + no_mesure + av_ap_id, subset(m2_bd, av_ap_id != 2), FUN = max, keep.names = T), 
                    all.x = T, suffixes = c("", ".x"))
        # attention, avec cette methode, il ne faut pas d'inventaire apres intervention par bandelette
      m2_bd_pl$no_inventaire <- ifelse(m2_bd_pl$av_ap_id %in% c(0,2), m2_bd_pl$no_inventaire + 1, 1)
      m2_bd_pl$av_ap_id      <- ifelse(m2_bd_pl$av_ap_id == 0, 0, 2)
      m2_bd_pl$suivi    <- "liste"               
      m2_bd_pl$surf_inv <- m2_bd_pl$surface           
      m2_bd_pl <- m2_bd_pl[,c("id_unite", "no_mesure", "an_vegetation", "av_ap_id", "no_inventaire", "surf_inv", "suivi", "age")] 
      
      # pour les suivi uniquement sur liste (placette sans transition et sous-placette) : intervention si au moins une tiges est coupee
    m2_pl_sp <- subset(merge(mesure, unite[,c("id_unite", "pl_sp", "surface")]), pl_sp == "sp" | (pl_sp == "pl" & type_suivi_id %in% c(1,5)))
      m2_pl_sp <- merge(m2_pl_sp, unique(subset(arbre, eclairci_realise_id %in% c(1,2), c("id_unite", "no_mesure", "eclairci_realise_id"))), all.x=T) 
      m2_pl_sp <- rbind(data.frame(subset(m2_pl_sp, is.na(eclairci_realise_id)), av_ap_id = 0, no_inventaire = 1),
                    data.frame(subset(m2_pl_sp, eclairci_realise_id == 1), av_ap_id = 1, no_inventaire = 1),
                    data.frame(subset(m2_pl_sp, eclairci_realise_id == 1), av_ap_id = 2, no_inventaire = 1))           
      m2_pl_sp$suivi    <- "liste"               
      m2_pl_sp$surf_inv <- m2_pl_sp$surface           
      m2_pl_sp <- m2_pl_sp[,c("id_unite", "no_mesure", "an_vegetation", "av_ap_id", "no_inventaire", "surf_inv", "suivi", "age")] 

      # fusion et attribution d'un numero d'inventaire depuis le debut      
    mesure2 <- rbind(m2_bd, m2_bd_pl, m2_pl_sp)
    mesure2 <- mesure2[with(mesure2, order(id_unite, no_mesure, av_ap_id, no_inventaire)),]
    mesure2[1, "ordre_inv"] <- 1
    for (val in 2:nrow(mesure2)) {
      mesure2[val, "ordre_inv"]  <- ifelse(mesure2[val, "id_unite"]       == mesure2[val - 1, "id_unite"], mesure2[val - 1, "ordre_inv"] + 1, 1)   
    }                    
# on refait list.tp
    list.tp <- list(unite = subset(unite, select = -id_placette), mesure = mesure, arbre = arbre, arbre_echantillon = arbre_echantillon, 
                    inventaire_somme = inventaire_somme, surface_somme = surface_somme, mesure2 = mesure2)      
  }  
  return(list.tp)
}