R software

PACKAGE G2Sd

NEW !!!!!!

La version 2.1 du package G2Sd "Grain-size Statistics and Description of Sediment" est disponible sur le CRAN.

Cette nouvelle version intègre la correction de bugs et UNE INTERFACE SIMPLIFIEE à l’aide du navigateur internet par défaut.
Le package peut être utilisé sans installer R en cliquant sur le lien : https://regisgallon.shinyapps.io/G2Sd/

Comment utiliser G2Sd?

G2Sd a été développé de manière à décrire le sédiment à partir de tamis utilisant la nomenclature métrique ou phi.

1. Mise en forme de la matrice granulo

Cette matrice doit comporter en colonne les stations échantillonnées et en ligne les ouvertures de tamis (microns ou phi)

>granulo
         Q1    Q2   Q3    Q4   Q5    Q6
25000  0.00  0.00 0.00  0.00 0.00  0.00
20000  0.00  0.00 0.00  0.00 0.00  0.00
16000  0.00  0.00 0.00  0.00 0.00  0.00
12500  0.00  0.00 0.00  0.00 0.00  0.00
10000  0.35  0.00 2.20  0.00 1.30  0.00
8000   0.00  0.25 0.00  0.00 0.30  0.00
6300   0.00  0.40 0.00  0.00 0.30  0.00
5000   0.00  0.15 0.00  0.00 0.20  0.00
4000   0.00  0.75 0.00  0.00 0.30  0.00
2500   1.00  0.75 0.30  3.30 1.40  0.00
2000   0.65  1.85 0.10  4.35 3.30  0.00
1600   0.60  2.70 0.40  5.80 5.20  0.00
1250   1.00  3.50 1.00  4.85 7.20  0.15
1000   0.70  3.60 1.30  2.60 6.50  0.10
800    0.80  3.30 1.60  2.60 5.60  0.10
630    0.70  2.70 1.70  2.25 4.30  0.20
500    1.00  3.15 2.45  3.05 4.60  0.30
400    0.85  2.30 2.10  2.50 2.95  0.35
315    1.10  2.30 2.30  2.75 2.55  0.50
250    1.70  2.40 2.70  3.15 2.30  0.95
200    2.40  1.90 2.60  2.70 1.90  1.40
160    3.15  1.50 2.70  2.60 1.95  2.50
125    3.10  1.10 2.40  2.00 1.65  3.30
100    2.95  1.00 2.10  1.75 1.65  5.25
80     3.40  1.65 2.10  1.90 1.80  7.10
63     3.45  2.60 1.30  1.80 1.60  4.80
50     2.15  1.40 0.70  1.00 0.90  2.20
40     0.15  0.05 0.10  0.10 0.00  0.20
0     18.65 19.75 1.90 10.35 5.85 21.80

Il est possible d’insérer une partie de la matrice granulométrique, elle sera intégrée à la matrice de calcul. La compatibilité se déroule en lançant la fonction granstat(), vous lirez alors "Compatibility progress…".  Votre matrice sera alors intégrée dans une matrice remplie de zéro avec 29 lignes et n colonnes (n étant le nombre de stations échantillonnées).

2. Utilisation de granstat()

Cette fonction est le cœur de calcul pour définir les caractéristiques des sédiments. Elle intègre trois statistiques descriptives :

  1. Méthodes arithmétiques
  2. Méthodes des moments
  3. Méthodes de Folk et Ward (1957)

Pour les détails des modes de calcul et l’interprétation des indices se référer à  Fournier et al. (2012).

NEW !!!! Une interface simplifiée est disponible en écrivant : granstat(web_interface=TRUE)

G2Sd_web

En utilisant aggr=F, vous pouvez sélectionner soit les statistiques ($stat), les indices ($index), les modes ($modes) ou les caractéristiques des sédiments ($sedim). Dans l’exemple, nous allons simplement nous intéresser à la méthode Folk et Ward (1957).

> granstat(granulo[,1:6],statistic="folk.ward", aggr=F)$index
                  Q1       Q2       Q3       Q4       Q5      Q6
D10(µm)        2.634    3.077   71.714    8.835   60.000   2.336
D50(µm)       82.805  237.828  275.271  422.006  748.405  67.546
D90(µm)      826.078 1701.610 1447.771 2161.076 2074.378 181.783
D90/D10      313.581  553.041   20.188  244.601   34.573  77.834
D90-D10      823.443 1698.533 1376.057 2152.241 2014.378 179.447
D75/D25       17.764   53.180    4.941   14.159    6.181  12.766
D75-D25      195.666  897.867  508.711 1349.931 1165.253 101.632
Trask(So)      4.215    7.292    2.223    3.763    2.486   3.573
Krumbein(Qd)   2.075    2.866    1.152    1.912    1.314   1.837
> granstat(granulo[,1:6],statistic="folk.ward", aggr=F)$sedim
                                   Q1                  Q2            Q3                  Q4                  Q5         Q6
Texture  Slightly Gravelly Muddy Sand Gravelly Muddy Sand Gravelly Sand Gravelly Muddy Sand Gravelly Muddy Sand Muddy Sand
Boulder                             0                   0             0                   0                   0          0
Gravel                          4.012               6.798         7.636              12.459              10.823          0
Sand                           53.962              58.477        84.435              68.893              78.887     52.734
mud                            42.026              34.726          7.93              18.648               10.29     47.266
boulder                             0                   0             0                   0                   0          0
vcgravel                            0                   0             0                   0                   0          0
cgravel                             0                   0             0                   0                   0          0
mgravel                         0.702                0.41         6.461                   0               2.439          0
fgravel                             0               2.129             0                   0                1.22          0
vfgravel                         3.31               4.259         1.175              12.459               7.165          0
vcsand                          4.614              16.052          7.93               21.58              28.811      0.488
csand                           5.015              14.988        16.887              12.866              22.104      1.172
msand                           7.322              11.466        20.852              13.681               11.89      3.516
fsand                          17.352               7.371        22.614              11.889               8.384     14.062
vfsand                         19.659                 8.6        16.153               8.876               7.698     33.496
vcsilt                          4.614               2.375         2.349               1.792               1.372      4.688
silt                           37.412              32.351          5.58              16.857               8.918     42.578

Le calcul des modes se fait à la main, vous devez définir modes=TRUE dans granstat(). Une fenêtre graphique va alors apparaitre, avec le pointeur de la souris vous devez définir les modes. Au maximum trois modes peuvent être définis. Si vous avez sélectionné moins de 3 modes vous devez appuyer sur "Echap" pour passer à la courbe suivante.

3. Fonctions graphiques granplot() et grandistrib()

La fonction granplot() permet de visualiser la distribution des sédiments.

 par(mfrow=c(2,2))
 granplot(granulo,xc=1,main=names(granulo)[1],hist=T,cum=T)
 granplot(granulo,xc=1,main=names(granulo)[1],hist=T,cum=F)
 granplot(granulo,xc=1,main=names(granulo)[1],hist=F,cum=T)

granplot_exemple1

En réalisant un boucle for, vous pouvez afficher la distribution des sédiments des stations échantillonnées :

par(mfrow=c(2,3))
for (i in 1:6)# on affichera que les 6 premières stations
{
granplot(granulo,xc=i,main=names(granulo)[i],hist=T,cum=T)
}

granplot2

La fonction grandistrib() permet de visualiser la composition du sédiment. Deux échelles peuvent être choisies :

 grandistrib(granulo[,1:6],scale="fine")
 grandistrib(granulo[,1:6],scale="large")

grandistrib

4. Fonction granmap() : non intégrée au package G2Sd

La fonction granmap() utilise le package AKIMA dont l’utilisation est restreinte par une licence. Cette fonction permet d’avoir une vision spatiale de la répartition des sédiments.

granmap <-  function(x,XY,points=FALSE,background=NULL,crs="+init=epsg:2154",main="",
pos.legend="bottomright",save=FALSE,filename="Rplot")
{
require(raster); require(sp); require(akima)
.mgrep <- function(mpattern,x,FUN="grep")
{
  select=NULL
  for (i in 1:length(mpattern))
  {
    if (FUN=="grep")
      values=grep(mpattern[i],x)
    if (FUN=="which")
      values=which(x==mpattern[i])  
    select=c(select,values)
  }
  return(sort(select))
}

XYcoord=SpatialPoints(XY,proj4string=CRS(crs))

Sedim <- granstat(x,aggr=F)$stat["Sediment",]
i=1
for (sedim in Sedim){
sedim=strsplit(sedim,NULL)
sedim <- paste0(sedim[[1]][1:(which(sedim[[1]]==",")[1])-1])
mot=NULL
for (l in sedim) mot=paste(mot,l,sep="")
Sedim[i] <- mot
i=i+1
}
color <- c("#330033","#660066","#990099","#CC00CC","#FF00FF",
"#330000","#660000","#990000","#CC0000","#FF0000","#FFCC00","#FFCC33","#FFCC66",
"#FFCC99","#FFCCCC","#006666")
Sedim_color <- cbind(description=c("Very Coarse Gravel" , "Coarse Gravel","Medium Gravel", "Fine Gravel",
"Very Fine Gravel","Very Coarse Sand","Coarse Sand","Medium Sand",
"Fine Sand","Very Fine Sand","Very Coarse Silt","Coarse Silt",
"Medium Silt","Fine Silt","Very Fine Silt","Clay"),1:16,col=color)

Sedim.df <- cbind(XY,t(Sedim[1,]))
Sedim.df$Sediment <- as.numeric(as.integer(Sedim.df$Sediment))
Sedim.df$Sediment <- as.numeric(Sedim_color[match(Sedim,Sedim_color[,1]),2])
out=interp(Sedim.df[,1],Sedim.df[,2],Sedim.df[,3],duplicate="mean")
out$z=round(out$z,0)
#   if (extract==FALSE)
#   {

classif <- cbind(t(Sedim[1,]),Sedim.df$Sediment)
classif <- unique(classif);classif=as.data.frame(classif)

classif <- cbind(classif,color=Sedim_color[.mgrep(classif[,1],Sedim_color[,1],FUN="which"),3])
classif[,2] <- as.numeric(as.character(classif[,2]))
labels.classif <- classif[order(classif[,2]),]
#       labels.classif[,2] <- 1:nrow(labels.classif)

if (save==FALSE){

layout(matrix(c(1,1,1,1,1,2), nrow = 1), widths = c(0.7, 0.3))
par(mar = c(5, 4, 4, 0))
image(out$x,out$y,out$z,col=as.character(classif[,3]),xlab="Longitude",ylab="Latitude",cex.lab=1,font.lab=2)
legend(pos.legend,legend=crs,bty="n",cex=0.8,text.font=3)
if (!is.null(background))
{
if (is.na(projection(background)))
{
warning("Please give projection to the background map.", call. = FALSE,immediate.=TRUE)
}

if (projection(XYcoord)!=projection(background))
{
warning("SpatialPoints and background have different projections.", call. = FALSE,immediate.=TRUE)
}
else{
plot(background,add=TRUE)
}
}
if (points==TRUE)
{
points(XY,pch=20)
text(XY,labels=row.names(XY),pos=4)}

par(mar = c(5, 0, 4, 0))
plot(1,1,type="n",axes=F,xlab="",ylab="")
legend("center", legend=labels.classif[,1], fill=as.character(classif[,3]), title="SEDIMENT",cex=1.1)

}
if (save==TRUE)
{
pdf(file=paste(filename,".pdf", sep = ""), width = 11.7, height = 8.3)

layout(matrix(c(1,1,1,1,1,2), nrow = 1), widths = c(0.7, 0.3))
par(mar = c(5, 4, 4, 0))
image(out$x,out$y,out$z,col=as.character(classif[,3]),xlab="Longitude",ylab="Latitude",cex.lab=1,font.lab=2)
legend(pos.legend,legend=crs,bty="n",cex=0.8,text.font=3)
if (!is.null(background))
{
if (is.na(projection(background)))
{
warning("Please give projection to the background map.", call. = FALSE,immediate.=TRUE)
}

if (projection(XYcoord)!=projection(background))
{
warning("SpatialPoints and background have different projections.", call. = FALSE,immediate.=TRUE)
}
else{

plot(background,add=TRUE)
}
}
if (points==TRUE)
{
points(XY,pch=20)
text(XY,labels=row.names(XY),pos=4)}

par(mar = c(5, 0, 4, 0))
plot(1,1,type="n",axes=F,xlab="",ylab="")
legend("center", legend=labels.classif[,1], fill=as.character(classif[,3]), title="SEDIMENT",cex=1.2)
dev.off()
}
}

Les arguments :

  • x : la matrice granulo
  • XY : les coordonnées spatiales des points échantillonnées
  • background : un objet de la classe raster (package raster)
  • crs : Système de projection (voir le global mapper)
  • pos.legend : Position de la légende
  • save : si TRUE sauve la carte au format pdf
  • filename : nom du fichier
granmap(granulo,coord_gran,points=TRUE,crs="+init=epsg:2154")

carto_sédiment

Pour tous renseignements ou remarques, n’hésitez pas à m’écrire.




Suivre

Recevez les nouvelles publications par mail.

%d bloggers like this: