User Request – Shepards Classification of Sediments

[This article was first published on ggtern: ternary diagrams in R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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).

plot

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")

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"))
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)
              )
            )
polygons$PointOrder <- 1:nrow(polygons) #IMPORTANT FOR CORRECT ORDERING.
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.

The post User Request – Shepards Classification of Sediments appeared first on ggtern: ternary diagrams in R.

To leave a comment for the author, please follow the link and comment on their blog: ggtern: ternary diagrams in R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)