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.

Note

This interface is very much a work in progress. There are some SCAN / SNOTEL site with multiple (above-ground) sensors per sensor prefix, which can lead to confusing results. For now the first above-ground sensor is selected for each above-ground sensor prefix: c('TAVG', 'PRCP', 'PREC', 'SNWD', 'WTEQ', 'WDIRV', 'WSPDV', 'LRADT').

There are some sites with multiple below-ground sensors installed at the same depth. A message is printed when this happens. You can use the sensor.id column to access specific sensors. We are working on a solution.

Installation

Install packages used in this tutorial from CRAN.

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

Sensor Metadata

Metadata for Select Stations

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

# load required packages
library(aqp)
library(soilDB)
library(sharpshootR)
library(reshape2)
library(latticeExtra)
library(tactile)
library(maps)

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

knitr::kable(m, row.names = FALSE)
Name Site State Network County Elevation_ft Latitude Longitude HUC climstanm upedonid pedlabsampnum
Blue Lakes 356 California SNOTEL Alpine 8067 38.60800 -119.92437 180400120101 NA NA NA
Deep Springs 2187 California SCAN Inyo 5399 37.37222 -117.97383 180902010102 Deep Springs S2012CA027002 13N97485
Eros Data Center 2072 South Dakota SCAN Minnehaha 1602 43.73835 -96.61455 101702031402 Eros Data Center S2003SD099001 03N0688
Jordan Valley Cwma 2148 Idaho SCAN Owyhee 4508 42.94863 -117.01060 170501080405 NA NA NA
# result is a list
x <- fetchSCAN(site.code=c(574), year=c(2016))

# print a list of sensor types
names(x[-1])
##  [1] "STO"      "SAL"      "TAVG"     "TMIN"     "TMAX"     "PRCP"     "PREC"     "SNWD"    
##  [9] "WTEQ"     "WDIRV"    "WSPDV"    "LRADT"    "metadata"

Metadata for all Stations

data('SCAN_SNOTEL_metadata', package = 'soilDB')

map('state', lwd=2)
title(main='SCAN / SNOTEL Sites')
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=Network == 'SCAN', bg='DarkRed', pch=21)
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=Network == 'SNOTEL', bg='RoyalBlue', pch=21)

legend('bottomleft', legend=c('SCAN', 'SNOTEL'), pch=21, pt.bg=c('DarkRed', 'RoyalBlue'), bty='n')

map('state', lwd=2)
title(main='SCAN / SNOTEL Sites')
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, bg='RoyalBlue', pch=21, cex=0.85)
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=!is.na(pedlabsampnum), bg='DarkRed', pch=22)

legend('bottomleft', legend=c('All Stations', 'Linked to Lab Data'), pch=c(21, 22), pt.bg=c('RoyalBlue', 'DarkRed'), bty='n')

Sensor and Site Data for 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(2014, 2015, 2016, 2017))

# check the results
str(x, 1)
## List of 14
##  $ SMS     :'data.frame':    6656 obs. of  8 variables:
##  $ STO     :'data.frame':    7296 obs. of  8 variables:
##  $ SAL     :'data.frame':    0 obs. of  0 variables
##  $ TAVG    :'data.frame':    1459 obs. of  8 variables:
##  $ TMIN    :'data.frame':    1459 obs. of  8 variables:
##  $ TMAX    :'data.frame':    1459 obs. of  8 variables:
##  $ PRCP    :'data.frame':    1319 obs. of  8 variables:
##  $ PREC    :'data.frame':    1323 obs. of  8 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
##  $ metadata:'data.frame':    1 obs. of  12 variables:
# tabulate number of soil moisture measurements per depth (cm)
table(x$SMS$depth)
## 
##    5   10   20   51  102 
## 1458 1458 1434 1429  877
# 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, 
  par.settings = tactile.theme(),
  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'
)

Visualizing sensor depths within the context of a soil profile

One way to get a visual diagram of where the SCAN soil moisture and temperature sensors occur within the context of the soil profile horizons is to use the linkage in the attached metadata to join the sensor depths to the profile horizonation.

# get more data
x <- fetchSCAN(site.code = 574, year = 2016)
s <- fetchKSSL(pedlabsampnum = na.omit(x$metadata$pedlabsampnum), returnMorphologicData = FALSE)

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

# plot profile of lab data horizons then over plot with location of sensor depths within them
par(mar = c(0, 0, 3, 1))
plotSPC(s, label = 'pedon_id', id.style = 'top', cex.names = 0.75, width = 0.15, depth.axis = list(line = -3.5))
title(paste('Sensor Depths for SCAN site', unique(x$SMS$Site), sep=" "), line = 0.5, cex.main = 1)

# over plot the sensor depths on the previous plot
lsp <- get("last_spc_plot", envir = aqp.env)

points(
  x = rep(lsp$x0, times = length(sensor.depths)), 
  y = sensor.depths, pch = 15, cex = 1.5, col = 'royalblue'
)

text(
  x = rep(lsp$x0, times = length(sensor.depths)), 
  y = sensor.depths, 
  labels = sensor.depths, 
  cex = 0.75, 
  pos = 3
)

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 14
##  $ SMS     :'data.frame':    5752 obs. of  8 variables:
##  $ STO     :'data.frame':    5843 obs. of  8 variables:
##  $ SAL     :'data.frame':    0 obs. of  0 variables
##  $ TAVG    :'data.frame':    1460 obs. of  8 variables:
##  $ TMIN    :'data.frame':    1460 obs. of  8 variables:
##  $ TMAX    :'data.frame':    1460 obs. of  8 variables:
##  $ PRCP    :'data.frame':    581 obs. of  8 variables:
##  $ PREC    :'data.frame':    1316 obs. of  8 variables:
##  $ SNWD    :'data.frame':    731 obs. of  8 variables:
##  $ WTEQ    :'data.frame':    713 obs. of  8 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
##  $ metadata:'data.frame':    2 obs. of  12 variables:
# deeper look
str(x$SMS)
## 'data.frame':    5752 obs. of  8 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" ...
##  $ Time      : chr  "12:00" "12:00" "12:00" "12:00" ...
##  $ water_year: num  2015 2015 2015 2015 2015 ...
##  $ water_day : int  93 94 95 96 97 98 99 100 101 102 ...
##  $ 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, 
  par.settings = tactile.theme(),
  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'
)

Example Figures

# 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, 
  par.settings = tactile.theme(),
  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 (deg. C)'
)

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, 
  par.settings = tactile.theme(),
  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 (deg. C)'
)

xyplot(
  value ~ Date, 
  groups = factor(Site), 
  data = x$PRCP, 
  as.table = TRUE, 
  par.settings = tactile.theme(),
  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 (inches)'
)

xyplot(
  value ~ Date, 
  groups = factor(Site), 
  data = x$PREC, 
  as.table = TRUE, 
  par.settings = tactile.theme(),
  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 (inches)'
)

xyplot(
  value ~ Date, 
  groups = factor(Site), 
  data = x$WTEQ, 
  as.table = TRUE, 
  par.settings = tactile.theme(),
  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 (inches)'
)

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

Estimating MAST, MSST, MWST, STR

This is a mostly-complete example of how to estimate mean anual soil temperature (MAST), mean summer soil temperature (MSST), mean winter soil temperature (MWST), and soil temperature regime (STR) using some of the internals from the fetchHenry function.

# get period of record for several sites
# note: multiple below-ground sensors
x <- fetchSCAN(site.code = c(2096, 2097, 2099, 2101), year = 2000:2017)

## 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)

# subset data at 50cm (51cm is close enough, rounding error)
soiltemp <- subset(x$STO, subset = depth == 51)

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


# nice colors
cols <- colorRampPalette(
  hcl.colors(n = 100, palette = 'spectral', rev = TRUE), 
  space = 'Lab', 
  interpolate = 'spline'
)


# plot soiltemp by year/day, panels by sensor ID
useOuterStrips(
  levelplot(
    value ~ doy * factor(year) | sensor.id + factor(Site), 
    main = 'Daily Mean Soil Temperature at 50cm (deg C)',
    data = soiltemp, col.regions = cols, xlim = c(1,365),
    colorkey = list(space='top'), 
    as.table = TRUE, 
    par.settings = tactile.theme(),
    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'
  )
)

# ... something isn't correct with sensor STO.I-2_20
# plot just STO.I_20
# ... better
levelplot(
  value ~ doy * factor(year) | factor(Site), 
  main = 'Daily Mean Soil Temperature at 50cm (deg C)',
  subset = sensor.id == 'STO.I_20',
  data = soiltemp, 
  col.regions = cols, 
  xlim = c(1, 365),
  colorkey = list(space='top'), 
  as.table = TRUE, 
  par.settings = tactile.theme(),
  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'
)

# subset to single sensor
soiltemp <- subset(soiltemp, subset = sensor.id == 'STO.I_20')


## TODO: generalize to accept SCAN / SNOTEL data as-is

# for now, make compatible
soiltemp$sid <- soiltemp$Site
soiltemp$sensor_value <- soiltemp$value
soiltemp$date_time <- soiltemp$Date

# prep dates and pad with missing DOY
soiltemp <- soilDB:::.formatDates(soiltemp, gran = 'day', pad.missing.days = TRUE)

# compute soil temperature summaries, assuming 50cm depth
d <- summarizeSoilTemperature(soiltemp)

knitr::kable(d)

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 = 2144, year = 2008:2015)

# check available data
str(x, 1)
## List of 14
##  $ SMS     :'data.frame':    14259 obs. of  8 variables:
##  $ STO     :'data.frame':    27933 obs. of  8 variables:
##  $ SAL     :'data.frame':    0 obs. of  0 variables
##  $ TAVG    :'data.frame':    2918 obs. of  8 variables:
##  $ TMIN    :'data.frame':    2917 obs. of  8 variables:
##  $ TMAX    :'data.frame':    2919 obs. of  8 variables:
##  $ PRCP    :'data.frame':    2918 obs. of  8 variables:
##  $ PREC    :'data.frame':    2927 obs. of  8 variables:
##  $ SNWD    :'data.frame':    1407 obs. of  8 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
##  $ metadata:'data.frame':    1 obs. of  12 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)
##             sensor.depth
## sensor.id       5   10   20   51  102
##   STO.I_2    2918    0    0    0    0
##   STO.I_4       0 1944    0    0    0
##   STO.I_8       0    0 2918    0    0
##   STO.I_20      0    0    0 2918    0
##   STO.I_40      0    0    0    0 2918
##   STO.I-2_2  2920    0    0    0    0
##   STO.I-2_4     0 2639    0    0    0
##   STO.I-2_8     0    0 2920    0    0
##   STO.I-2_20    0    0    0 2920    0
##   STO.I-2_40    0    0    0    0 2918
# 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)
##   Site       Date  Time water_year water_day value depth sensor.id year doy
## 1 2144 2008-01-01 12:00       2008        93   2.1     5   SMS.I_2 2008   1
## 2 2144 2008-01-02 12:00       2008        94   1.8     5   SMS.I_2 2008   2
## 3 2144 2008-01-03 12:00       2008        95   2.6     5   SMS.I_2 2008   3
## 4 2144 2008-01-04 12:00       2008        96   3.3     5   SMS.I_2 2008   4
## 5 2144 2008-01-05 12:00       2008        97   5.1     5   SMS.I_2 2008   5
## 6 2144 2008-01-06 12:00       2008        98  12.1     5   SMS.I_2 2008   6
head(sms.wide)
##   Site       Date    5   10   20  51 102
## 1 2144 2008-01-01  2.1  4.4  8.8 9.2 2.0
## 2 2144 2008-01-02  1.8  4.1  8.1 9.3 1.6
## 3 2144 2008-01-03  2.6  5.2  8.2 9.2 1.6
## 4 2144 2008-01-04  3.3  6.2  8.2 9.0 NaN
## 5 2144 2008-01-05  5.1  7.7  8.2 9.1 2.0
## 6 2144 2008-01-06 12.1 18.8 12.7 9.1 2.1
# make some pretty colors for plotting daily VWC
cols <- colorRampPalette(
  hcl.colors(n = 100, palette = 'spectral', rev = TRUE), 
  space = 'Lab', 
  interpolate = 'spline'
)

# 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, 
  par.settings = tactile.theme(),
  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, 
  par.settings = tactile.theme(),
  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)))
)

Water Retention Curve Development

Get KSSL data for this site, via user pedon ID. This will include estimated parameters for the van Genuchten model which can be used to derive a water retention curve.

Previously the VG parameters had to be obtained manually from these reports.

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

# VG parameters now returned by fetchKSSL() as of 2016-11-17
knitr::kable(horizons(s)[, c('hzn_desgn', 'hzn_top', 'hzn_bot', 'theta_r', 'theta_s', 'alpha', 'npar')])
hzn_desgn hzn_top hzn_bot theta_r theta_s alpha npar
A 0 9 0.0337216 0.4864061 -1.581517 0.1227247
ABk 9 31 0.0445697 0.4701070 -1.605478 0.1267858
Btk1 31 48 0.0516780 0.4647200 -1.795169 0.1317818
Btk2 48 83 0.0562656 0.4179194 -2.058838 0.1377132
2Bk 83 152 NA NA NA NA

Conversion of VWC to Matric Potential: Single Depth

Generate water retention curve for sensor at 20cm. This requires extracting VG parameters from the KSSL data from the associated horizon. See ?KSSL_VG_model or this tutorial for details.

# get VG parameters for sensor at 20cm depth
s.vg.hz <- dice(s, 20 ~ theta_r + theta_s + alpha + npar, SPC = FALSE)

# 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 = tactile.theme(plot.line = list(lwd = 2))
)

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

Lets add some points on the water retention curve from the measured VWC values.

# 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'
)

Add permanent wilting point, field capacity, and saturation.

# PWP, FC, SAT (close to 0 suction)
sms.q <- data.frame(phi = c(1500, 33.33, 1), na.rm = TRUE)
sms.q$theta <- vg.hz$VG_function(sms.q$phi)
row.names(sms.q) <- c('PWP', 'FC', 'SAT')

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 = 'Critical Points, 20cm', 
  sub = 'van Genuchten Model Parameters fit by USDA-ARS Rosetta'
)

Estimate soil moisture state (Newhall model conventions) using water retention curve benchmarks above.

# VWC are encoded as percent (0-100)
# must convert to fraction (0-1)
# using water retention curve benchmarks from above
sms$state <- estimateSoilMoistureState(
  sms$value / 100, 
  U = 0, 
  sat = 0.45, 
  fc = 0.25, 
  pwp = 0.10, 
  style = 'newhall'
)

table(sms$state)
## 
##       dry     moist saturated 
##      9138      4958       163
# colors / style
ll <- levels(sms$state)
n.states <- length(ll)
ms.colors <- hcl.colors(n = n.states, palette = 'spectral', rev = TRUE)

# color palette function
ms.colors.f <- colorRampPalette(
  ms.colors,
  space = 'Lab',
  interpolate = 'spline'
)

suppressWarnings(
  trellis.par.set(list(superpose.polygon = list(col = ms.colors, border = ms.colors)))
)

sK <- simpleKey(text = ll, space = 'top', columns = n.states, rectangles = TRUE, points = FALSE, cex = 1)

levelplot(
  state ~ doy * factor(year) | factor(depth), main='Daily Soil Moisture State',
  data=sms, layout=c(1,5), 
  col.regions = ms.colors.f,
  key = sK,
  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'
)

Conversion of VWC to Matric Potential: All Depths

This time we will iterate over all sensors, develop water retention curves, and convert VWV to water potential values. Note that this approach will only work with a SoilProfileCollection object that contains a single profile.

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

# convert VWC -> phi (KPa)
sms.phi <- lapply(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 <- dice(s, fm, SPC = FALSE)

  # 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))
    # convert percent -> fraction
    i$phi <- vg.hz$VG_inverse_function(i$value / 100.0)
  else
    i$phi <- NA
  
  # done
  return(i)
})

# flatten
sms.phi <- do.call('rbind', sms.phi)


# 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(
  hcl.colors(n = 11, palette = 'spectral', rev = TRUE), 
  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')


This document is based on aqp version 2.0.2 and soilDB version 2.7.8.