Mercurial > repos > mnhn65mo > butterfly_analysis
changeset 2:c29354d16967 draft
Uploaded
| author | mnhn65mo |
|---|---|
| date | Mon, 13 Aug 2018 04:55:07 -0400 |
| parents | c5a10acd34e3 |
| children | 9daba7983d38 |
| files | butterfly_crossplot.R butterfly_crossplot.xml clim_data.R code_couleurs.csv raster_getdata.xml stat_bag.r |
| diffstat | 6 files changed, 759 insertions(+), 344 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/butterfly_crossplot.R Mon Aug 13 04:55:07 2018 -0400 @@ -0,0 +1,556 @@ +################################################################## +#### Script generique pour realiser les figures en croix ###### +#### a partir des donnees brut ###### +################################################################## + +### Version V1.2 _ 2018-07-31 + +library(ggplot2) +library(RColorBrewer) + +args <- commandArgs(trailingOnly = TRUE) + +### importation code +sourcefunctions<-args[1] +source(sourcefunctions) + +## fonction d'importation des fichier des donnes +### fonction d'importation, de concatenation des fichiers +### verification des nom de colonnes +### verification des doublon de ligne +read.data <- function(file=NULL,decimalSigne=".") { +# cat("1) IMPORTATION \n--------------\n") +# cat("<--",file,"\n") + data <- read.table(file,sep="\t",stringsAsFactors=FALSE,header=TRUE,dec=decimalSigne) + ## verification qu'il y a plusieur colonnes et essaye different separateur + if(ncol(data)==1) { + data <- read.table(file,sep=";",stringsAsFactors=FALSE,header=TRUE,dec=decimalSigne) + if(ncol(data)==1) { + data <- read.table(file,sep=",",stringsAsFactors=FALSE,header=TRUE,dec=decimalSigne) + if(ncol(data)==1) { + data <- read.table(file,sep=" ",stringsAsFactors=FALSE,header=TRUE,dec=decimalSigne) + if(ncol(data)==1) { + stop("!!!! L'importation a echoue\n les seperatateurs de colonne utilise ne sont pas parmi ([tabulation], ';' ',' [espace])\n -> veuillez verifier votre fichier de donnees\n") + } + } + } + } + return(data) +} + + + + + +filtre1niveau <- function(func, + nom_fichier = filename, + dec=".", + nom_fichierCouleur= color_filename, + col_abscisse = "AB_MOYENNE", + figure_abscisse = "Abondance", + col_ordonnee = "DIVERSITE_MOYENNE", + figure_ordonnee = "Diversite", + nomGenerique="GLOBAL", + vec_figure_titre = c("Les Papillons"), + colourProtocole = TRUE, + nomProtocole = "Papillons", + vec_col_filtre = vec_col_filtre_usr, + col_sousGroup = NULL,# + val_filtre = NULL,# + figure_nom_filtre = NULL,# + bagplot = TRUE, + bagProp=c(.05,.5,.95), + seuilSegment=30, + segmentSousSeuil=TRUE, + forcageMajusculeFiltre=TRUE, + forcageMajusculeSousGroupe=TRUE){ + + dCouleur <- read.data(file=nom_fichierCouleur) + d <- read.data(file=nom_fichier,decimalSigne=dec) + if(colourProtocole & !is.null(nomProtocole)) colourProtocole_p <- as.character(dCouleur[dCouleur[,2]==nomProtocole,3]) else colourProtocole_p <- NULL + + for(f in 1:length(vec_col_filtre)) { + if(length(vec_figure_titre)==1){ + figure_titre_f <- vec_figure_titre + }else{ + figure_titre_f <- vec_figure_titre[f] + } + col_filtre_f <- vec_col_filtre[f] + cat(col_sousGroup) #Just to check + if(func=="ggfiltre1niveau"){ + cat("ggfiltre1niveau") + ggfiltre1niveau(d, + col_abscisse, + figure_abscisse, + col_ordonnee, + figure_ordonnee, + figure_titre = figure_titre_f, + col_filtre = col_filtre_f, + nomGenerique, + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= subset(dCouleur,Filtre==col_filtre_f), + colourProtocole = colourProtocole_p, + nomProtocole, + bagplot, + bagProp=c(.05,.5,.95), + seuilSegment, + segmentSousSeuil, + forcageMajusculeFiltre) + }else if(func=="gglocal"){ + cat("gglocal") + gglocal(d, + col_abscisse, + figure_abscisse, + col_ordonnee, + figure_ordonnee, + figure_titre = figure_titre_f, + col_filtre = col_filtre_f, + nomGenerique = nomGenerique, + col_sousGroup = col_sousGroup, + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= subset(dCouleur,Filtre==col_filtre_f), + colourProtocole = colourProtocole_p, + nomProtocole, + couleurLocal="#f609c1", + bagplot, + bagProp, + seuilSegment, + segmentSousSeuil, + forcageMajusculeFiltre, + forcageMajusculeSousGroupe) + }else{ + cat("ggCompareLevel") + ggCompareLevel(d, + col_abscisse, + figure_abscisse, + col_ordonnee, + figure_ordonnee, + figure_titre = figure_titre_f, + col_filtre = col_filtre_f, + nomGenerique = nomGenerique, + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= subset(dCouleur,Filtre==col_filtre_f), + colourProtocole = colourProtocole_p, + nomProtocole, + bagplot, + bagProp, + seuilSegment, + segmentSousSeuil, + forcageMajusculeFiltre) + } + } +} + +ggfiltre1niveau <- function(d, + col_abscisse = "AB_MOYENNE", + figure_abscisse = "Abondance", + col_ordonnee = "DIVERSITE_MOYENNE", + figure_ordonnee = "Diversite", + figure_titre = "Referentiel papillon", + col_filtre = "nom_reseau", + nomGenerique = "Global", + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= NULL, + colourProtocole = NULL, + nomProtocole = NULL, + bagplot = TRUE, + bagProp=c(.05,.5,.95), + seuilSegment=30, + segmentSousSeuil=TRUE, + forcageMajusculeFiltre=TRUE, + result_dir="resultats/") { + + d$groupe <- as.character(d[,col_filtre]) + d$abscisse <- d[,col_abscisse] + d$ordonnee <- d[,col_ordonnee] + d$groupe <-gsub("/","_",d$groupe) + d$groupe <-gsub("!","",d$groupe) + + if(forcageMajusculeFiltre){ + d$groupe <- toupper(d$groupe)} + + d <- subset(d,!(is.na(groupe)) & !(is.na(abscisse)) & !(is.na(ordonnee)) & groupe != "") + + if(is.null(val_filtre)){ + lesModalites <- unique(d$groupe) + }else{ + lesModalites <- val_filtre + } + +# repResult <- dir(result_dir) +# current_dir<-getwd() +# dir.create(file.path(current_dir,result_dir)) +# +# if(!(col_filtre %in% repResult)){ +# dir.create(file.path(".",paste(result_dir,col_filtre,sep="")))} +# +# nomRep1 <- paste(result_dir,col_filtre,"/",sep="") + + d.autre <- d + d.autre$groupe <- nomGenerique + + for(m in lesModalites) { + d.reseau <- subset(d,groupe==m) + d.reseau$groupe <- m + ggTable <- rbind(d.autre,d.reseau) + + seuilResum <- nrow(d.reseau) >= seuilSegment + + ggTableResum <- aggregate(cbind(ordonnee, abscisse) ~ groupe, data = ggTable,quantile, c(.25,.5,.75)) + ggTableResum <- data.frame(ggTableResum[,1],ggTableResum[,2][,1:3],ggTableResum[,3][,1:3]) + colnames(ggTableResum) <- c("groupe","ordonnee.inf","ordonnee.med","ordonnee.sup","abscisse.inf","abscisse.med","abscisse.sup") + + if(ggTableResum$groupe[2]==nomGenerique){ + ggTableResum <- ggTableResum[c(2,1),]} + + if(!(is.null(tab_figure_couleur))) { + if(m %in% tab_figure_couleur$Modalite) { + figure_couleur <- setNames(c(as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == nomGenerique]), + as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == m])), + c(nomGenerique,m)) + }else{ + figure_couleur <- setNames(c(as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == nomGenerique]), + as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == ""])), + c(nomGenerique,m)) + } + } + +# repResult <- dir(nomRep1) +# if(!(m %in% repResult)){ +# dir.create(paste(nomRep1,m,sep=""))} +# nomRep <- paste(nomRep1,m,"/",sep="") +# +# +# if(!is.null(nomProtocole)){ +# repResult <- dir(nomRep) +# if(!(nomProtocole %in% repResult)){ +# dir.create(paste(nomRep,nomProtocole,sep=""))} +# nomRep <- paste(nomRep,nomProtocole,"/",sep="") +# } + + + gg <- ggplot(ggTable,aes(x=abscisse,y=ordonnee,colour=groupe,fill=groupe)) + if(bagplot){ + gg <- gg + stat_bag(data=d.autre,prop=bagProp[1],colour=NA,alpha=.7) + stat_bag(data=d.autre,prop=bagProp[2],colour=NA,alpha=.4) + stat_bag(data=d.autre,prop=bagProp[3],colour=NA,alpha=.2) } + else { + gg <- gg + geom_point(alpha=.2) + } + gg <- gg + geom_hline(data=subset(ggTableResum,groupe== nomGenerique),aes(yintercept = ordonnee.med,colour=groupe),size=.5,linetype="dashed") + geom_vline(data=subset(ggTableResum,groupe==nomGenerique),aes(xintercept = abscisse.med,colour=groupe),size=.5,linetype="dashed") + if(segmentSousSeuil) { + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.8,size=2.5) + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.8,size=2.5) + if(!(seuilResum)) { + gg <- gg + geom_segment(data=subset(ggTableResum,groupe!=nomGenerique),aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.5,size = 1.5,colour="white") + gg <- gg + geom_segment(data=subset(ggTableResum,groupe!=nomGenerique),aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.5,size = 1.5,colour="white") + } + } else { + gg <- gg + geom_segment(data=subset(ggTableResum,groupe==nomGenerique),aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.8,size = 2.5) + gg <- gg + geom_segment(data=subset(ggTableResum,groupe==nomGenerique),aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.8,size = 2.5) + } + + gg <- gg + geom_point(data=d.reseau,size=2) + gg <- gg + labs(list(title=figure_titre,x=figure_abscisse,y=figure_ordonnee)) + + if(!is.null(colourProtocole)){ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA), axis.ticks = element_line(colour = colourProtocole, size = 1), axis.ticks.length = unit(0.3, "cm"),plot.title = element_text(colour = colourProtocole)) + }else{ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA)) + } + + if(!(is.null(tab_figure_couleur))){ + gg <- gg + scale_colour_manual(values = figure_couleur,name = "") + scale_fill_manual(values = figure_couleur,name = "",guide=FALSE)} + + ggfile <- paste(nomRep,nomProtocole,"_",m,".png",sep="") + cat("Check",ggfile,":") + ggsave(ggfile,gg) + cat("\n") + flush.console() + } +} + + +############################################################## +gglocal <- function(d, + col_abscisse = "AB_MOYENNE", + figure_abscisse = "Abondance", + col_ordonnee = "DIVERSITE_MOYENNE", + figure_ordonnee = "Diversite", + figure_titre = "Graphe referentiel", + col_filtre = "NOM_RESEAU", + nomGenerique = "GLOBAL", + col_sousGroup = "PARCELLEID", + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= NULL, + colourProtocole = NULL, + nomProtocole = NULL, + couleurLocal="#f609c1", + bagplot = TRUE, + bagProp=c(.05,.5,.95), + seuilSegment=30, + segmentSousSeuil=TRUE, + forcageMajusculeFiltre=TRUE, + forcageMajusculeSousGroupe=TRUE) { + + d$groupe <- d[,col_filtre] + d$abscisse <- d[,col_abscisse] + d$ordonnee <- d[,col_ordonnee] + d$sousGroup <- d[,col_sousGroup] + d$groupe <-gsub("/","_",d$groupe) + d$groupe <-gsub("!","",d$groupe) + d$sousGroup <-gsub("/","_",d$sousGroup) + d$sousGroup <-gsub("!","",d$sousGroup) + if(forcageMajusculeFiltre){ + d$groupe <- toupper(d$groupe)} + if(forcageMajusculeSousGroupe){ + d$sousGroup <- toupper(d$sousGroup)} + d <- subset(d,!(is.na(groupe)) & !(is.na(sousGroup)) & !(is.na(abscisse)) & !(is.na(ordonnee)) & groupe != "") + vecSousGroup <- as.character(unique(d$sousGroup)) + if(is.null(val_filtre)){ + lesModalites <- unique(d$groupe)} + else{ lesModalites <- val_filtre} + repResult <- dir("resultats/") +# if(!(col_filtre %in% repResult)){ +# dir.create(paste("resultats/",col_filtre,sep=""))} +# nomRep1 <- paste("resultats/",col_filtre,"/",sep="") + d.autre <- d + d.autre$groupe <- nomGenerique + for(m in lesModalites) { + d.reseau <- subset(d,groupe==m) + d.reseau$groupe <- m + ggTable <- rbind(d.autre,d.reseau) + seuilResum <- nrow(d.reseau) >= seuilSegment + ggTableResum <- aggregate(cbind(ordonnee, abscisse) ~ groupe, data = ggTable,quantile, c(.25,.5,.75)) + ggTableResum <- data.frame(ggTableResum[,1],ggTableResum[,2][,1:3],ggTableResum[,3][,1:3]) + colnames(ggTableResum) <- c("groupe","ordonnee.inf","ordonnee.med","ordonnee.sup","abscisse.inf","abscisse.med","abscisse.sup") + if(ggTableResum$groupe[2]==nomGenerique){ + ggTableResum <- ggTableResum[c(2,1),]} + if(!(is.null(tab_figure_couleur))) { + if(m %in% tab_figure_couleur$Modalite) { + figure_couleur <- setNames(c(as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == nomGenerique]), + as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == m]),couleurLocal), + c(nomGenerique,m,"")) + } else { + figure_couleur <- setNames(c(as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == nomGenerique]), + as.character(tab_figure_couleur$couleur[tab_figure_couleur$Modalite == ""]),couleurLocal), + c(nomGenerique,m,"")) + } + } +# repResult <- dir(nomRep1) +# if(!(m %in% repResult)){ +# dir.create(paste(nomRep1,m,sep=""))} +# nomRep <- paste(nomRep1,m,"/",sep="") +# if(!is.null(nomProtocole)) { +# repResult <- dir(nomRep) +# if(!(nomProtocole %in% repResult)){ +# dir.create(paste(nomRep,nomProtocole,sep=""))} +# nomRep <- paste(nomRep,nomProtocole,"/",sep="") +# } + d.reseau <- subset(d.reseau, !(is.na(sousGroup))) + figure_size<- setNames(c(1,3,2.5), c(nomGenerique,m,"")) + figure_shape<- setNames(c(16,16,20), c(nomGenerique,m,"")) + vecSousGroup <- as.character(unique(d.reseau$sousGroup)) + for(p in vecSousGroup) { + dp <- subset(d.reseau,sousGroup == p) + dp$groupe <- dp$sousGroup + ggTableSous <- rbind(d.reseau,dp) + ggTableSous <- rbind(d.autre,d.reseau,dp) + names(figure_couleur)[3] <- p + names(figure_shape)[3] <- p + names(figure_size)[3] <- p + gg <- ggplot(ggTableSous,aes(x=abscisse,y=ordonnee,colour=groupe,fill=groupe,shape=groupe,size=groupe)) + if(bagplot){ + gg <- gg + stat_bag(data=d.autre,prop=bagProp[1],colour=NA,alpha=.7) + stat_bag(data=d.autre,prop=bagProp[2],colour=NA,alpha=.4) + stat_bag(data=d.autre,prop=bagProp[3],colour=NA,alpha=.2) + }else{ + gg <- gg + geom_point(alpha=.2)} + gg <- gg + geom_hline(data=subset(ggTableResum,groupe == nomGenerique),aes(yintercept = ordonnee.med,colour=groupe),size=.5,linetype="dashed") + gg <- gg + geom_vline(data=subset(ggTableResum,groupe == nomGenerique),aes(xintercept = abscisse.med,colour=groupe),size=.5,linetype="dashed") + if(segmentSousSeuil) { + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.8,size=2.5) + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.8,size=2.5) + if(!(seuilResum)) { + gg <- gg + geom_segment(data=subset(ggTableResum,groupe!=nomGenerique),aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.5,size = 1.5,colour="white") + gg <- gg + geom_segment(data=subset(ggTableResum,groupe!=nomGenerique),aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.5,size = 1.5,colour="white") + } + } else { + gg <- gg + geom_segment(data=subset(ggTableResum,groupe==nomGenerique),aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.8,size = 2.5) + gg <- gg + geom_segment(data=subset(ggTableResum,groupe==nomGenerique),aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.8,size = 2.5) + } + gg <- gg + geom_point(data=subset(ggTableSous,groupe != nomGenerique)) + if(!(is.null(tab_figure_couleur))){ + gg <- gg + scale_colour_manual(values = figure_couleur,name = "") + scale_fill_manual(values = figure_couleur,name = "",guide=FALSE)} + gg <- gg + scale_shape_manual(values = figure_shape,name = "",guide=FALSE) + scale_size_manual(values = figure_size,guide=FALSE) + gg <- gg + labs(list(title=figure_titre,x=figure_abscisse,y=figure_ordonnee)) + if(!is.null(colourProtocole)){ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA), axis.ticks = element_line(colour = colourProtocole, size = 1), axis.ticks.length = unit(0.3, "cm"),plot.title = element_text(colour = colourProtocole)) } + else{ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA))} + ggfile <- paste(nomRep,nomProtocole,"_",m,"-",p,".png",sep="") + cat("Check",ggfile,":") + ggsave(ggfile,gg) + cat("\n") + flush.console() + } + } +} + + + +##################################################### +ggCompareLevel <- function(d, + col_abscisse = "abond_moyenne", + figure_abscisse = "Abondance", + col_ordonnee = "diversite_moyenne", + figure_ordonnee = "Diversite", + figure_titre = "Rhooo il dechire ce graphe", + col_filtre = "nom_reseau", + nomGenerique = "Global", + val_filtre = NULL, + figure_nom_filtre = NULL, + tab_figure_couleur= NULL, + colourProtocole = NULL, + nomProtocole = NULL, + bagplot = TRUE, + bagProp=c(.05,.5,.95), + seuilSegment=30, + segmentSousSeuil=FALSE, + forcageMajusculeFiltre=TRUE){ + + d$groupe <- d[,col_filtre] + d$abscisse <- d[,col_abscisse] + d$ordonnee <- d[,col_ordonnee] + d$groupe <-gsub("/","_",d$groupe) + d$groupe <-gsub("!","",d$groupe) + + if(forcageMajusculeFiltre){ + d$groupe <- toupper(d$groupe)} + d <- subset(d,!(is.na(groupe)) & !(is.na(abscisse)) & !(is.na(ordonnee)) & groupe != "") + if(is.null(val_filtre)){ + lesModalites <- unique(d$groupe) + }else{ + lesModalites <- val_filtre + } +# repResult <- dir("resultats/") +# if(!(col_filtre %in% repResult)){ +# dir.create(paste("resultats/",col_filtre,sep="")) +# } +# if(!is.null(nomProtocole)){ +# repResult <- dir(paste("resultats/",col_filtre,sep="")) +# if(!(nomProtocole %in% repResult)){ +# dir.create(paste("resultats/",col_filtre,"/",nomProtocole,sep=""))} +# nomRep <- paste("resultats/",col_filtre,"/",nomProtocole,"/",sep="") +# }else{ +# nomRep <- paste("resultats/",col_filtre,"/",sep="") +# } + d.autre <- d + d.autre$groupe <- nomGenerique + d.reseau <- subset(d,groupe %in% lesModalites) + ggTable <- rbind(d.autre,d.reseau) + ggTableResum <- aggregate(cbind(ordonnee, abscisse) ~ groupe, data = ggTable,quantile, c(.25,.5,.75)) + ggTableResum <- data.frame(ggTableResum[,1],ggTableResum[,2][,1:3],ggTableResum[,3][,1:3]) + colnames(ggTableResum) <- c("groupe","ordonnee.inf","ordonnee.med","ordonnee.sup","abscisse.inf","abscisse.med","abscisse.sup") + ggSeuil <- aggregate(ordonnee ~ groupe, data=ggTable,length) + ggSeuil$seuilResum <- ggSeuil$ordonnee >= seuilSegment + colnames(ggSeuil)[ncol(ggSeuil)] <- "seuil" + ggTableResum <- merge(ggTableResum,ggSeuil,by="groupe") + t_figure_couleur <- subset(tab_figure_couleur,Modalite %in% c(nomGenerique,lesModalites)) + modaliteSansCouleur <- lesModalites[(!(lesModalites %in% t_figure_couleur$Modalite))] + nbNxCol <- length(modaliteSansCouleur) + mypalette<-brewer.pal(nbNxCol,"YlGnBu") + figure_couleur <- setNames(c(as.character(t_figure_couleur$couleur),mypalette),c(as.character(t_figure_couleur$Modalite),modaliteSansCouleur)) + tab_coul <- data.frame(groupe=names(figure_couleur),couleur=figure_couleur) + tab_coul <- merge(tab_coul,ggTableResum,"groupe") + tab_coul$nom <- paste(tab_coul$groupe," (",tab_coul$ordonnee,")",sep="") + figure_couleur <- setNames(as.character(tab_coul$couleur),tab_coul$groupe) + figure_couleur_nom<- tab_coul$nom + gg <- ggplot(ggTable,aes(x=abscisse,y=ordonnee,colour=groupe,fill=groupe)) + if(bagplot){ + gg <- gg + stat_bag(data=d.autre,prop=bagProp[1],colour=NA,alpha=.7) + stat_bag(data=d.autre,prop=bagProp[2],colour=NA,alpha=.4) + stat_bag(data=d.autre,prop=bagProp[3],colour=NA,alpha=.2) + }else{ + gg <- gg + geom_point(alpha=.2) + } + gg <- gg + geom_hline(data=subset(ggTableResum,groupe=="Autre"),aes(yintercept = ordonnee.med,colour=groupe),size=.5,linetype="dashed") + geom_vline(data=subset(ggTableResum,groupe=="Autre"),aes(xintercept = abscisse.med,colour=groupe),size=.5,linetype="dashed") + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.7,size = 2.5) + gg <- gg + geom_segment(data=ggTableResum,aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.7,size = 2.5) + if(any(ggTableResum$seuil)){ + gg <- gg + geom_segment(data=subset(ggTableResum,!(seuil)),aes(x = abscisse.med, y = ordonnee.inf, xend = abscisse.med, yend = ordonnee.sup),alpha=.5,size = 1.5,colour="white") + gg <- gg + geom_segment(data=subset(ggTableResum,!(seuil)),aes(x = abscisse.inf, y = ordonnee.med, xend = abscisse.sup, yend = ordonnee.med),alpha=.5,size = 1.5,colour="white") + } + + #browser() # gg <- gg + geom_point(data=d.reseau,size=2) + gg <- gg + scale_colour_manual(values = figure_couleur,name = "",labels = figure_couleur_nom) + scale_fill_manual(values = figure_couleur,name = "",guide=FALSE) + gg <- gg + labs(list(title=figure_titre,x=figure_abscisse,y=figure_ordonnee)) + if(!is.null(colourProtocole)){ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA), axis.ticks = element_line(colour = colourProtocole, size = 1), axis.ticks.length = unit(0.3, "cm"),plot.title = element_text(colour = colourProtocole)) + }else{ + gg <- gg + theme(legend.justification=c(1,0), legend.position=c(1,0),legend.text = element_text(size = 7),legend.background = element_rect(fill=NA)) + } + ggfile <- paste(nomRep,nomProtocole,"_",col_filtre,"_","comparaison.png",sep="") + cat("Check",ggfile,":") + ggsave(ggfile,gg) + cat("\n") +flush.console() +} + + +######################################### + +#Lancement des fonctions : + + #Variables a definir : + +#filename="BDD_PAPILLONS_2016.txt" +#color_filename<-"code_couleurs.csv" + + #func +#func="ggCompareLevel" +#func="ggfiltre1niveau" +#func="gglocal" + + #colSousGroupe +#col_sousGroup_usr = NULL #ggfiltre #ggCompareLevel +#col_sousGroup_usr = "PARCELLENOM" #gglocal + + #vec_col_filtre_usr +#vec_col_filtre_usr = c("CONDUITEPARCELLE") #ggCompareLevel +#vec_col_filtre_usr = c("REGION") #ggfiltre +#vec_col_filtre_usr = c("NOM_RESEAU") #gglocal + + + +#Exe fonction : + +#filtre1niveau(func=func,nom_fichier=filename,nom_fichierCouleur=color_filename,col_sousGroup=NULL) #ggfiltre ou ggCompareLevel, depend de func et de vec_col_filtre_usr +#filtre1niveau(func=func,nom_fichier=filename,nom_fichierCouleur=color_filename,col_sousGroup = col_sousGroup_usr,vec_col_filtre=vec_col_filtre_usr) ## ==local + +######################################################## + +filename=args[2] +color_filename=args[3] +func=args[4] + +if(func=="ggCompareLevel"){ +col_sousGroup_usr=NULL +vec_col_filtre_usr=c("CONDUITEPARCELLE") +}else if(func=="ggfiltre1niveau"){ +col_sousGroup_usr=NULL +vec_col_filtre_usr=c("REGION") +}else if(func=="gglocal"){ +col_sousGroup_usr="PARCELLENOM" +vec_col_filtre_usr=c("NOM_RESEAU") +}else{ +#sortie erreur +write("Error, unknown function. Exit(1).", stderr()) +q('no') +} + +#create result dir +nomRep="resultats/" +dir.create(file.path(".", nomRep), showWarnings = FALSE) + + +filtre1niveau(func=func,nom_fichier=filename,nom_fichierCouleur=color_filename,col_sousGroup=col_sousGroup_usr,vec_col_filtre=vec_col_filtre_usr)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/butterfly_crossplot.xml Mon Aug 13 04:55:07 2018 -0400 @@ -0,0 +1,34 @@ +<tool id="Butterfly_crossplot" name="Butterfly data analysis" version="0.1.0"> + <requirements> + <requirement type="package" version="2.2.1">r-ggplot2</requirement> + <requirement type="package" version="1.1_2">r-rcolorbrewer</requirement> + </requirements> + <command detect_errors="exit_code"><![CDATA[ + Rscript '$__tool_directory__/butterfly_crossplot.R' '$__tool_directory__/stat_bag.r' '$butterfly_db' '$__tool_directory__/code_couleurs.csv' '$function' + ##'$__tool_directory__/BDD_PAPILLONS_2016.txt' + ]]></command> + <inputs> + <param name="function" type="select" label="Chose your analyse"> + <option value="ggCompareLevel">Compare production methods</option> + <option value="ggfiltre1niveau">Compare regions</option> + <option value="gglocal">Compare local networks</option> + </param> + <param name="butterfly_db" type="data" format="csv,tabular" label="Butterfly datafile"/> + </inputs> + <outputs> + <data format="png" name="ggCompareLevel" label="Compare production methods" from_work_dir="resultats/Papillons_CONDUITEPARCELLE_comparaison.png"> + <filter>function=='ggCompareLevel'</filter> + </data> + <collection type="list" name="output_ggfiltre" label="Compare regions"> + <discover_datasets pattern="__designation_and_ext__" visible="false" format="png" directory="resultats"/> + <filter>function=='ggfiltre1niveau'</filter> + </collection> + <collection type="list" name="output_gglocal" label="Compare local network"> + <discover_datasets pattern="__designation_and_ext__" visible="false" format="png" directory="resultats"/> + <filter>function=='gglocal'</filter> + </collection> + </outputs> + <help><![CDATA[ + TODO: Fill in help. + ]]></help> +</tool>
--- a/clim_data.R Mon Aug 13 04:47:57 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,136 +0,0 @@ -#!/usr/bin/env Rscript - - -args <- commandArgs(trailingOnly = TRUE) - -#Rscript clim_data.R 'worldclim' 'var' 'resolution' 'OutputFormat' #'FRA' #'prec1' -if (length(args)==0 | args[1]=="-h" | args[1]=="--help"){ - print ("general script execution : Rscript clim_data.R \'worldclim\' \'var\' resolution \'OutputFormat\' #\'variable-to-plot\' #\'region_code\'") - print ("eg : Rscript clim_data.R \'worldclim\' \'prec\' 10 \'raster\' #\'prec1' #\'FRA\'") - q('no') -} - - -#Climatic variables dictionaries -months<-c("January","February","March","April","May","June","July","August","September","October","November","December") -bioclimatic_vars<-c("Annual Mean Temperature", "Mean Diurnal Range (Mean of monthly (max temp - min temp))", "Isothermality (BIO2/BIO7) (x 100)","Temperature Seasonality (standard deviation x100)","Max Temperature of Warmest Month","Min Temperature of Coldest Month","Temperature Annual Range (BIO5-BIO6)","Mean Temperature of Wettest Quarter","Mean Temperature of Driest Quarter","Mean Temperature of Warmest Quarter","Mean Temperature of Coldest Quarter","Annual Precipitation","Precipitation of Wettest Month","Precipitation of Driest Month","Precipitation Seasonality (Coefficient of Variation)","Precipitation of Wettest Quarter","Precipitation of Driest Quarter","Precipitation of Warmest Quarter","Precipitation of Coldest Quarter") - -#Function to create custom plot title -get_plot_title<-function(usr_var,usr_var_to_plot){ - match<-str_extract(usr_var_to_plot,"[0-9]+") - if(usr_var %in% c("prec","tmin","tmax")){ - printable_var<-(months[as.integer(match)]) - if(usr_var=="prec"){ - printable_var<-paste(printable_var," precipitations (mm)",sep="") - }else if(usr_var=="tmin"){ - printable_var<-paste(printable_var," minimum temperature (°C *10)",sep="") - }else{ - printable_var<-paste(printable_var," maximum temperature (°C *10)",sep="") - } - }else if(usr_var=="bio"){ - printable_var<-(bioclimatic_vars[as.integer(match)]) - printable_var<-paste("Bioclimatic variable - ",printable_var,sep="") - } - title<-paste("Worldclim data - ",printable_var,".",sep="") - return(title) -} - - -#Call libraries -library('raster',quietly=TRUE) -library(sp,quietly = TRUE, warn.conflicts = FALSE) -library(ncdf4,quietly = TRUE, warn.conflicts = FALSE) -#library(rgdal,quietly = TRUE, warn.conflicts = FALSE) #To save as geotif -library(stringr) - - -#Get args -usr_data=args[1] -usr_var=args[2] -usr_res=as.numeric(args[3]) -usr_of=args[4] - - -# Retrieve 'var' data from WorldClim -global.var <- getData(usr_data, download = TRUE, var = usr_var, res = usr_res) - -# Check if we actualy get some -if (length(global.var)==0){ - cat("No data found.") -}else{ - writeRaster(global.var, "output_writeRaster", format=usr_of,overwrite=TRUE) - final_msg<-paste("WorldClim data for ", usr_var, " at resolution ", usr_res, " in ", usr_of, " format\n", sep="") - cat(final_msg) -} - - - - - - -################# -##Visualisation## -################# - -#Get args -if(length(args[5])>=0 && length(args[6])>=0){ - usr_var_to_plot=args[5] - usr_plot_region=args[6] -}else{q('no')} - -list_region_mask<-c("FRA","DEU","GBR","ESP","ITA") - -if(usr_plot_region %in% list_region_mask){ -#Country mask - region <- getData("GADM",country=usr_plot_region,level=0) - region_mask <- mask(global.var, region) - region_var_to_plot_expression<-paste("region_mask$",usr_var_to_plot,sep="") -}else{ #All map and resize manualy - region_var_to_plot_expression<-paste("global.var$",usr_var_to_plot,sep="") -} - - -region_var_to_plot<-eval(parse(text=region_var_to_plot_expression)) - -#PLotmap -jpeg(file="worldclim_plot_usr_region.jpeg",bg="white") - -title<-get_plot_title(usr_var,usr_var_to_plot) - - -if(usr_plot_region=="FRA"){ - #FRA - plot(region_var_to_plot, xlim = c(-7, 12), ylim = c(40, 52), axes=TRUE,xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="GBR"){ - #GBR - plot(region_var_to_plot, xlim = c(-10, 5), ylim = c(46, 63), axes=TRUE,xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="NA"){ - #North America : - plot(region_ar_to_plot,xlim=c(-180,-50),ylim=c(10,75),xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="EU"){ - #Europe - plot(region_var_to_plot,xlim=c(-28,48),ylim=c(34,72),xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="DEU"){ - #DEU - plot(region_var_to_plot, xlim = c(5, 15), ylim = c(45, 57),axes=TRUE,xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="ESP"){ - #ESP - plot(region_var_to_plot, xlim = c(-10, 6), ylim = c(35, 45), axes=TRUE,xlab="Longitude",ylab="Latitude", main=title) -}else if(usr_plot_region=="ITA"){ - #ITA - plot(region_var_to_plot, xlim = c(4, 20), ylim = c(35, 48), axes=TRUE,xlab="Longitude",ylab="Latitude", main=title) -}else if(usr_plot_region=="WM"){ - #Worldmap - plot(region_var_to_plot,xlab="Longitude",ylab="Latitude",main=title) -}else if(usr_plot_region=="AUS"){ - #AUS - plot(region_var_to_plot,xlim=c(110,155),ylim=c(-45,-10),xlab="Longitude",ylab="Latitude",axes=TRUE,main=title) -}else{ - write("Error with country code.", stderr()) - q('no') -} - -garbage_output<-dev.off - -#Exit -q('no')
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/code_couleurs.csv Mon Aug 13 04:55:07 2018 -0400 @@ -0,0 +1,24 @@ +Filtre;Modalite;couleur +TYPECULTURE;ARBORICULTURE;"#006600" +TYPECULTURE;AUTRE_CULTURE_PERENNE;"#FF5050" +TYPECULTURE;GRANDE_CULTURE;"#FFC000" +TYPECULTURE;MARAICHAGE;"#FF6600" +TYPECULTURE;PRAIRIE;"#97B314" +TYPECULTURE;VITICULTURE;"#660066" +NOM_RESEAU;;"#31849B" +Protocole;Abeilles;"#EE7F00" +Protocole;Invertebres;"#97B314" +Protocole;Papillons;"#009EE0" +Protocole;Vers de terre;"#666666" +TYPECULTURE;GLOBAL;"#9f0d0d" +NOM_RESEAU;GLOBAL;"#9f0d0d" +REGIONNAME;;"#31849B" +REGIONNAME;GLOBAL;"#9f0d0d" +ANNEE;;"#31849B" +ANNEE;GLOBAL;"#9f0d0d" +REGION;GLOBAL;"#9f0d0d" +REGION;;"#666666" +CONDUITEPARCELLE;;"#31849B" +CONDUITEPARCELLE;GLOBAL;"#9f0d0d" +TRAVAILSOL;;"#31849B" +TRAVAILSOL;GLOBAL;"#9f0d0d" \ No newline at end of file
--- a/raster_getdata.xml Mon Aug 13 04:47:57 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,208 +0,0 @@ -<tool id="raster_getdata" name="Get climatic data from Worldclim" version="0.1.0"> - <requirements> - -<!-- - <requirement type="package" >libgcc-ng</requirement> - <requirement type="package" >libiconv</requirement> - <requirement type="package" >jpeg</requirement> - <requirement type="package" >libstdcxx-ng</requirement> - <requirement type="package" >libkml</requirement> - <requirement type="package" version="1.2_18">r-rgdal</requirement> ---> - - <requirement type="package" version="1.20.2">r-getopt</requirement> - <requirement type="package" version="1.1.6">r-ncdf4</requirement> - <requirement type="package" version="2.6_7">r-raster</requirement> - <requirement type="package" version="1.3.1">r-stringr</requirement> - - </requirements> - <command detect_errors="exit_code"><![CDATA[ - Rscript - '$__tool_directory__/clim_data.R' - 'worldclim' - '$usr_var' - $usr_res - '$usr_of' - - #if $condi_plot.plot=='yes_plot' - '$usr_var$condi_plot.var_to_plot' - '$condi_plot.usr_country' - #end if - - ]]></command> - <inputs> - <param name="usr_var" type="select" label="Variable" help="Temperature minimum, maximum, precipitations and bioclimatic variables."> - <option value="tmin">tmin</option> - <option value="tmax">tmax</option> - <option value="prec">prec</option> - <option value="bio">bio</option> - </param> - <param name="usr_res" type="select" label="Resolution" help="In minutes of a degree."> - <option value="2.5">2.5</option> - <option value="5">5</option> - <option value="10">10</option> - </param> - <param name="usr_of" type="select" label="Output raster format" > - <option value="raster">Native raster package format .grd</option> - <option value="CDF">netCDF .nc</option> -<!-- - <option value="GTiff">GeoTiff .tif</option> - <option value="ascii">ESRI Ascii .asc</option> - <option value="SAGA">SAGA GIS .sdat</option> - <option value="IDRISI">IDRISI .rst</option> - <option value="ENVI">ENVI .hdr Labelled .envi</option> - <option value="EHdr">ESRI .hdr Labelled .bil</option> - <option value="HFA">Erdas Imagine Images .img</option> ---> - - </param> - <conditional name="condi_plot"> - <param name="plot" type="select" label="Return a layers visualisation"> - <option value="yes_plot">Yes</option> - <option value="no_plot">No</option> - </param> - <when value="yes_plot"> - <param name="var_to_plot" type="select" label="Select what you want to visualize" help="Be careful to select a layer that is present in the dowloaded raster. eg : bioclimatic variable (bio) and layer bio:Isothermality."> - <option value="1">[tmin|tmax|prec]:January</option> - <option value="2">[tmin|tmax|prec]:February</option> - <option value="3">[tmin|tmax|prec]:March</option> - <option value="4">[tmin|tmax|prec]:April</option> - <option value="5">[tmin|tmax|prec]:May</option> - <option value="6">[tmin|tmax|prec]:June</option> - <option value="7">[tmin|tmax|prec]:July</option> - <option value="8">[tmin|tmax|prec]:August</option> - <option value="9">[tmin|tmax|prec]:September</option> - <option value="10">[tmin|tmax|prec]:October</option> - <option value="11">[tmin|tmax|prec]:November</option> - <option value="12">[tmin|tmax|prec]:December</option> - <option value="1">bio:Annual Mean Temperature</option> - <option value="2">bio:Mean Diurnal Range (Mean of monthly (max temp - min temp))</option> - <option value="3">bio:Isothermality (BIO2/BIO7) (* 100)</option> - <option value="4">bio:Temperature Seasonality (standard deviation *100)</option> - <option value="5">bio:Max Temperature of Warmest Month</option> - <option value="6">bio:Min Temperature of Coldest Month</option> - <option value="7">bio:Temperature Annual Range (BIO5-BIO6)</option> - <option value="8">bio:Mean Temperature of Wettest Quarter</option> - <option value="9">bio:Mean Temperature of Driest Quarter</option> - <option value="10">bio:Mean Temperature of Warmest Quarter</option> - <option value="11">bio:Mean Temperature of Coldest Quarter</option> - <option value="12">bio:Annual Precipitation</option> - <option value="13">bio:Precipitation of Wettest Month</option> - <option value="14">bio:Precipitation of Driest Month</option> - <option value="15">bio:Precipitation Seasonality (Coefficient of Variation)</option> - <option value="16">bio:Precipitation of Wettest Quarter</option> - <option value="17">bio:Precipitation of Driest Quarter</option> - <option value="18">bio:Precipitation of Warmest Quarter</option> - <option value="19">bio:Precipitation of Coldest Quarter</option> - </param> - - <param name="usr_country" type="select" label="Select a worldmap vision or a more precise region" > - <option value="WM">Worldmap</option> - <option value="NA">North America</option> - <option value="EU">Europe</option> - <option value="AUS">Australia</option> - <option value="FRA">France</option> - <option value="DEU">Germany</option> - <option value="GBR">United Kingdom</option> - <option value="ESP">Spain</option> - <option value="ITA">Italy</option> - </param> - </when> - <when value="no_plot"> - </when> - </conditional> - </inputs> - <outputs> - <data name="output_df" from_work_dir="output_writeRaster.nc" label="WorldClim data - NetCDF format"> - <filter> usr_of=='CDF'</filter> - </data> - <data name="output_grd" from_work_dir="output_writeRaster.grd" label="WorldClim data - Native raster pkg format .grd"> - <filter> usr_of=='raster'</filter> - </data> - <data name="output_gri" from_work_dir="output_writeRaster.gri" label="WorldClim data - Native raster pkg format .gri"> - <filter> usr_of=='raster'</filter> - </data> - <data name="plot_jpeg" from_work_dir="worldclim_plot_usr_region.jpeg" format="jpg" label="WorldClim data - Visualisation"> - <filter> plot=='yes_plot'</filter> - </data> - </outputs> - <help><![CDATA[ - -.. class:: infomark - -======================================= -Get global climate data from Worldclim -======================================= - -| - -**What it does** - -This tool retrieve global climate data from Wordclim using the R raster package. - -The climatic data available are : - - * Minimum temperature (°C x10) - - * Maximum temperature (°C x10) - - * Precipitation (mm) - - * Bioclimatic variables : - - * BIO1 = Annual Mean Temperature - - * BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp)) - - * BIO3 = Isothermality (BIO2/BIO7) (x 100) - - * BIO4 = Temperature Seasonality (standard deviation x100) - - * BIO5 = Max Temperature of Warmest Month - - * BIO6 = Min Temperature of Coldest Month - - * BIO7 = Temperature Annual Range (BIO5-BIO6) - - * BIO8 = Mean Temperature of Wettest Quarter - - * BIO9 = Mean Temperature of Driest Quarter - - * BIO10 = Mean Temperature of Warmest Quarter - - * BIO11 = Mean Temperature of Coldest Quarter - - * BIO12 = Annual Precipitation - - * BIO13 = Precipitation of Wettest Month - - * BIO14 = Precipitation of Driest Month - - * BIO15 = Precipitation Seasonality (Coefficient of Variation) - - * BIO16 = Precipitation of Wettest Quarter - - * BIO17 = Precipitation of Driest Quarter - - * BIO18 = Precipitation of Warmest Quarter - - * BIO19 = Precipitation of Coldest Quarter - -| - -**How to use it** - -Simply select the variable to return and the resolution : 2.5, 5 or 10 minutes of a degree. - -Available output format are : netcdf and grd (Native raster package format). - -| - -**Sources** - -Worldclim : http://www.worldclim.org/ - -Raster package : https://cran.r-project.org/web/packages/raster/index.html - - ]]></help> -</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stat_bag.r Mon Aug 13 04:55:07 2018 -0400 @@ -0,0 +1,145 @@ +library(ggplot2) +StatBag <- ggproto("Statbag", Stat, + compute_group = function(data, scales, prop = 0.5) { + + ################################# + ################################# + # originally from aplpack package, plotting functions removed + plothulls_ <- function(x, y, fraction, n.hull = 1, + col.hull, lty.hull, lwd.hull, density=0, ...){ + # function for data peeling: + # x,y : data + # fraction.in.inner.hull : max percentage of points within the hull to be drawn + # n.hull : number of hulls to be plotted (if there is no fractiion argument) + # col.hull, lty.hull, lwd.hull : style of hull line + # plotting bits have been removed, BM 160321 + # pw 130524 + if(ncol(x) == 2){ y <- x[,2]; x <- x[,1] } + n <- length(x) + if(!missing(fraction)) { # find special hull + n.hull <- 1 + if(missing(col.hull)) col.hull <- 1 + if(missing(lty.hull)) lty.hull <- 1 + if(missing(lwd.hull)) lwd.hull <- 1 + x.old <- x; y.old <- y + idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx] + for( i in 1:(length(x)/3)){ + x <- x[-idx]; y <- y[-idx] + if( (length(x)/n) < fraction ){ + return(cbind(x.hull,y.hull)) + } + idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx]; + } + } + if(missing(col.hull)) col.hull <- 1:n.hull + if(length(col.hull)) col.hull <- rep(col.hull,n.hull) + if(missing(lty.hull)) lty.hull <- 1:n.hull + if(length(lty.hull)) lty.hull <- rep(lty.hull,n.hull) + if(missing(lwd.hull)) lwd.hull <- 1 + if(length(lwd.hull)) lwd.hull <- rep(lwd.hull,n.hull) + result <- NULL + for( i in 1:n.hull){ + idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx] + result <- c(result, list( cbind(x.hull,y.hull) )) + x <- x[-idx]; y <- y[-idx] + if(0 == length(x)) return(result) + } + result + } # end of definition of plothulls + ################################# + + + # prepare data to go into function below + the_matrix <- matrix(data = c(data$x, data$y), ncol = 2) + + # get data out of function as df with names + setNames(data.frame(plothulls_(the_matrix, fraction = prop)), nm = c("x", "y")) + # how can we get the hull and loop vertices passed on also? + }, + + required_aes = c("x", "y") +) + +#' @inheritParams ggplot2::stat_identity +#' @param prop Proportion of all the points to be included in the bag (default is 0.5) +stat_bag <- function(mapping = NULL, data = NULL, geom = "polygon", + position = "identity", na.rm = FALSE, show.legend = NA, + inherit.aes = TRUE, prop = 0.5, alpha = 0.3, ...) { + layer( + stat = StatBag, data = data, mapping = mapping, geom = geom, + position = position, show.legend = show.legend, inherit.aes = inherit.aes, + params = list(na.rm = na.rm, prop = prop, alpha = alpha, ...) + ) +} + + +geom_bag <- function(mapping = NULL, data = NULL, + stat = "identity", position = "identity", + prop = 0.5, + alpha = 0.3, + ..., + na.rm = FALSE, + show.legend = NA, + inherit.aes = TRUE) { + layer( + data = data, + mapping = mapping, + stat = StatBag, + geom = GeomBag, + position = position, + show.legend = show.legend, + inherit.aes = inherit.aes, + params = list( + na.rm = na.rm, + alpha = alpha, + prop = prop, + ... + ) + ) +} + +#' @rdname ggplot2-ggproto +#' @format NULL +#' @usage NULL +#' @export +GeomBag <- ggproto("GeomBag", Geom, + draw_group = function(data, panel_scales, coord) { + n <- nrow(data) + if (n == 1) return(zeroGrob()) + + munched <- coord_munch(coord, data, panel_scales) + # Sort by group to make sure that colors, fill, etc. come in same order + munched <- munched[order(munched$group), ] + + # For gpar(), there is one entry per polygon (not one entry per point). + # We'll pull the first value from each group, and assume all these values + # are the same within each group. + first_idx <- !duplicated(munched$group) + first_rows <- munched[first_idx, ] + + ggplot2:::ggname("geom_bag", + grid:::polygonGrob(munched$x, munched$y, default.units = "native", + id = munched$group, + gp = grid::gpar( + col = first_rows$colour, + fill = alpha(first_rows$fill, first_rows$alpha), + lwd = first_rows$size * .pt, + lty = first_rows$linetype + ) + ) + ) + + + }, + + default_aes = aes(colour = "NA", fill = "grey20", size = 0.5, linetype = 1, + alpha = NA, prop = 0.5), + + handle_na = function(data, params) { + data + }, + + required_aes = c("x", "y"), + + draw_key = draw_key_polygon +)