Introduction

This document demonstrates how to use the soilDB package to download climate data from the SCAN/SNOTEL network. Stay tuned for updates and more detailed examples.

# run these commands in the R console
install.packages('latticeExtra', dep=TRUE)
install.packages('plyr', dep=TRUE)
install.packages('rvest', dep=TRUE)
install.packages('httr', dep=TRUE)
install.packages('reshape2', dep=TRUE)
install.packages('soilDB', dep=TRUE)

install.packages('devtools', dep=TRUE)
devtools::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade_dependencies=FALSE)
devtools::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)

Quick Example: getting sensor metadata

Each SCAN / SNOTEL station may have a different set of sensors.

# load required packages
library(aqp)
library(soilDB)
library(reshape2)
library(plyr)
library(latticeExtra)

# get basic sensor metadata for several SCAN/SNOTEL sites
m <- SCAN_sensor_metadata(site.code=c(2072,356,2148,2187))

# print a list of sensor types by station
dlply(m, 'site.code', function(i) unique(i$Ecode))

## $`356`
##  [1] "WTEQ" "PREC" "TOBS" "TMAX" "TMIN" "TAVG" "SNWD" "SMS"  "STO"  "SAL"  "RDC"  "BATT"
## 
## $`2072`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC" "LRADT"
## 
## $`2148`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC"
## 
## $`2187`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC" "LRADT"
## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   site.code
## 1       356
## 2      2072
## 3      2148
## 4      2187

Getting Data from a Single Station

Lets get some recent (daily) climate data from the Rogers Farm #1 SCAN station; site number 2001. The fetchSCAN() function does all of the work. Notice that the resulting object is a list, each element is a suite of sensor data. For example, the SMS element contains soil moisture from several depths

# fetchSCAN can accept multiple years
x <- fetchSCAN(site.code=2001, year=c(2013,2014,2015,2016))

# check the results
str(x, 1)

## List of 11
##  $ SMS  :'data.frame':   7325 obs. of  5 variables:
##  $ STO  :'data.frame':   7325 obs. of  5 variables:
##  $ SAL  :'data.frame':   0 obs. of  0 variables
##  $ TAVG :'data.frame':   1465 obs. of  5 variables:
##  $ PRCP :'data.frame':   1465 obs. of  5 variables:
##  $ PREC :'data.frame':   1465 obs. of  5 variables:
##  $ SNWD :'data.frame':   0 obs. of  0 variables
##  $ WTEQ :'data.frame':   0 obs. of  0 variables
##  $ WDIRV:'data.frame':   0 obs. of  0 variables
##  $ WSPDV:'data.frame':   0 obs. of  0 variables
##  $ LRADT:'data.frame':   0 obs. of  0 variables

# tabulate number of soil moisture measurements per depth (cm)
table(x$SMS$depth)

# get unique set of soil moisture sensor depths
sensor.depths <- unique(x$SMS$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(x$SMS$Date)), as.Date(max(x$SMS$Date)), by='3 months')

# plot soil moisture time series, panels are depth
xyplot(value ~ Date | factor(depth), data=x$SMS, as.table=TRUE, type=c('l','g'), strip=strip.custom(bg=grey(0.80)), layout=c(1,length(sensor.depths)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), ylab='Volumetric Water Content', main='Soil Moisture at Several Depths')

Getting Data from Multiple Stations

That was interesting, but most of the time we want to make comparisons between stations. Note that the format of data returned makes it possible "stack" data from several stations that do not share the same collection of sensors.

# get more data
x <- fetchSCAN(site.code=c(356, 2072), year=c(2015, 2016))

# same format as before, sensor data are "stacked" within each list element
str(x, 1)

## List of 11
##  $ SMS  :'data.frame':   5864 obs. of  5 variables:
##  $ STO  :'data.frame':   5864 obs. of  5 variables:
##  $ SAL  :'data.frame':   0 obs. of  0 variables
##  $ TAVG :'data.frame':   1466 obs. of  5 variables:
##  $ PRCP :'data.frame':   733 obs. of  5 variables:
##  $ PREC :'data.frame':   1466 obs. of  5 variables:
##  $ SNWD :'data.frame':   733 obs. of  5 variables:
##  $ WTEQ :'data.frame':   733 obs. of  5 variables:
##  $ WDIRV:'data.frame':   0 obs. of  0 variables
##  $ WSPDV:'data.frame':   0 obs. of  0 variables
##  $ LRADT:'data.frame':   0 obs. of  0 variables

# deeper look
str(x$SMS)

## 'data.frame':    5864 obs. of  5 variables:
##  $ Site     : int  356 356 356 356 356 356 356 356 356 356 ...
##  $ Date     : Date, format: "2015-01-01" "2015-01-02" "2015-01-03" ...
##  $ value    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ depth    : num  5 5 5 5 5 5 5 5 5 5 ...
##  $ sensor.id: Factor w/ 5 levels "SMS.I_2","SMS.I_8",..: 1 1 1 1 1 1 1 1 1 1 ...

# get unique set of soil moisture sensor depths
sensor.depths <- unique(x$SMS$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(x$SMS$Date)), as.Date(max(x$SMS$Date)), by='3 months')

# plot soil moisture time series, panels are depth
xyplot(value ~ Date | factor(depth), groups=factor(Site), data=x$SMS, as.table=TRUE, type=c('l','g'), auto.key=list(columns=2, lines=TRUE, points=FALSE), strip=strip.custom(bg=grey(0.80)), layout=c(1,length(sensor.depths)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), ylab='Volumetric Water Content', main='Soil Moisture at Several Depths')

# combine sensors
g <- make.groups('soil moisture'=x$SMS, 'soil temperature'=x$STO)

# get unique set of soil moisture sensor depths
sensor.depths <- unique(g$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(g$Date)), as.Date(max(g$Date)), by='3 months')

# plot soil moisture time series, panels are depth
p <- xyplot(value ~ Date | factor(Site) + which, groups=factor(depth), data=g, as.table=TRUE, type=c('l','g'), auto.key=list(columns=length(sensor.depths), lines=TRUE, points=FALSE), strip=strip.custom(bg=grey(0.80)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y"), y=list(relation='free', rot=0)), ylab='', main='Soil Moisture and Temperature at Several Depths')

useOuterStrips(p)

# plot a single suite of sensors
# compare sites / depths
# no SAL sensors for these sites
# xyplot(value ~ Date | factor(depth), groups=factor(Site), data=x$SAL, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), layout=c(1,length(sensor.depths)), scales=list(alternating=3, y=list(relation='free', rot=0)))

xyplot(value ~ Date, groups=factor(Site), data=x$TAVG, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Average Air Temperature by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$PRCP, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Daily Precipitation by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$PREC, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Accumulated Precipitation by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$WTEQ, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Daily Snow Water Equivalent by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$SNWD, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Daily Snow Depth (in) by Site')

Case Study: Pine Nut SCAN Site

Get SCAN data from the Pine Nut site(2144) for years 2008 through 2015. Note that this site has multiple soil temperature sensors installed at each depth.

x <- fetchSCAN(site.code=c(2144), year=2008:2015)

# check available data
str(x, 1)

## List of 11
##  $ SMS  :'data.frame':   14650 obs. of  5 variables:
##  $ STO  :'data.frame':   29300 obs. of  5 variables:
##  $ SAL  :'data.frame':   0 obs. of  0 variables
##  $ TAVG :'data.frame':   2930 obs. of  5 variables:
##  $ PRCP :'data.frame':   2930 obs. of  5 variables:
##  $ PREC :'data.frame':   2930 obs. of  5 variables:
##  $ SNWD :'data.frame':   1465 obs. of  5 variables:
##  $ WTEQ :'data.frame':   0 obs. of  0 variables
##  $ WDIRV:'data.frame':   0 obs. of  0 variables
##  $ WSPDV:'data.frame':   0 obs. of  0 variables
##  $ LRADT:'data.frame':   0 obs. of  0 variables

## still working on a solution for this issue: https://github.com/ncss-tech/soilDB/issues/14
# ! multiple soil temperature sensors installed at each depth
table(sensor.id=x$STO$sensor.id, sensor.depth=x$STO$depth)

# extract SMS data to a new data.frame
sms <- x$SMS

# add parts of date for plotting
sms$year <- format(sms$Date, "%Y")
sms$doy <- as.numeric(format(sms$Date, "%j"))

# cast long -> wide format for later use
sms.wide <- dcast(sms, Site + Date ~ depth, value.var = 'value', mean, na.rm=TRUE)

# check:
head(sms)

head(sms.wide)

# make some pretty colors for plotting daily VWC
cols <- colorRampPalette(brewer.pal(11, 'RdYlBu'), space='Lab', interpolate='spline', bias=2)

# plot VWC by year/day, panels by depth
levelplot(value ~ doy * factor(year) | factor(depth), main='Daily Mean VWC (%)',
          data=sms, layout=c(1,5), col.regions=cols,
          colorkey=list(space='top'), as.table=TRUE, scales=list(alternating=3, cex=0.75, x=list(tick.number=20)), 
          par.strip.text=list(cex=0.85), strip=strip.custom(bg='grey'), 
          xlab='Julian Day', ylab='Year')

# plot VWC by year/day, panels by year with depths stacked
levelplot(value ~ doy * rev(factor(depth)) | factor(year), main='Daily Mean VWC (%)',
          data=sms, layout=c(1,8), col.regions=cols, 
          colorkey=list(space='top'), as.table=TRUE, scales=list(alternating=3, cex=0.75, x=list(tick.number=20)), 
          par.strip.text=list(cex=0.85), strip=strip.custom(bg='grey'), 
          xlab='Julian Day', ylab='Sensor Depth (cm)', ylim=rev(levels(factor(sms$depth))))

# get KSSL data
s <- fetchKSSL(pedon_id = 'S08NV003003')

# VG parameters now returned by fetchKSSL() as of 2016-11-17
(horizons(s)[, c('hzn_desgn', 'hzn_top', 'hzn_bot', 'theta_r', 'theta_s', 'alpha', 'npar')])

# get VG parameters for sensor at 20cm depth
s.vg.hz <- slice(s, 20 ~ theta_r + theta_s + alpha + npar, just.the.data = TRUE)

# get VG curve and inverse function
# this is based on several assumptions about the source data... details pending
vg.hz <- KSSL_VG_model(VG_params = s.vg.hz)

# check: looks good
p.model <- xyplot(phi ~ theta, data=vg.hz$VG_curve, type=c('l', 'g'), scales=list(alternating=3, x=list(tick.number=10), y=list(log=10, tick.number=10)), yscale.components=yscale.components.logpower, ylab=expression("Suction " (kPa)), xlab=expression("Volumetric Water Content " (cm^3/cm^3)), par.settings = list(plot.line=list(col='RoyalBlue', lwd=2)))

update(p.model, main='Estimated Water Retention Curve at 20cm', sub='van Genuchten Model Parameters fit by USDA-ARS Rosetta')

# add quantiles from measured values
sms.q <- data.frame(theta=quantile(sms.wide$`20`/100, probs=c(0, 0.05, 0.25, 0.5, 0.75, 0.95, 1), na.rm=TRUE))
sms.q$phi <- vg.hz$VG_inverse_function(sms.q$theta)

p.points <- xyplot(phi ~ theta, data=sms.q, type='p', par.settings = list(plot.symbol=list(col='orange', pch=16, cex=1.25)), scales=list(alternating=3, x=list(tick.number=10), y=list(log=10, tick.number=10)), yscale.components=yscale.components.logpower, ylab=expression("Suction " (kPa)), xlab=expression("Volumetric Water Content " (cm^3/cm^3)), panel=function(x, y, ...) {
  panel.xyplot(x, y, ...)
  panel.text(x, y, row.names(sms.q), adj=c(0,-0.75), cex=0.65)
})

update(p.model + p.points, main='Select Percentiles of Data Measured at 20cm', sub='van Genuchten Model Parameters fit by USDA-ARS Rosetta')

# split soil moisture data by sensor depth
sms.list <- split(sms, sms$depth)

# convert VWC -> phi (KPa)
sms.phi <- ldply(sms.list, function(i) {
  # i is a DF associated with a single sensor depth
  this.depth <- unique(i$depth)
  # get VG params for this depth
  fm <- as.formula(paste0(this.depth, ' ~ theta_r + theta_s + alpha + npar'))
  s.vg.hz <- slice(s, fm, just.the.data = TRUE)

  # get VG curve and inverse function
  # this is based on several assumptions about the source data...
  vg.hz <- KSSL_VG_model(VG_params = s.vg.hz)
  
  # estimate phi
  if(!is.null(vg.hz$VG_inverse_function))
    i$phi <- vg.hz$VG_inverse_function(i$value / 100.0)
  else
    i$phi <- NA
  
  # done
  return(i)
})


# classify each day as <= 1500 kPa (15 bar)
sms.phi$phi_leq_15_bar <- factor(sms.phi$phi <= 1500, levels=c('FALSE', 'TRUE'))

# new color palette
cols.phi <- colorRampPalette(rev(brewer.pal(11, 'RdYlBu')), space='Lab', interpolate='spline', bias=1.5)

# plot daily mean water potential in log10-kPa
levelplot(log(phi, base=10) ~ doy * factor(year) | factor(depth), main='Daily Mean Water Potential log10(-kPa)',
          data=sms.phi, layout=c(1,5), col.regions=cols.phi,
          colorkey=list(space='top'), as.table=TRUE, scales=list(alternating=3, cex=0.75, x=list(tick.number=20)), 
          par.strip.text=list(cex=0.85), strip=strip.custom(bg='grey'), 
          xlab='Julian Day', ylab='Year')

# plot daily mean water potential in log10-kPa by year and sensor stack
levelplot(log(phi, base=10) ~ doy * rev(factor(depth)) | factor(year), main='Daily Mean Water Potential log10(-kPa)',
          data=sms.phi, layout=c(1,8), col.regions=cols.phi,
          colorkey=list(space='top'), as.table=TRUE, scales=list(alternating=3, cex=0.75, x=list(tick.number=20)), 
          par.strip.text=list(cex=0.85), strip=strip.custom(bg='grey'), 
          xlab='Julian Day', ylab='Sensor Depth (cm)', ylim=rev(levels(factor(sms.phi$depth))))

levelplot(phi_leq_15_bar ~ doy * factor(year) | factor(depth), main='Daily Mean Water Potential <= 15 bar (1500 kPa)',
          data=sms.phi, layout=c(1,5), col.regions=c(grey(0.85), 'RoyalBlue'), cuts=1,
          colorkey=FALSE, as.table=TRUE, scales=list(alternating=3, cex=0.75, x=list(tick.number=20)), 
          par.strip.text=list(cex=0.85), strip=strip.custom(bg='grey'), 
          xlab='Julian Day', ylab='Year')