# 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

Dear Nicholas,

I am trying to get a plot with 8 areas.

I can manage to get 7 areas but the last one doesn’t work.

Here is my script:

library(ggtern)

library(ggplot2)

library(plyr)

#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.950,0.050,0.000),

c( 3,0.900,0.050,0.050),

c( 4,0.950,0.000,0.050),

c( 5,0.800,0.100,0.100),

c( 6,0.750,0.250,0.000),

c( 7,0.500,0.250,0.250),

c( 8,0.750,0.000,0.250),

c( 9,0.000,0.100,0.000),

c(10,0.000,0.900,0.100),

c(11,0.000,0.500,0.500),

c(12,0.000,0.100,0.900),

c(13,0.000,0.000,0.100),

c(14,0.650,0.250,0.100),

c(15,0.650,0.100,0.250)

)

)

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

#Give each Polygon a number

polygon.labels <- data.frame(Label=c(“Quartzarenit”,”Subarkose”,”Silty Clay”,

“Sublitharenite”,”Lithicsubarkose”,”Arkose”,

“Lithicarkose”,”Litharenite”,”Sandy Silt”,”Silt”))

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

#Create a map of polygons to points

polygons <- data.frame(

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

c(2,2),c(2,3),c(2,5),c(2,14),c(2,6),

c(4,4),c(4,3),c(4,5),c(4,15),c(4,8),

c(5,5),c(5,14),c(5,7),c(5,15),

c(6,6),c(6,14),c(6,10),c(6,9),

c(8,8),c(8,15),c(8,12),c(8,13),

c(7,7),c(7,11),c(7,10),c(7,14)

)

)

polygons$PointOrder <- 1:nrow(polygons)

colnames(polygons) = c(“IDLabel”,”IDPoint”,”PointOrder”)

#Merge the three sets together to create a master set.

df <- merge(polygons,points)

df <- merge(df,polygon.labels)

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

#Determine the Labels Data

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

colnames(Labs) = c(“Label”,”T”,”L”,”R”)

#Build the final plot

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=”McBride 1963″,

fill = “Legende”,

T=”Quartz”,

L=”Feldspat”,

R=”Geisteinsfragmente”)

print(base) #to console

#Render to file.

png(“plot.png”,width=800,height=600)

print(base)

dev.off()

Thanks for the help

The best in class, Depoxito present you high-end experience that lecture to the see

and tone of valid VIP standarts, we have enough money you the best fascinating

to high-level experience of VIPs expect in any top stop casino, grand enliven casino royale present you the further

studio design element including the grand blackjack, offering our VIP

Customer the best experience of a Salon privee table.

New style table also feature across the room similar to

grand roulette upgraded upon our provider playtechs mini

prestige roulette which delivering more fascinating and richer playing experience.

The new experience contains a sum of seven tables including five blackjack tables, one roulette table and

one baccarat table. Grand liven up casino royale has been high hand-engineered to fit the

needs of our customer to using it, and contains unique elements

that is specially expected to maximize the impact value we got

from our customers and diversify it to the existing network.

Soon, Depoxito will fabricate an greater than before realism

technology on breathing casino for our VIP member, these most ahead of its time technology ever seen in conscious casino including this enlarged reality.

Which permit players to experience products on an entire new level which is never seen before literally leaping

out of the game and taking the blackjack, baccarat, roulette and

further game into the total entire level.

Depoxito VIP Baccarat, we come up with the money for you the entirely exclusive bring to life VIP Baccarat that is played gone

happening to 7 players at the thesame table and our extremely

trained lovely conscious baccarat dealer. And of course our VIP aficionado will mood

as if they were really sitting at one of the top casino baccarat table.

This immersive gaming experience creates a hugely thrill-seeking vent

that our VIP players will locate hard to surpass.

Here is the list of alive casino game that depoxito provide, we

present the widest range of sentient casino games on the shout from the rooftops including :

blackjack unlimited, blackjack prestige, roulette,

baccarat, poker, hi-lo, sic bo, and grand alive casino royale such as Grand Baccarat, Grand Blackjack and Grand Roulette for our

VIP member. And of course as a believer of Depoxito you can enjoy

all the games that we provide to you, every you craving to

get is just visit our site depoxito and register it without help takes in the works to

3 minutes and later youre welcome to deed any game that you want.

Be our VIP, living thing our VIP enthusiast of course

approved you the best further you can acquire from us all you infatuation to

be a VIP supporter is very easy. all you infatuation is just keep playing upon our

site, growth and discharge duty later a VIP gone the amount that our company had written, keep

playing and our customer relieve will approach you that you are promoted

to become a VIP advocate on our site.