Plotting XRD (X-Ray Diffraction) Data
Premise:
Some examples on how to prepare and present data collected from an XRD analysis. The clay fraction from seven horizons was analyzed by XRD, using the five common treatments: potassium saturation (K), potassium saturation heated to 350 Deg C (K 350), potassium saturation heated to 550 Deg C (K 550), magnesium saturation (Mg), and magnesium + glycerin saturation (Mg+GLY). These data files have been attached, and can be found near the bottom of the page.
Plotting the entire data set with lattice graphics:
## load libs
require(lattice)
require(reshape)
## read the composite data in as a table
## format is 2theta,MG,MG+GLY,K,K350,K550
h1 <- read.csv("tioga1_0-8.csv", header=FALSE)
h2 <- read.csv("tioga1_8-15.csv", header=FALSE)
h3 <- read.csv("tioga1_15-35.csv", header=FALSE)
h4 <- read.csv("tioga1_35-65.csv", header=FALSE)
h5 <- read.csv("tioga1_65-90.csv", header=FALSE)
h6 <- read.csv("tioga1_90-120.csv", header=FALSE)
h7 <- read.csv("tioga1_120-150.csv", header=FALSE)
## load some common d-spacings:
d_spacings <- c(0.33,0.358,0.434,0.482,0.717,1,1.2,1.4,1.8)
d_spacing_labels <- c(".33", ".36", ".43", ".48", ".7","1.0","1.2","1.4","1.8")
## combine horizons, and
xrd <- make.groups(h1, h2, h3, h4, h5, h6, h7)
names(xrd) <- c('x', 'MG', 'MG+GLY', 'K', 'K 350', 'K 550', 'horizon')
## convert data into long format
xrd.long <- melt(data=xrd, id.var=c('x', 'horizon'), measure.var=c('K','K 350', 'K 550', 'MG', 'MG+GLY'), variable_name='treatment')
## set a better ordering of the treatments
xrd.long$treatment <- ordered(xrd.long$treatment, c('MG', 'MG+GLY', 'K', 'K 350', 'K 550'))
## change the strip background colors
## trellis.par.set(list(strip.background = list(col = grey(c(0.9,0.8)) )))
## plot the data along with some common d-spacings:
xyplot(value ~ x | treatment + horizon , data=xrd.long, as.table=TRUE, ylim=c(0,500), xlab='Deg 2Theta', ylab='Counts', panel=function(x, y, ...) {panel.abline(v=(asin(0.154/(2*d_spacings)) * 180/pi * 2), col=grey(0.9)) ; panel.xyplot(x, y, ..., type='l', col='black')} )
## another approach: labels on the side
xyplot(value ~ x | horizon + treatment , data=xrd.long, as.table=TRUE, ylim=c(0,500), xlab='Deg 2Theta', ylab='Counts', panel=function(x, y, ...) {panel.abline(v=(asin(0.154/(2*d_spacings)) * 180/pi * 2), col=grey(0.9)) ; panel.xyplot(x, y, ..., type='l', col='black')}, strip.left=TRUE, strip=FALSE)
require(lattice)
require(reshape)
## read the composite data in as a table
## format is 2theta,MG,MG+GLY,K,K350,K550
h1 <- read.csv("tioga1_0-8.csv", header=FALSE)
h2 <- read.csv("tioga1_8-15.csv", header=FALSE)
h3 <- read.csv("tioga1_15-35.csv", header=FALSE)
h4 <- read.csv("tioga1_35-65.csv", header=FALSE)
h5 <- read.csv("tioga1_65-90.csv", header=FALSE)
h6 <- read.csv("tioga1_90-120.csv", header=FALSE)
h7 <- read.csv("tioga1_120-150.csv", header=FALSE)
## load some common d-spacings:
d_spacings <- c(0.33,0.358,0.434,0.482,0.717,1,1.2,1.4,1.8)
d_spacing_labels <- c(".33", ".36", ".43", ".48", ".7","1.0","1.2","1.4","1.8")
## combine horizons, and
xrd <- make.groups(h1, h2, h3, h4, h5, h6, h7)
names(xrd) <- c('x', 'MG', 'MG+GLY', 'K', 'K 350', 'K 550', 'horizon')
## convert data into long format
xrd.long <- melt(data=xrd, id.var=c('x', 'horizon'), measure.var=c('K','K 350', 'K 550', 'MG', 'MG+GLY'), variable_name='treatment')
## set a better ordering of the treatments
xrd.long$treatment <- ordered(xrd.long$treatment, c('MG', 'MG+GLY', 'K', 'K 350', 'K 550'))
## change the strip background colors
## trellis.par.set(list(strip.background = list(col = grey(c(0.9,0.8)) )))
## plot the data along with some common d-spacings:
xyplot(value ~ x | treatment + horizon , data=xrd.long, as.table=TRUE, ylim=c(0,500), xlab='Deg 2Theta', ylab='Counts', panel=function(x, y, ...) {panel.abline(v=(asin(0.154/(2*d_spacings)) * 180/pi * 2), col=grey(0.9)) ; panel.xyplot(x, y, ..., type='l', col='black')} )
## another approach: labels on the side
xyplot(value ~ x | horizon + treatment , data=xrd.long, as.table=TRUE, ylim=c(0,500), xlab='Deg 2Theta', ylab='Counts', panel=function(x, y, ...) {panel.abline(v=(asin(0.154/(2*d_spacings)) * 180/pi * 2), col=grey(0.9)) ; panel.xyplot(x, y, ..., type='l', col='black')}, strip.left=TRUE, strip=FALSE)
Example XRD plot with lattice graphics: 7 horizons and 5 treatments
Find peaks in an XRD dataset
Locating relevant peaks in an X-ray diffractogram is an important step in identifying phyllosilicate minerals in soils. An automated approach to finding peaks in any dataset was presented by Martin Maechler, contributed to the R-Help mailing list Nov 25, 2005. paste these functions into an R session to use them
peaks <- function(series, span = 3, do.pad = TRUE) {
if((span <- as.integer(span)) %% 2 != 1) stop("'span' must be odd")
s1 <- 1:1 + (s <- span %/% 2)
if(span == 1) return(rep.int(TRUE, length(series)))
z <- embed(series, span)
v <- apply(z[,s1] > z[, -s1, drop=FALSE], 1, all)
if(do.pad) {
pad <- rep.int(FALSE, s)
c(pad, v, pad)
} else v
}
peaksign <- function(series, span = 3, do.pad = TRUE)
{
## Purpose: return (-1 / 0 / 1) if series[i] is ( trough / "normal" / peak )
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 25 Nov 2005
if((span <- as.integer(span)) %% 2 != 1 || span == 1)
stop("'span' must be odd and >= 3")
s1 <- 1:1 + (s <- span %/% 2)
z <- embed(series, span)
d <- z[,s1] - z[, -s1, drop=FALSE]
ans <- rep.int(0:0, nrow(d))
ans[apply(d > 0, 1, all)] <- as.integer(1)
ans[apply(d < 0, 1, all)] <- as.integer(-1)
if(do.pad) {
pad <- rep.int(0:0, s)
c(pad, ans, pad)
} else ans
}
check.pks <- function(y, span = 3)
stopifnot(identical(peaks( y, span), peaksign(y, span) == 1),
identical(peaks(-y, span), peaksign(y, span) == -1))
for(y in list(1:10, rep(1,10), c(11,2,2,3,4,4,6,6,6))) {
for(sp in c(3,5,7))
check.pks(y, span = sp)
stopifnot(peaksign(y) == 0)
}
if((span <- as.integer(span)) %% 2 != 1) stop("'span' must be odd")
s1 <- 1:1 + (s <- span %/% 2)
if(span == 1) return(rep.int(TRUE, length(series)))
z <- embed(series, span)
v <- apply(z[,s1] > z[, -s1, drop=FALSE], 1, all)
if(do.pad) {
pad <- rep.int(FALSE, s)
c(pad, v, pad)
} else v
}
peaksign <- function(series, span = 3, do.pad = TRUE)
{
## Purpose: return (-1 / 0 / 1) if series[i] is ( trough / "normal" / peak )
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 25 Nov 2005
if((span <- as.integer(span)) %% 2 != 1 || span == 1)
stop("'span' must be odd and >= 3")
s1 <- 1:1 + (s <- span %/% 2)
z <- embed(series, span)
d <- z[,s1] - z[, -s1, drop=FALSE]
ans <- rep.int(0:0, nrow(d))
ans[apply(d > 0, 1, all)] <- as.integer(1)
ans[apply(d < 0, 1, all)] <- as.integer(-1)
if(do.pad) {
pad <- rep.int(0:0, s)
c(pad, ans, pad)
} else ans
}
check.pks <- function(y, span = 3)
stopifnot(identical(peaks( y, span), peaksign(y, span) == 1),
identical(peaks(-y, span), peaksign(y, span) == -1))
for(y in list(1:10, rep(1,10), c(11,2,2,3,4,4,6,6,6))) {
for(sp in c(3,5,7))
check.pks(y, span = sp)
stopifnot(peaksign(y) == 0)
}
Commands to find and plot the peaks, based on suggestions by peaks() function author.
## load some sample data
d <- read.csv("tioga1_35-65.csv", header=FALSE)
## name the columns
names(d) <- c('x', 'MG', 'MG+GLY', 'K', 'K 350', 'K 550')
## locate peaks in the 'K' signal
## the second argument is the "sensitivity" of the peak finding algorithm
d.peaks <- peaks(d$K, 35)
## save a vector of the positions in the K signal where the peaks were identified
peak_idx <- which(d.peaks)
## simple plot of raw K signal
plot(K ~ x, data=d, type="l", cex=.25, main="Tioga1 35-65cm\nK Treatment", xlab="Deg. 2 Theta", ylab="Intensity", ylim=c(0,max(d$K) + 50))
## add peaks: note that we are sub0setting the original data by the peak location index
points(K ~ x, data=d[peak_idx, ], col = 2, cex = 1.5)
## compute peak d-spacings
peak_d_spacings <- signif( (0.154/ (sin(d$x[peak_idx] * pi/180))), 3)
## annotate d-spacings / or peak index
## text(d$x[peak_idx], d$K[peak_idx] + 20, peak_d_spacings, col='blue', cex=0.75 )
text(d$x[peak_idx], d$K[peak_idx] + 20, 1:length(peak_d_spacings), col='blue', cex=0.75 )
## print a simple table of d-spacings by index
data.frame(peak_d_spacings)
peak_d_spacings
1 1.370
2 1.040
3 0.751
4 0.638
5 0.563
6 0.488
7 0.452
8 0.435
9 0.367
10 0.347
d <- read.csv("tioga1_35-65.csv", header=FALSE)
## name the columns
names(d) <- c('x', 'MG', 'MG+GLY', 'K', 'K 350', 'K 550')
## locate peaks in the 'K' signal
## the second argument is the "sensitivity" of the peak finding algorithm
d.peaks <- peaks(d$K, 35)
## save a vector of the positions in the K signal where the peaks were identified
peak_idx <- which(d.peaks)
## simple plot of raw K signal
plot(K ~ x, data=d, type="l", cex=.25, main="Tioga1 35-65cm\nK Treatment", xlab="Deg. 2 Theta", ylab="Intensity", ylim=c(0,max(d$K) + 50))
## add peaks: note that we are sub0setting the original data by the peak location index
points(K ~ x, data=d[peak_idx, ], col = 2, cex = 1.5)
## compute peak d-spacings
peak_d_spacings <- signif( (0.154/ (sin(d$x[peak_idx] * pi/180))), 3)
## annotate d-spacings / or peak index
## text(d$x[peak_idx], d$K[peak_idx] + 20, peak_d_spacings, col='blue', cex=0.75 )
text(d$x[peak_idx], d$K[peak_idx] + 20, 1:length(peak_d_spacings), col='blue', cex=0.75 )
## print a simple table of d-spacings by index
data.frame(peak_d_spacings)
peak_d_spacings
1 1.370
2 1.040
3 0.751
4 0.638
5 0.563
6 0.488
7 0.452
8 0.435
9 0.367
10 0.347
Peaks found with two different tolerance settings.
Some ideas on annotating common d-spacings
## note that we need to define this function before we can use it
## put all of the plotting commands into a wrapper function:
plot_xrd <- function(d)
{
## plot the difractograms offset by 200
plot(MG ~ x, data=d, type="l", cex=.25, main="Tioga1 35-65cm", xlab="2 Theta", ylab="Intensity", xlim=c(2,32), ylim=c(0,1600), xaxt='none')
## add the other treatments
lines(MG_GLY + 200 ~ x, data=d)
lines(K + 400 ~ x, data=d)
lines(K_350 + 600 ~ x, data=d)
lines(K_550 + 800 ~ x, data=d)
## label the lines
text(31.5, c(30,230,430,630,830), c("MG","MG+GLY","K25","K350","K550"))
## plot the zero-line on each graph:
abline(h=c(0,200,400,600,800), lty=2, col="gray")
## plot the established boundaries to these common spacings
abline(v=(asin(0.154/(2*c(0.71,0.73,0.72,0.75,0.99,1.01,1.24,1.28,1.4,1.5,1.77,1.8))) * 180/pi * 2), col="green", lty=2)
## plot some common d-spacings
# abline(v=(asin(0.154/(2*c(.715,.73,1,1.2,1.4,1.8))) * 180/pi * 2), col="blue")
## annotate the lines, recall that d-spacing is in reverse order with respect to Two_theta
text((asin(0.154/(2*c(.715,.73,1,1.2,1.4,1.8))) * 180/pi * 2), c(1300,1400,1400,1400,1400,1400), c(".715",".73","1.0","1.2","1.4","1.8"), col=1, cex=1)
## add the axis
axis(1, 1:30)
}
##
##
##
## load some sample data
d <- read.csv("tioga1_35-65.csv", header=FALSE)
## name the columns
names(d) <- c('x', 'MG', 'MG_GLY', 'K', 'K_350', 'K_550')
## run the wrapper function to plot the Tioga1 0-8cm data:
plot_xrd(d)
## put all of the plotting commands into a wrapper function:
plot_xrd <- function(d)
{
## plot the difractograms offset by 200
plot(MG ~ x, data=d, type="l", cex=.25, main="Tioga1 35-65cm", xlab="2 Theta", ylab="Intensity", xlim=c(2,32), ylim=c(0,1600), xaxt='none')
## add the other treatments
lines(MG_GLY + 200 ~ x, data=d)
lines(K + 400 ~ x, data=d)
lines(K_350 + 600 ~ x, data=d)
lines(K_550 + 800 ~ x, data=d)
## label the lines
text(31.5, c(30,230,430,630,830), c("MG","MG+GLY","K25","K350","K550"))
## plot the zero-line on each graph:
abline(h=c(0,200,400,600,800), lty=2, col="gray")
## plot the established boundaries to these common spacings
abline(v=(asin(0.154/(2*c(0.71,0.73,0.72,0.75,0.99,1.01,1.24,1.28,1.4,1.5,1.77,1.8))) * 180/pi * 2), col="green", lty=2)
## plot some common d-spacings
# abline(v=(asin(0.154/(2*c(.715,.73,1,1.2,1.4,1.8))) * 180/pi * 2), col="blue")
## annotate the lines, recall that d-spacing is in reverse order with respect to Two_theta
text((asin(0.154/(2*c(.715,.73,1,1.2,1.4,1.8))) * 180/pi * 2), c(1300,1400,1400,1400,1400,1400), c(".715",".73","1.0","1.2","1.4","1.8"), col=1, cex=1)
## add the axis
axis(1, 1:30)
}
##
##
##
## load some sample data
d <- read.csv("tioga1_35-65.csv", header=FALSE)
## name the columns
names(d) <- c('x', 'MG', 'MG_GLY', 'K', 'K_350', 'K_550')
## run the wrapper function to plot the Tioga1 0-8cm data:
plot_xrd(d)
Example XRD plot 2: illustrating common d-spacings
Two-page display of XRD data for an entire soil profile
Sample 1
Sample 2