Plotting Shovel Test Unit Profiles in R with the aqp package

showing the inclusion of densities.
Examples of STU profiles

The code in the following post was first generated to find a quick way to create a large number of Shovel Test Unit (STU) soil profiles with some consistency and a few added features. This methodology uses the aqp package in R as the base function to generate the profiles and data in the form of a *.csv file. This package was originally developed to work with geology or soil science data, but we adapted it a bit for our purpose. Subsequent to developing this initial code, others in the office further developed the function to fit their work flow and profiles. The basics are here and adapting it is pretty straight forward if you dive into the guts of the amp package and read the docs.

Features shown here include converting Munsell colors to hex, artifact densities, exporting as PDF and a few other tweaks. Hopefully, you will find it useful!

aqp package http://aqp.r-forge.r-project.org/

first thing is to load the packages. The first time you do this you will need to install the packages. After that, you can just import them.

# this is the main package for profiles
install.packages("aqp") 
# these may come in handy
install.packages("RColorBrewer")
install.packages("latticeExtra")
install.packages("plyr")
install.packages("reshape")

Once the packages are installed, we need to import tham into the R environment as such

require(aqp)
require(RColorBrewer)
require(latticeExtra)
require(plyr)
require(reshape)

Now these packages are imported. That means we can use all the functions and data the come with these packages. Google searches will return documentation on these packages and tell you want functions they contain and some examples.

Data

time to load the data

# this sets the working directory to a folder that contains your data and will
# contain your output.  Setting this is not mandatory, but makes it easier to get
# and put data with less typing [Notice that the slashes in the URL as forward /,
# not the standard back slash \]
setwd("[Your File Location or Working Directory]")

# in this working directory, I made a folder called data to hold the 
# CSV of STU profiles
# I read this using the read.csv function and assign it to the variable "dat"
dat <- read.csv(paste(c(getwd(), "/data/Soil_Profiles.csv"), collapse=''), 
                stringsAsFactors = FALSE)[,-1]

# that "paste" fucntion concatenates the working directory and the 
# file location and name
# set "stringsAsFactors = FALSE" tells the functin to leave text fields as text 
# in the data
# the "[,-1]" at the end drops the index column that is automatically added

Let’s see out data

# this will print the whole thing.  use head(dat,#) where # is replaced with the
# number of rows you want to print
print(dat)
##         id top bottom name    p1 soil_color
## 1  Pole 13   0     19   Ap 6.963    #5E4423
## 2  Pole 13  19     55    B 5.951    #815A21
## 3  Pole 14   0     24   Ap 5.245    #3D2F21
## 4  Pole 14  24     45   Ab 5.069    #5E4423
## 5  Pole 14  45     50    E 5.069    #BEAD97
## 6  Pole 14  50     60    B 5.338    #B88E50
## 7  Pole 22   0     19   Ap 6.454    #7A5C36
## 8  Pole 22  19     30   Ab 5.844    #554636
## 9  Pole 22  30     38    B 5.437    #94764F
## 10 Pole 22  38     50    C 5.181    #9D7337
## 11 Pole 37   0     31   Ap 5.019    #7A5C36
## 12 Pole 37  31     63  Ap1 5.283    #554636
## 13 Pole 37  63     90  Ap2 6.440    #5E4423
## 14 Pole 37  90    117    B 6.670    #94764F
# the "str" function reports the structure of the data
str(dat)
## 'data.frame':    14 obs. of  6 variables:
##  $ id        : chr  "Pole 13" "Pole 13" "Pole 14" "Pole 14" ...
##  $ top       : int  0 19 0 24 45 50 0 19 30 38 ...
##  $ bottom    : int  19 55 24 45 50 60 19 30 38 50 ...
##  $ name      : chr  "Ap" "B" "Ap" "Ab" ...
##  $ p1        : num  6.96 5.95 5.24 5.07 5.07 ...
##  $ soil_color: chr  "#5E4423" "#815A21" "#3D2F21" "#5E4423" ...
# class tells us what class the object "dat" is
class(dat)
## [1] "data.frame"
# summary gives us a look into the range of the data
summary(dat)
##       id                 top            bottom          name          
##  Length:14          Min.   : 0.00   Min.   : 19.0   Length:14         
##  Class :character   1st Qu.: 4.75   1st Qu.: 30.2   Class :character  
##  Mode  :character   Median :27.00   Median : 47.5   Mode  :character  
##                     Mean   :29.21   Mean   : 49.4                     
##                     3rd Qu.:43.25   3rd Qu.: 58.8                     
##                     Max.   :90.00   Max.   :117.0                     
##        p1        soil_color       
##  Min.   :5.02   Length:14         
##  1st Qu.:5.20   Class :character  
##  Median :5.39   Mode  :character  
##  Mean   :5.71                     
##  3rd Qu.:6.32                     
##  Max.   :6.96

Label

Add a label with hoirzon and munsell color to data.frame

dat <- within(dat,  label <- paste(name, soil_color, sep=" - "))

print(dat)
##         id top bottom name    p1 soil_color         label
## 1  Pole 13   0     19   Ap 6.963    #5E4423  Ap - #5E4423
## 2  Pole 13  19     55    B 5.951    #815A21   B - #815A21
## 3  Pole 14   0     24   Ap 5.245    #3D2F21  Ap - #3D2F21
## 4  Pole 14  24     45   Ab 5.069    #5E4423  Ab - #5E4423
## 5  Pole 14  45     50    E 5.069    #BEAD97   E - #BEAD97
## 6  Pole 14  50     60    B 5.338    #B88E50   B - #B88E50
## 7  Pole 22   0     19   Ap 6.454    #7A5C36  Ap - #7A5C36
## 8  Pole 22  19     30   Ab 5.844    #554636  Ab - #554636
## 9  Pole 22  30     38    B 5.437    #94764F   B - #94764F
## 10 Pole 22  38     50    C 5.181    #9D7337   C - #9D7337
## 11 Pole 37   0     31   Ap 5.019    #7A5C36  Ap - #7A5C36
## 12 Pole 37  31     63  Ap1 5.283    #554636 Ap1 - #554636
## 13 Pole 37  63     90  Ap2 6.440    #5E4423 Ap2 - #5E4423
## 14 Pole 37  90    117    B 6.670    #94764F   B - #94764F

Conversion

Now we need to put these data in the format that is used by the aqp package to do the cool stuff. This “depths” function takes the csv data and converts is into a object of class “SoilProfileCollection” for future use. By doing so, the aqp package conforms that data to a type it knows well and can use in the functions the package provides.

# the depth function uses the id field in the data to set the STU number
# it uses the top to set the top of the STU and bottom for the... well the bottom.
depths(dat) <- id ~ top + bottom

# now we see that the class of dat has been changed from data.frame to "SoilProfileCollection"
str(dat) # you can see that the original data is now within this object with other attributes
## Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##   ..@ idcol     : chr "id"
##   ..@ depthcols : chr [1:2] "top" "bottom"
##   ..@ metadata  :'data.frame':   1 obs. of  1 variable:
##   .. ..$ depth_units: chr "cm"
##   ..@ horizons  :'data.frame':   14 obs. of  7 variables:
##   .. ..$ id        : chr [1:14] "Pole 13" "Pole 13" "Pole 14" "Pole 14" ...
##   .. ..$ top       : int [1:14] 0 19 0 24 45 50 0 19 30 38 ...
##   .. ..$ bottom    : int [1:14] 19 55 24 45 50 60 19 30 38 50 ...
##   .. ..$ name      : chr [1:14] "Ap" "B" "Ap" "Ab" ...
##   .. ..$ p1        : num [1:14] 6.96 5.95 5.24 5.07 5.07 ...
##   .. ..$ soil_color: chr [1:14] "#5E4423" "#815A21" "#3D2F21" "#5E4423" ...
##   .. ..$ label     : chr [1:14] "Ap - #5E4423" "B - #815A21" "Ap - #3D2F21" "Ab - #5E4423" ...
##   ..@ site      :'data.frame':   4 obs. of  1 variable:
##   .. ..$ id: chr [1:4] "Pole 13" "Pole 14" "Pole 22" "Pole 37"
##   ..@ sp        :Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   .. .. ..@ coords     : num [1, 1] 0
##   .. .. ..@ bbox       : logi [1, 1] NA
##   .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
##   .. .. .. .. ..@ projargs: chr NA
##   ..@ diagnostic:'data.frame':   0 obs. of  0 variables
class(dat)
## [1] "SoilProfileCollection"
## attr(,"package")
## [1] "aqp"
summary(dat)
##                Length                 Class                  Mode 
##                     4 SoilProfileCollection                    S4

Plotting

Whats left to do is to plot! The plot function taks the “SoilProfileCollection” object called dat and uses the “name” field in the data to assign soil horizons… simple as that!

plot(dat, name = 'label')
Basic profile with correct munsell colors
Basic profile with correct munsell colors

In R studo this can be save to a PDF in the plot window, but in basic R, you can use the following

# this will save a pdf called "profiles" to your working directory 
# after you run "dev.off()"
pdf("profiles.pdf")
plot(dat, name = 'label')
dev.off()
## pdf 
##   2

Adding Munsell colors from field notation (e.g. “10YR3/2”)

# lets make up some example data.  Put the Hue, Value, and Chroma into 
# seperate colums
# this should be included in the original CSV probably
hue <- c("10YR", "10YR", "5YR", "2.5YR", "2.5YR")
valu <- c(3, 5, 4, 6, 2)
chroma <- c(2, 4, 4, 2, 3)
# We join these together into a data.frame called munsell 
# note "stringsAsFactors = FALSE"
munsell <- data.frame(hue, valu, chroma, stringsAsFactors = FALSE)
print(munsell)
##     hue valu chroma
## 1  10YR    3      2
## 2  10YR    5      4
## 3   5YR    4      4
## 4 2.5YR    6      2
## 5 2.5YR    2      3
# now use the "munsell2rgb" function to convert by giving it the three column 
# of hue, value, chroma.  This is added as a field called "color" to the munsell data.frame
munsell$color <- munsell2rgb(munsell$hue, munsell$valu, munsell$chroma)

#see!
print(munsell)
##     hue valu chroma     color
## 1  10YR    3      2 #554636FF
## 2  10YR    5      4 #947650FF
## 3   5YR    4      4 #805840FF
## 4 2.5YR    6      2 #A79086FF
## 5 2.5YR    2      3 #482A22FF

Artifact Counts, fractions, and depth of cultural material

Reloading the original csv data here to add some fake artifact counts to it

# reload data
dat <- read.csv(paste(c(getwd(), "/data/Soil_Profiles.csv"), collapse=''), 
                stringsAsFactors = FALSE)[,-1]
# add field for count of jasper flakes
dat$jasper <- sample(1:20, nrow(dat), replace=TRUE)
# a little function to compute the proportion of flakes in each horizon 
# (e.g 50% in Ap)
dat$jasperpcnt <- unlist(sapply(split(dat, f = dat$id), 
                  function(x) (x$jasper/sum(x$jasper)*100), simplify=TRUE))
# add a label field that has the horizon and jasper flake count
dat <- within(dat, jasperlabel    &lt;- paste(name, jasper, sep=" - "))

# see the data
print(dat)
##         id top bottom name    p1 soil_color jasper jasperpcnt jasperlabel
## 1  Pole 13   0     19   Ap 6.963    #5E4423      3     42.857      Ap - 3
## 2  Pole 13  19     55    B 5.951    #815A21      4     57.143       B - 4
## 3  Pole 14   0     24   Ap 5.245    #3D2F21     18     42.857     Ap - 18
## 4  Pole 14  24     45   Ab 5.069    #5E4423      5     11.905      Ab - 5
## 5  Pole 14  45     50    E 5.069    #BEAD97      7     16.667       E - 7
## 6  Pole 14  50     60    B 5.338    #B88E50     12     28.571      B - 12
## 7  Pole 22   0     19   Ap 6.454    #7A5C36      3      6.977      Ap - 3
## 8  Pole 22  19     30   Ab 5.844    #554636     20     46.512     Ab - 20
## 9  Pole 22  30     38    B 5.437    #94764F     13     30.233      B - 13
## 10 Pole 22  38     50    C 5.181    #9D7337      7     16.279       C - 7
## 11 Pole 37   0     31   Ap 5.019    #7A5C36     20     38.462     Ap - 20
## 12 Pole 37  31     63  Ap1 5.283    #554636      2      3.846     Ap1 - 2
## 13 Pole 37  63     90  Ap2 6.440    #5E4423     18     34.615    Ap2 - 18
## 14 Pole 37  90    117    B 6.670    #94764F     12     23.077      B - 12
# convert to apq object 
depths(dat) <- id ~ top + bottom

# plot with jasper quantity as color and lable with counts
plot(dat, name = "jasperlabel", color = "jasper")

# add overlay symbolizing the percent of jasper in each horizon
addVolumeFraction(dat, "jasperpcnt")
crazy soil colors and artifact densities
Color by artifact densities
# add in brackets to show depth of cultural material (fake in this example case)
fake.tops <- rep(0, times=length(dat))
fake.bottoms <- runif(n=length(dat), min=0, max=20)
par(mar=c(0,0,3,0)) # tighter figure margins
plot(dat, name = "jasperlabel", color = "soil_color")
addVolumeFraction(dat, "jasperpcnt")
addBracket(fake.tops, fake.bottoms, col='red') # add depth brackets

showing the inclusion of densities.
showing the inclusion of densities.
Advertisements
Plotting Shovel Test Unit Profiles in R with the aqp package

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s