ncss-tech / soilDB

soilDB: Simplified Access to National Cooperative Soil Survey Databases
http://ncss-tech.github.io/soilDB/
GNU General Public License v3.0
81 stars 20 forks source link

Can I use fetchKSSL to get trace elements, like Lead, Copper, Zinc, Mercury, etc. #118

Closed Kramdog closed 4 years ago

Kramdog commented 4 years ago

Hey github.
I'm trying to incorporate trace elements that the Lincoln lab will get for me and putt it into R Studio to be able to compare those results to XRF data. Specifically heavy metal data.

Thanks, Ben

dylanbeaudette commented 4 years ago

Thanks Ben, this is a good idea. There are several parts that require work:

brownag commented 4 years ago

Key attributes for "flattening" functions in each of the new tables

geochemical

prep_code, major_element_method, trace_element_method

glass

prep_code, analyzed_size_frac, glass_count_method

xray_thermal

prep_code, analyzed_size_frac, x_ray_method, differential_scanning_calorimeter_method, differential_thermal_analysis_method, thermal_gravimetric_method

brownag commented 4 years ago

Summary attributes for each of the new tables

geochemical

iron_oxide_total, aluminum_oxide_total, potassium_oxide_total

glass

total_grains_counted, resistant_minerals_total_mineral_soil, glass_count_mineral_interpretation

xray_thermal

clay_mineral_interpretation, coarse_silt_mineral_interpretation, fine_sand_mineral_interpretation, very_fine_sand_mineral_interpretation

brownag commented 4 years ago

@Kramdog

4bc6d32 - adds basic functionality for getting geochemical, optical and x-ray lab data; install latest version of soilDB off GitHub and you can go and look for your trace elements in the geochem table.

Now, you can do many things! I would be interested to see what you can pull out of the database in your area.

for instance:

# get geochemical and morphologic data for all KSSL pedons in MLRA 18
res <- fetchKSSL(mlra="18", returnGeochemicalData = TRUE, returnMorphologicData = TRUE)

# do some simple analyses on each of the new tables

# look at density plot of Fe2O3 %
plot(density(res$geochem$iron_oxide_total, na.rm=TRUE))

# get labsampnums for horizons with >90% resistant minerals
unique(res$optical$labsampnum[which(res$optical$resistant_minerals_total_mineral_soil > 90)])

# look at xrd peak heights for kaolinite
table(res$xrd_thermal$kk_kaolinite_x_ray)

Soon to come: functions to facilitate manipulating the individual tables and getting them into a format that is 1:1 with SPC horizon records.

brownag commented 4 years ago

@Kramdog I am starting to play around with workflows that use these new data tables. Here is a quick demo you might find useful. Needs latest aqp and soilDB.

# fetchKSSL heavy metals demo by mlra (all horizons with optical glass counts by 7B1a2)
# demonstration of new extended geochemical data option in fetchKSSL
# andrew brown; 2020/3/19

library(dplyr)
library(aqp)
library(soilDB)

# select some MLRAs to compare
mlras <- list('127','147','148')

# further select some elements (patterns to search in column names)
element.list <- list("lead", "copper", "zinc", "mercury")

# use lapply to iterate over MLRAs and do fetchKSSL
# with returnGeochemicalData = TRUE the result is a list, not just a SPC
# For each MLRA:
#  - SPC is in the $SPC element
#  - geochem is in $geochem
#  - optical in $optical
#  - xray in $xrd_thermal

f <- lapply(mlras, function(a.mlra) soilDB::fetchKSSL(mlra = a.mlra, returnGeochemicalData = TRUE))

# extract just SPC component from each MLRA and combine them with union
spc <- aqp::union(lapply(f, function(x) x$SPC))

# extract geochem table from each MLRA and combine with rbind
geochem <- do.call('rbind', lapply(f, function(x) x$geochem))

# keep just rows with trace elemental analysis counts
geochem <- dplyr::filter(geochem, trace_element_method == "4H1a")

# remove non-target elements, keep identifier columns
keep.idx <- do.call('c', lapply(element.list, 
                               function(element) grep(element, colnames(geochem))))
geochem<- geochem[, c(1:3,keep.idx)]

# transfer mlra labels to geochem table
spc$hzmlra <- denormalize(spc, 'mlra')
h <- horizons(spc)
geochem$mlra <- h[match(geochem$labsampnum, h$labsampnum), 'hzmlra']

# do some simple boxplots, one for each element
lapply(as.list(colnames(geochem[,-c(1:3, length(colnames(geochem)))])), function(element) {
  print(sprintf("%s ~ mlra", element))
  boxplot(data=geochem, as.formula(sprintf("%s ~ mlra", element)), xlab="MLRA")
})

Produces plots like this:

lead_by_mlra

brownag commented 4 years ago

I worked up a more complicated example involving glass count comparisons (also by MLRA).

I'm thinking these will make a nice sequence of tutorials. The geochem is perhaps the most straightforward table to work with because the analyses all tend to use common methods and prep code, with no difference in size fraction.

Glass/grain counts are more complicated -- in that you need to deal with multiple size fractions, which requires a linkage back to the regular KSSL PSDA data. A bit more involved.

I still haven't come up with any demos for XRD/thermal data, but I was thinking of trying to do something with Vertisols, or Vertic properties.

brownag commented 4 years ago

glass_by_mlra

brownag commented 4 years ago

Some demo scripts live here: https://github.com/brownag/sANDREWbox/tree/master/KSSLExtended

Kramdog commented 4 years ago

Good stuff Andrew . Thanks