User Request – Shepards Classification of Sediments

by Nicholas Hamilton 26/07/2014 | 2:25 18 Posted in Examples, Featured, User Request

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.

Tags:featured

Nick Hamilton

Author of the ggtern package. PhD candidate at the University of New South Wales in the School of Materials Science and Engineering. My field is in the simulations of metallic glasses via several advanced simulation techniques such as QMD and RMC.

Comments

Eduardo Richard 28/07/2014 | 12:13

Thanks for the help this is really useful to me

Reply
macroecology 31/07/2014 | 9:41

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?

Reply
Nicholas Hamilton 31/07/2014 | 10:02

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

Reply
Hy 16/09/2014 | 2:36

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

Reply

Nicholas Hamilton 16/09/2014 | 8:57

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

Reply

Rob & Pam 22/10/2014 | 7:39

Very cool, Nick. Congrats from Rob & Pam

Reply
Leandro 14/03/2016 | 10:50

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

Reply

Ugo 4/01/2018 | 12:01

Hi, I would appreciate it too 🙂
Thanks for the help.

Reply
Nicholas Hamilton 4/01/2018 | 12:13

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.

R

newPoints = data.frame(LNew=runif(10),TNew=runif(10),RNew=runif(10)) base + geom_point(data=newPoints,aes(LNew,TNew,RNew), size=5,colour='red')

1
2

newPoints = data.frame(LNew=runif(10),TNew=runif(10),RNew=runif(10))
base + geom_point(data=newPoints,aes(LNew,TNew,RNew), size=5,colour='red')

Reply

Ugo 8/01/2018 | 5:58

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?

Reply
Ugo 8/01/2018 | 6:08

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 :

R

<span class="GNKRCKGCMRB ace_keyword">base + geom_point(data=F,aes(P,O,N)) + stat_density_tern(aes(fill=..level.., alpha=..level..),...</span>

1

<span class="GNKRCKGCMRB ace_keyword">base + geom_point(data=F,aes(P,O,N)) + stat_density_tern(aes(fill=..level.., alpha=..level..),...</span>

Reply

Degang Wu 23/07/2016 | 7:59

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.

Reply

Nicholas Hamilton 4/08/2016 | 12:59

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

Reply

Kevin 28/02/2018 | 4:48

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)

Reply
Ugo 9/03/2018 | 8:21

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.

Reply
Ugo 9/03/2018 | 8:23

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

Reply
David Scheer 28/10/2018 | 10:49

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

Reply
Depoxitovip 8/12/2020 | 7:37

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.

Reply