This project is about carrying out sediment grain size analysis is R.

The code below is a function that converts string input to a dataframe table.

text_to_table <- function(text, ...)
{
  dfr <- read.table(tc <- textConnection(text), ...)
  close(tc)
  dfr
}

We now proceed to read in the dataset.

Sieves<-text_to_table('Phi  wt
-1.00 0.00
-0.24   23.71
0.00    0.01
0.49    0.18
0.74    2.35
1.00    2.51
1.25    0.53
1.74    6.43
2.00    6.22
3.16    45.5
3.74    9.92
3.99    1.7
5.00    0.77', header=T)
head(Sieves)
##     Phi    wt
## 1 -1.00  0.00
## 2 -0.24 23.71
## 3  0.00  0.01
## 4  0.49  0.18
## 5  0.74  2.35
## 6  1.00  2.51

Create other data columns: weight percentage(W_p), cumulative weight(Cum_w), and cumulative weight percentage(Cum_w_P ).

Sieves$W_p<-(Sieves$wt *100)/sum(Sieves$wt)
Sieves$Cum_w<-cumsum(Sieves$wt)
Sieves$Cum_w_P<-cumsum(Sieves$W_p)

We are going to fit a brezier function to the generate smooth cumulatve frequence graph rather than using simple points and connecting line. These brezier function is provided below

bezierCurve <- function(x, y, n=10)
{
  outx <- NULL
  outy <- NULL
  
  i <- 1
  for (t in seq(0, 1, length.out=n))
  {
    b <- bez(x, y, t)
    outx[i] <- b$x
    outy[i] <- b$y
    
    i <- i+1
  }
  
  return (list(x=outx, y=outy))
}

bez <- function(x, y, t)#smoothing function
{
  outx <- 0
  outy <- 0
  n <- length(x)-1
  for (i in 0:n)
  {
    outx <- outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1]
    outy <- outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1]
  }
  
  return (list(x=outx, y=outy))
}

The axis text of the bargraph showing the particle size distribution will be customized, and the Hmisc package will be used to place tick marks on the plots. The code chunck below implement it.

#install.packages('Hmisc')
#library(Hmisc)
par(pty="s") #Square plots
Xaxis_label<-c(-1,-0.24, 0.0,  0.49,  0.74,  1.00,  1.25, 1.74,  2.0,  3.16,  3.74,  3.99, "pan")
barplot(Sieves$wt,names=Xaxis_label, 
  xlab= expression("Sieve size  " *phi), 
  space = 0,
  ylim = c(0,50),
  yaxt='n',
  las=2,
  ylab = "% weight retained",
  col ="rosybrown" ,
  main = "Histogram showing particule size distribution")
axis(2, at=seq(0,40,10),labels=seq(0,40,10), col.axis="black", las=2)

Below is the code to plot the cumulative frequency curve. The phi values used to compute Folk and Ward (1957) parameters are measured by eyeballing to inserting lines; where horizontal and vertical lines meet at the cumulative frequency curve.

plot(Sieves$Phi, Sieves$Cum_w_P, type="n", 
     lty=1, pch=20 , cex=2,
     yaxt='n', 
     xaxt='n',
     xlim = c(-1, 5),
     ylim = c(0, 100),
     main = "Cumulative frequency curve",
     ylab = "Cumulative Weight Percentage (%)",
     xlab= expression("Sieve size  " *phi))
points(bezierCurve(Sieves$Phi,Sieves$Cum_w_P,20), type="l", col="black", lty = 1, lwd=2)
axis(2, at=seq(0,100,10), labels =seq(0,100,10), lty=2, las=2)
Xaxis_label<-c(-1,-0.5, 0.0,  0.5,  1.0,    1.5,  2.0,  2.5,  3.0,3.5,4.0, 4.5,  "pan")
att<-c(-1,-0.5, 0.0,  0.5,  1.0,    1.5,  2.0,  2.5,  3.0,3.5,4.0, 4.5, 5.0)
axis(1, at=att, labels =Xaxis_label, lty=2,las=2)
grid(col = 'grey70')
clip(-1,-0.81, 0, 5)
abline(h=5, v=-0.81, lw=1.5, lty=2, col='red')

clip(-1,-0.34, 0, 16)
abline(h=16, v=-0.34, lw=1.5, lty=2, col='red')      

clip(-1,0.52, 0, 25)
abline(h=25, v=0.52, lw=1.5, lty=2, col='red')      

clip(-1,2.14, 0, 50)
abline(h=50, v=2.14, lw=1.5, lty=2, col='red')  

clip(-1,3.05, 0, 75)
abline(h=75, v=3.05, lw=1.5, lty=2, col='red') 

clip(-1,3.40, 0, 84)
abline(h=84, v=3.40, lw=1.5, lty=2, col='red') 

clip(-1,4.0, 0, 95)
abline(h=95, v=4.0, lw=1.5, lty=2, col='red')

Finally the measure critical phi values is read from the v paramters of the abline function and summaried in a dataframe and exported for further graphical analysis.

Crit_Phi_values<-c(-0.81,-0.34,0.52,2.14,3.05,3.40,4.0)
Phi<-c('Phi5','Phi16','Phi25','Phi50','Phi75','Phi84','Phi95')
Result<- data.frame(Phi,Crit_Phi_values)
#write.table(Result, 'Result.txt',  sep = "\t")

I will appreciate contributions that help to automate this project.