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

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.