The goal of xrftools is to provide tools to read, plot, and interpret X-Ray fluorescence spectra.
You can install xrftools from github with:
# install.packages("remotes")
remotes::install_github("paleolimbot/xrftools")Read in a Panalytical XRF spectrum and plot it.
library(tidyverse)
library(xrftools)
#>
#> Attaching package: 'xrftools'
#> The following object is masked from 'package:stats':
#>
#> filter
pan_example_dir <- system.file("spectra_files/Panalytical", package = "xrftools")
pan_files <- list.files(pan_example_dir, ".mp2", full.names = TRUE)
specs <- read_xrf_panalytical(pan_files)
specs %>%
xrf_despectra() %>%
unnest(.spectra) %>%
ggplot(aes(x = energy_kev, y = cps, col = SampleIdent)) +
geom_line() +
facet_wrap(vars(ConditionSet), scales = "free_y")
The xrf package can use several existing methods for estimating “background” or “baseline” values. The most useful of these for XRF spectra is the Sensitive Nonlinear Iterative Peak (SNIP) method, implemented in the Peaks package.
specs %>%
slice(3) %>%
xrf_add_baseline_snip(iterations = 20) %>%
xrf_despectra() %>%
unnest() %>%
filter(energy_kev <= 15) %>%
ggplot(aes(x = energy_kev)) +
geom_line(aes(y = cps, col = "raw")) +
geom_line(aes(y = baseline, col = "baseline")) +
geom_line(aes(y = cps - baseline, col = "cps - baseline"))
#> Warning: `cols` is now required when using unnest().
#> Please use `cols = c(.spectra)`
specs %>%
slice(3) %>%
xrf_add_baseline_snip(iterations = 20) %>%
xrf_add_smooth_filter(filter = xrf_filter_gaussian(alpha = 2.5), .iter = 5) %>%
xrf_despectra() %>%
unnest() %>%
filter(energy_kev <= 15) %>%
ggplot(aes(x = energy_kev)) +
geom_line(aes(y = cps, col = "raw"), alpha = 0.3) +
geom_line(aes(y = smooth - baseline, col = "smooth"))
#> Warning: `cols` is now required when using unnest().
#> Please use `cols = c(.spectra)`
deconvoluted <- specs %>%
filter(ConditionSet %in% c("Omnian", "Omnian2")) %>%
xrf_add_smooth_filter(filter = xrf_filter_gaussian(width = 5), .iter = 20) %>%
xrf_add_baseline_snip(.values = .spectra$smooth, iterations = 20) %>%
xrf_add_deconvolution_gls(
.spectra$energy_kev,
.spectra$smooth - .spectra$baseline,
energy_max_kev = kV * 0.75, peaks = xrf_energies("major")
)
certified_vals <- system.file("spectra_files/oreas_concentrations.csv", package = "xrftools") %>%
read_csv(col_types = cols(standard = col_character(), value = col_double(), .default = col_guess())) %>%
filter(method == "4-Acid Digestion") %>%
select(SampleIdent = standard, element, certified_value = value)
deconv <- deconvoluted %>%
unnest(.deconvolution_peaks) %>%
select(SampleIdent, ConditionSet, kV, element, height, peak_area, peak_area_se)
deconv %>%
left_join(certified_vals, by = c("SampleIdent", "element")) %>%
ggplot(aes(x = certified_value, y = peak_area, col = SampleIdent, shape = ConditionSet)) +
geom_errorbar(aes(ymin = peak_area - peak_area_se, ymax = peak_area + peak_area_se)) +
geom_point() +
facet_wrap(~element, scales = "free") +
theme_bw(10)
#> Warning: Removed 22 rows containing missing values (geom_point).
deconv_element <- deconvoluted %>%
unnest(.deconvolution_components)
ggplot() +
geom_line(
aes(x = energy_kev, y = response),
data = deconvoluted %>% unnest(.deconvolution_response),
size = 0.2
) +
geom_area(
aes(x = energy_kev, y = response_fit, fill = element),
data = deconvoluted %>% unnest(.deconvolution_components),
alpha = 0.5
) +
facet_wrap(~ConditionSet + SampleIdent, scales = "free") +
scale_y_sqrt() +
theme_bw(10)
#> Warning in self$trans$transform(x): NaNs produced
#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 4504 rows containing missing values (position_stack).
spec <- specs %>%
slice(7) %>%
xrf_add_baseline_snip(iterations = 20) %>%
xrf_add_smooth_filter(filter = xrf_filter_gaussian(alpha = 2.5), .iter = 5) %>%
pull(.spectra) %>%
first()
sigma_index <- 6
energy_kev <- spec$energy_kev
values <- spec$fit - spec$background
energy_res <- mean(diff(energy_kev))
search <- Peaks::SpectrumSearch(values, threshold = 0.01, sigma = sigma_index)
g <- function(x, mu = 0, sigma = 1, height = 1) height * exp(-0.5 * ((x - mu) / sigma) ^ 2)
tbl <- tibble::tibble(
peak_index = search$pos,
peak_energy_kev = energy_kev[peak_index],
peak_sigma = sigma_index * energy_res,
peak_height = values[peak_index],
peak_area = peak_height * peak_sigma * sqrt(2 * pi)
)
peak_response = pmap(list(mu = tbl$peak_energy_kev, sigma = tbl$peak_sigma, height = tbl$peak_height), g, energy_kev)
spec$deconv <- do.call(cbind, peak_response) %>%
rowSums()
spec$deconv2 <- search$y
ggplot(spec, aes(energy_kev)) +
geom_line(aes(y = fit - background, col = "original")) +
geom_line(aes(y = deconv, col = "deconv")) +
xlim(0, 10) +
stat_xrf_peaks(aes(y = fit - background), epsilon = 10)oreas22d <- specs %>%
filter(SampleIdent == "oreas 22d") %>%
unnest(.spectra) %>%
mutate(cps = counts / LiveTime)
xrf_en <- x_ray_xrf_energies %>%
crossing(tibble(ConditionSet = unique(oreas22d$ConditionSet))) %>%
group_by(ConditionSet) %>%
mutate(
data = list(oreas22d[oreas22d$ConditionSet == ConditionSet[1],]),
counts = approx(data[[1]]$energy_kev, data[[1]]$counts, energy_kev)$y,
fit = approx(data[[1]]$energy_kev, data[[1]]$fit, energy_kev)$y,
background = approx(data[[1]]$energy_kev, data[[1]]$background, energy_kev)$y,
cps = approx(data[[1]]$energy_kev, data[[1]]$cps, energy_kev)$y
) %>%
select(-data)
library(plotly)
plot_ly() %>%
add_lines(x = ~energy_kev, y = ~cps, color = ~ConditionSet, hoverinfo = "none",
data = oreas22d) %>%
add_markers(x = ~energy_kev, y = ~cps, text = ~element, color = ~ConditionSet, data = xrf_en)