User Request – Shepards Classification of Sediments

July 25, 2014
By

(This article was first published on ggtern: ternary diagrams in R, and kindly contributed to R-bloggers)

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 his blog: ggtern: ternary diagrams in R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.