A Unified Interface to SCAN/SNOTEL Data
D.E. Beaudette and J. Skovlin
2024-09-18
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(terra)
library(spData)
# get basic station metadata for several SCAN/SNOTEL sites
m <- SCAN_site_metadata(site.code = c(2072, 356, 2148, 2187))
# select columns from station metadata
knitr::kable(m[, c('Site', 'Name', 'Network', 'HUC', 'upedonid')], row.names = FALSE)| Site | Name | Network | HUC | upedonid |
|---|---|---|---|---|
| 2072 | Eros Data Center | SCAN | 101702031402 | S2003SD099001 |
| 356 | Blue Lakes | SNTL | 180400120101 | NA |
| 2148 | Jordan Valley Cwma | SCAN | 170501080405 | NA |
| 2187 | Deep Springs | SCAN | 180902010102 | S2012CA027002 |
# result is a list
x <- fetchSCAN(site.code = c(574), year = c(2016))
# print a list of sensor types and associated metadata
names(x)## [1] "SMS" "STO" "SAL" "TAVG" "TMIN" "TMAX" "PRCP" "PREC"
## [9] "SNWD" "WTEQ" "WDIRV" "WSPDV" "LRADT" "metadata"Metadata for all Stations
data('SCAN_SNOTEL_metadata', package = 'soilDB')
# number of stations by network
table(SCAN_SNOTEL_metadata$Network)##
## CSCAN SCAN SNTL SNTLT
## 21 210 910 45data("us_states")
us_states <- vect(us_states)
us_states <- project(us_states, 'epsg:5070')
s <- vect(SCAN_SNOTEL_metadata, geom = c('Longitude', 'Latitude'), crs = 'epsg:4269')
s <- project(s, 'epsg:5070')plot(us_states, axes = FALSE, main = 'SCAN, SNOTEL, CSCAN, SNOWLITE Stations')
points(s[s$Network %in% c('SCAN', 'CSCAN')], col = 'firebrick')
points(s[s$Network %in% c('SNTL', 'SNTLT')], col = 'royalblue')
legend(x = -2029093, y = 801043.2, legend = c('SCAN / CSCAN', 'SNOTEL / SNOWLITE'), pch = 21, pt.bg = c('firebrick', 'royalblue'), bty = 'n')
plot(us_states, axes = FALSE, main = 'SCAN, SNOTEL, CSCAN, SNOWLITE Stations')
points(s, col = 'royalblue')
points(s[!is.na(s$pedlabsampnum), ], col = 'firebrick')
legend(x = -2029093, y = 801043.2, legend = c('All Stations', 'Linked to Lab Data'), pch = 21, pt.bg = c('royalblue', 'firebrick'), 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 9 variables:
## $ STO :'data.frame': 7296 obs. of 9 variables:
## $ SAL :'data.frame': 0 obs. of 0 variables
## $ TAVG :'data.frame': 1459 obs. of 9 variables:
## $ TMIN :'data.frame': 1459 obs. of 9 variables:
## $ TMAX :'data.frame': 1459 obs. of 9 variables:
## $ PRCP :'data.frame': 1319 obs. of 9 variables:
## $ PREC :'data.frame': 1323 obs. of 9 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 19 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 9 variables:
## $ STO :'data.frame': 5843 obs. of 9 variables:
## $ SAL :'data.frame': 0 obs. of 0 variables
## $ TAVG :'data.frame': 1460 obs. of 9 variables:
## $ TMIN :'data.frame': 1460 obs. of 9 variables:
## $ TMAX :'data.frame': 1460 obs. of 9 variables:
## $ PRCP :'data.frame': 581 obs. of 9 variables:
## $ PREC :'data.frame': 1316 obs. of 9 variables:
## $ SNWD :'data.frame': 731 obs. of 9 variables:
## $ WTEQ :'data.frame': 713 obs. of 9 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 19 variables:# deeper look
str(x$SMS)## 'data.frame': 5752 obs. of 9 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 ...
## $ datetime : POSIXct, format: "2015-01-01 12:00:00" "2015-01-02 12:00:00" "2015-01-03 12:00:00" ...# 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 9 variables:
## $ STO :'data.frame': 27933 obs. of 9 variables:
## $ SAL :'data.frame': 0 obs. of 0 variables
## $ TAVG :'data.frame': 2918 obs. of 9 variables:
## $ TMIN :'data.frame': 2917 obs. of 9 variables:
## $ TMAX :'data.frame': 2919 obs. of 9 variables:
## $ PRCP :'data.frame': 2918 obs. of 9 variables:
## $ PREC :'data.frame': 2927 obs. of 9 variables:
## $ SNWD :'data.frame': 1407 obs. of 9 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 19 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 datetime year doy
## 1 2144 2008-01-01 12:00 2008 93 2.1 5 SMS.I_2 2008-01-01 12:00:00 2008 1
## 2 2144 2008-01-02 12:00 2008 94 1.8 5 SMS.I_2 2008-01-02 12:00:00 2008 2
## 3 2144 2008-01-03 12:00 2008 95 2.6 5 SMS.I_2 2008-01-03 12:00:00 2008 3
## 4 2144 2008-01-04 12:00 2008 96 3.3 5 SMS.I_2 2008-01-04 12:00:00 2008 4
## 5 2144 2008-01-05 12:00 2008 97 5.1 5 SMS.I_2 2008-01-05 12:00:00 2008 5
## 6 2144 2008-01-06 12:00 2008 98 12.1 5 SMS.I_2 2008-01-06 12:00:00 2008 6head(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 <- fetchLDM(x = 'S08NV003003', what = 'upedonid')
# 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.4 and
soilDB version 2.8.5.