soils plot basic
Matt Harris
Wednesday, February 25, 2015
Making soil profiles V1
Using the 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 wnat 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("C:/Users/matthew_harris/Dropbox/R/Soil_profiles")
# 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 + 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 addedLet’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.96Label
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 - #94764FConversion
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 variablesclass(dat)## [1] "SoilProfileCollection"
## attr(,"package")
## [1] "aqp"summary(dat)## Length Class Mode
## 4 SoilProfileCollection S4Plotting
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')
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
## 2Adding 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 #482A22FFArtifact 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 <- 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")
# 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