Making Soil Property vs. Depth Plots
Example with randomly generated data
Generate some data
## generate some profile depths: 0 - 150, in 10 cm increments
depth <- seq(0,150, by=10)
## generate some property: random numbers in this case
prop <- rnorm(n=length(depth), mean=15, sd=2)
## since the 0 is not a depth, and we would like the graph to start from 0
## make the first property row (associated with depth 0) the same as the second
## property row
prop[1] <- prop[2]
## combine into a table: data read in from a spread sheet would already be in this format
soil <- data.frame(depth=depth, prop=prop)
depth <- seq(0,150, by=10)
## generate some property: random numbers in this case
prop <- rnorm(n=length(depth), mean=15, sd=2)
## since the 0 is not a depth, and we would like the graph to start from 0
## make the first property row (associated with depth 0) the same as the second
## property row
prop[1] <- prop[2]
## combine into a table: data read in from a spread sheet would already be in this format
soil <- data.frame(depth=depth, prop=prop)
The dataframe 'soil' looks like this:
depth prop 1 0 13.80257 ** note that these are the same 2 12 13.80257 ** note that these are the same 3 24 18.40298 4 36 13.37446 5 48 13.27973 6 60 14.65288 7 72 16.07339 8 84 15.97451 9 96 16.29970 10 108 16.32155 11 120 14.63699 12 132 13.26486 13 144 13.81730
Plot the data:
## note the reversal of the y-axis with ylim=c(150,0)
plot(depth ~ prop, data=soil, ylim=c(150,0), type='s', ylab='Depth', xlab='Property', main='Property vs. Depth Plot')
plot(depth ~ prop, data=soil, ylim=c(150,0), type='s', ylab='Depth', xlab='Property', main='Property vs. Depth Plot')
Additional Example Using Lattice Graphics
Examples with Some Real Data
Notes:
- See attached files at bottom of page
- Helper function could use some generalization. Until then, your data will need to have the columns:
- pedon_id
- top
- bottom
- These examples require a recent version of R and Lattice (>= 2.5.1)
Helper Function (copy this into your R session first)
## function to add a repeat top horizon
## for correct step-plot
## assumes that there are columns named pedon_id, bottom, top
profile_fix <- function(d)
{
## init some vars
p <- levels(d$pedon_id)
idx <- array()
i <- 1
## loop over pedon ids
for(p.id in p)
{
## extract one at a time
p.row <- subset(d, subset=pedon_id == p.id)
## make a list of the positions where bottom horizons occur
idx[i] <- which(d$top == min(p.row$top) & d$pedon_id == p.id)
## increment counter
i <- i+1
}
## extract out bottom horizons
d.temp <- d[idx,]
## set the top of these to the bottom boundary
d.temp$bottom <- d.temp$top
## add duplicate bottom horizon records, with top set to bottom
d.new <- rbind(d, d.temp)
return(d.new)
}
## for correct step-plot
## assumes that there are columns named pedon_id, bottom, top
profile_fix <- function(d)
{
## init some vars
p <- levels(d$pedon_id)
idx <- array()
i <- 1
## loop over pedon ids
for(p.id in p)
{
## extract one at a time
p.row <- subset(d, subset=pedon_id == p.id)
## make a list of the positions where bottom horizons occur
idx[i] <- which(d$top == min(p.row$top) & d$pedon_id == p.id)
## increment counter
i <- i+1
}
## extract out bottom horizons
d.temp <- d[idx,]
## set the top of these to the bottom boundary
d.temp$bottom <- d.temp$top
## add duplicate bottom horizon records, with top set to bottom
d.new <- rbind(d, d.temp)
return(d.new)
}
Load Data and Packages
## load libs
library(lattice)
## read in the first example
x <- read.csv('psa.csv')
## convert pedon_id to a factor:
x <- transform(x, pedon_id = factor(pedon_id))
## add extra top horizon
x.new <- profile_fix(x)
##
##
## read in the second example
y <- read.csv('example_data.csv')
## add the extra top horizon
y.new <- profile_fix(y)
library(lattice)
## read in the first example
x <- read.csv('psa.csv')
## convert pedon_id to a factor:
x <- transform(x, pedon_id = factor(pedon_id))
## add extra top horizon
x.new <- profile_fix(x)
##
##
## read in the second example
y <- read.csv('example_data.csv')
## add the extra top horizon
y.new <- profile_fix(y)
Example 1
## plot using step function
## note special syntax: horizontal=TRUE
xyplot(bottom ~ sand + silt + clay | pedon_id, horizontal=TRUE,
data=x.new, ylim=c(160,-5), type='s', auto.key=TRUE,
col=c('Orange', 'RoyalBlue', 'Dark Green'), lty=c(2,2,1), lwd=c(1,1,2),
ylab='Depth (cm)', xlab='Percent Sand, Silt, Clay',
key=list(
lines=list(col=c('Orange', 'RoyalBlue', 'Dark Green'), lwd=c(1,1,2), lty=c(2,2,1)
),
text=list(
c('Sand', 'Silt', 'Clay')
)
)
)
## note special syntax: horizontal=TRUE
xyplot(bottom ~ sand + silt + clay | pedon_id, horizontal=TRUE,
data=x.new, ylim=c(160,-5), type='s', auto.key=TRUE,
col=c('Orange', 'RoyalBlue', 'Dark Green'), lty=c(2,2,1), lwd=c(1,1,2),
ylab='Depth (cm)', xlab='Percent Sand, Silt, Clay',
key=list(
lines=list(col=c('Orange', 'RoyalBlue', 'Dark Green'), lwd=c(1,1,2), lty=c(2,2,1)
),
text=list(
c('Sand', 'Silt', 'Clay')
)
)
)
Depth Profile Example 1: sand, silt, and clay vs. depth for three pedons
Example 2
## plot with step function
xyplot(bottom ~ field_pct_clay | pedon_id, horizontal=TRUE,
data=y.new, ylim=c(250,-5), type='s', as.table=TRUE,
ylab='Depth (cm)', xlab='Percent Clay', lwd=2, col='black'
)
xyplot(bottom ~ field_pct_clay | pedon_id, horizontal=TRUE,
data=y.new, ylim=c(250,-5), type='s', as.table=TRUE,
ylab='Depth (cm)', xlab='Percent Clay', lwd=2, col='black'
)
Depth Profile Example 2: clay vs. depth for 9 pedons