# User Request – Shepards Classification of Sediments

I received a request overnight on how to render the Shepard’s classification diagram, which is an alternative to the USDA’s textural soil classification. This is quite simple to produce (albeit a little tedious), however, before I walk through the script, immediately below, please see the final result (which you can compare to an original).

The diagram consists of 21 points, and 10 regions, in the following codes, we shall create a library of points and then map them to the polygons. In order to view how this was produced, please continue reading below…

### Preparing the Data.

Firstly, we need to create the dictionary of points.

#Build a library of points, left to right, top to bottom... points <- data.frame( rbind(c( 1,1.000,0.000,0.000), c( 2,0.750,0.250,0.000), c( 3,0.750,0.125,0.125), c( 4,0.750,0.000,0.250), c( 5,0.600,0.200,0.200), c( 6,0.500,0.500,0.000), c( 7,0.500,0.000,0.500), c( 8,0.400,0.400,0.200), c( 9,0.400,0.200,0.400), c(10,0.250,0.750,0.000), c(11,0.250,0.000,0.750), c(12,0.200,0.600,0.200), c(13,0.200,0.400,0.400), c(14,0.200,0.200,0.600), c(15,0.125,0.750,0.125), c(16,0.125,0.125,0.750), c(17,0.000,1.000,0.000), c(18,0.000,0.750,0.250), c(19,0.000,0.500,0.500), c(20,0.000,0.250,0.750), c(21,0.000,0.000,1.000) ) ) colnames(points) = c("IDPoint","T","L","R")

To visualise the points, the following code can be used:

base <- ggtern(data=points,aes(L,T,R)) + theme_bw() + theme_hidetitles() + theme_hidearrows() + geom_point(shape=21,size=10,color="blue",fill="white") + geom_text(aes(label=IDPoint),color="blue") print(base)

Which produces the following:

Assign each polygon a unique number and respective label.

#Give each Polygon a number polygon.labels <- data.frame( Label=c("Clay", "Sandy Clay", "Silty Clay", "Sand + Silt + Clay", "Clayey Sand", "Clayey Silt", "Sand", "Silty Sand", "Sandy Silt", "Silt")) #Assign each label an index polygon.labels$IDLabel=1:nrow(polygon.labels)

Create the map between the polygon numbers and the points which make up those numbers. Make sure they are in clockwise or anticlockwise order (but not mixed)

#Create a map of polygons to points polygons <- data.frame( rbind(c(1,1),c(1,2),c(1,4), c(2,6),c(2,2),c(2,3),c(2,5),c(2,8), c(3,3),c(3,4),c(3,7),c(3,9),c(3,5), c(4,5),c(4,14),c(4,12), c(5,6),c(5,8),c(5,12),c(5,15),c(5,10), c(6,7),c(6,11),c(6,16),c(6,14),c(6,9), c(7,17),c(7,10),c(7,18), c(8,15),c(8,12),c(8,13),c(8,19),c(8,18), c(9,13),c(9,14),c(9,16),c(9,20),c(9,19), c(10,11),c(10,21),c(10,20) ) ) #IMPORTANT FOR CORRECT ORDERING. polygons$PointOrder <- 1:nrow(polygons) #Rename the columns colnames(polygons) = c("IDLabel","IDPoint","PointOrder")

Now we merge the polygons, points and polygon labels to create a master dataframe.

#Merge the three sets together to create a master set. df <- merge(polygons,points) df <- merge(df,polygon.labels) df <- df[order(df$PointOrder),]

We also create a separate data frame for the labels positioned at the centroid of each polygon.

#Determine the Labels Data library(plyr) Labs = ddply(df,"Label",function(x){c(c(mean(x$T),mean(x$L),mean(x$R)))}) colnames(Labs) = c("Label","T","L","R")

This concludes the data preparation step.

### Constructing the Actual Plot

Now we can build the final plot, which employs the **geom_polygon(…)** and **geom_text(…)** geometries and the above data-sets, we apply some transparency so the grid can be seen through the polygons, and base the drawing of the simple **theme_bw(…)** arrangement.

#Build the final plot library(ggtern) base <- ggtern(data=df,aes(L,T,R)) + geom_polygon(aes(fill=Label,group=Label),color="black",alpha=0.25) + geom_text(data=Labs,aes(label=Label),size=4,color="black") + theme_bw() + custom_percent("Percent") + labs(title="Shepard Sediment Classification Diagram", fill = "Classification", T="Clay", L="Sand", R="Silt") print(base) #to console

Finally, if one likes, we can also render it directly to an image.

#Render to file. png("plot.png",width=800,height=600) print(base) dev.off()

To download the full script, please click HERE.

Thanks for the help this is really useful to me

Hm… package ‘gtern’ is not available (for R version 3.1.1). ^_^

Anyway, thanks. Seems like I could use this. Have to look at your site to understand how to put int he values in the end. Never considered plotting my soil data into this via R, for some reason. Strange, isn’t it?

Clearly the binary hasn’t been built for your R version, try: install.packages(“ggtern”,type=”source”)

Is there any simple way to change the filling colors manually?

Hy, use the scale_fill_manual(…) function, see here: http://docs.ggplot2.org/0.9.3.1/scale_manual.html

Very cool, Nick. Congrats from Rob & Pam

I’m new to ggtern. I’ve managed to make the diagram following your code (great stuff!!!) but I’m having trouble plotting my points on top of the diagram. I have a series of soils for which I have the clay, sand and silt %’s but can’t figure out how to plot them on top of this diagram. If you could lend a hand or point me into the right direction I would appreciate it a lot.

Thanks

Hi, I would appreciate it too 🙂

Thanks for the help.

Leandro, most of the basics for ggplot2 work for ggtern.

Using ‘base’ as the last plot from the script above, add some random points via local data.frame to a new geometry layer.

Dear Nicholas,

Thanks, it works great.

I try to apply a ‘stat_density_tern()’ script to ‘newPoints’ but it’s applied to ‘base’ instead.

How can I fix this?

Dear Nicholas,

Thanksn it works great.

I try to apply a ‘stat_density_tern()’ script to ‘newPoints’ but it’s applied to ‘base’.

How can I fix this?

Here’s my script :

Very useful example!

However, in ggtern 2.1.4, running the code for the first figure results in a figure that looks different: points on the edges of the equilateral triangle are cropped.

Turn the clipping mask off (theme_nomask()) or put the clipping mask (geom_mask()) UNDERNEATH the points layer.

Dear Nicholas,

I’m relative new with rstudio and ggtern and was wondering if you could help me figure out how to plot my own points on your Diagram. I have my own data set that I would really like to plot, but I keep getting error massages about not recognizing the Z aesthetics. Any help will be greatly appreciated.

Here’s my code:

setwd(“~/rminid/Data”)

dat <- read.csv("tertiary.csv", header = T)

head(dat)

install.packages("ggplot2")

install.packages("ggtern")

install.packages("plyr")

install.packages("grid")

library(ggplot2)

library(ggtern)

library(plyr)

library(grid)

points <- data.frame(rbind(c( 1,1.000,0.000,0.000),

+ c( 2,0.750,0.250,0.000),

+ c( 3,0.750,0.125,0.125),

+ c( 4,0.750,0.000,0.250),

+ c( 5,0.600,0.200,0.200),

+ c( 6,0.500,0.500,0.000),

+ c( 7,0.500,0.000,0.500),

+ c( 8,0.400,0.400,0.200),

+ c( 9,0.400,0.200,0.400),

+ c(10,0.250,0.750,0.000),

+ c(11,0.250,0.000,0.750),

+ c(12,0.200,0.600,0.200),

+ c(13,0.200,0.400,0.400),

+ c(14,0.200,0.200,0.600),

+ c(15,0.125,0.750,0.125),

+ c(16,0.125,0.125,0.750),

+ c(17,0.000,1.000,0.000),

+ c(18,0.000,0.750,0.250),

+ c(19,0.000,0.500,0.500),

+ c(20,0.000,0.250,0.750),

+ c(21,0.000,0.000,1.000)))

colnames(points) = c("IDPoint","T","L","R")

polygon.labels <- data.frame(Label=c("Silica dominated lithotype","Clay-rich siliceous mudstone","Carbonate-rich siliceous mudstone","Mixed mudstone","Silica-rich argillaceous mudstone","Silica-rich calcareous mudstone","Clay dominated lithotype","Carbonate-rich argillaceous mudstone","Clay-rich calcareous mudstone","Carbonate dominated lithotype"))

polygon.labels$IDLabel=1:nrow(polygon.labels)

polygons <- data.frame(rbind(c(1,1),c(1,2),c(1,4),c(2,6),c(2,2),c(2,3),c(2,5),c(2,8),c(3,3),c(3,4),c(3,7),c(3,9),c(3,5),c(4,5),c(4,14),c(4,12),c(5,6),c(5,8),c(5,12),c(5,15),c(5,10),c(6,7),c(6,11),c(6,16),c(6,14),c(6,9),c(7,17),c(7,10),c(7,18),c(8,15),c(8,12),c(8,13),c(8,19),c(8,18),c(9,13),c(9,14),c(9,16),c(9,20),c(9,19),c(10,11),c(10,21),c(10,20)))

polygons$PointOrder <- 1:nrow(polygons)

colnames(polygons) = c("IDLabel","IDPoint","PointOrder")

df <- merge(polygons,points)

df <- merge(df,polygon.labels)

df <- df[order(df$PointOrder),]

Labs = ddply(df,"Label",function(x){c(c(mean(x$T),mean(x$L),mean(x$R)))})

colnames(Labs) = c("Label","T","L","R")

base <- ggtern(data=df,aes(L,T,R)) + geom_polygon(aes(fill=Label,group=Label),color="black",alpha=0.25) + geom_text(data=Labs,aes(label=Label),size=2,color="black") + theme_bw() + custom_percent("Percent") + labs(title="Shepard Sediment Classification Diagram", fill = "Classification", T="Quartz", L="Total Clay",R="Calcium")

print(base)

Dear Nicholas,

I’m Trying to get a plot with only 4 areas (triangles).

But I can’t manage to get the center area.

Here is my script:

points <- data.frame(

rbind(c( 1,1.000,0.000,0.000),

c( 2,0.500,0.500,0.000),

c( 3,0.500,0.000,0.500),

c( 4,0.000,1.000,0.000),

c( 5,0.000,0.500,0.500),

c( 6,0.000,0.000,1.000)

)

)

colnames(points) = c(“IDPoint”,”T”,”L”,”R”)

#Give each Polygon a number

polygon.labels <- data.frame(Label=c(“O”,”ONP”,”N”,”P”))

polygon.labels$IDLabel=1:nrow(polygon.labels)

#Create a map of polygons to points

polygons <- data.frame(

rbind(c(1,1),c(1,2),c(1,3),

c(2,2),c(2,4),c(2,5),

c(3,3),c(3,5),c(3,6)

)

)

Thanks for the help.

Dear Nicholas,

I am trying to get a plot with 4 areas (triangles).

I can’t manage to get the center area.

Here is my script:

points <- data.frame(

rbind(c( 1,1.000,0.000,0.000),

c( 2,0.500,0.500,0.000),

c( 3,0.500,0.000,0.500),

c( 4,0.000,1.000,0.000),

c( 5,0.000,0.500,0.500),

c( 6,0.000,0.000,1.000)

)

)

colnames(points) = c(“IDPoint”,”T”,”L”,”R”)

#Give each Polygon a number

polygon.labels <- data.frame(Label=c(“O”,”ONP”,”N”,”P”))

polygon.labels$IDLabel=1:nrow(polygon.labels)

#Create a map of polygons to points

polygons <- data.frame(

rbind(c(1,1),c(1,2),c(1,3),

c(2,2),c(2,4),c(2,5),

c(3,3),c(3,5),c(3,6)

)

)

Thanks for the help