Soil Series Climate Summary Evaluation
D.E. Beaudette
2019-02-08
Introduction
Setup
With a recent version of R (>= 2.15), it is possible to get all of the packages that this tutorial depends on via:
# stable packages from CRAN
install.packages('aqp', dep=TRUE)
install.packages('soilDB', dep=TRUE)
install.packages('sharpshootR', dep=TRUE)
install.packages('latticeExtra', dep=TRUE)
install.packages('ggplot2', dep=TRUE)
install.packages('raster', dep=TRUE)
# latest versions from GitHub
devtools::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)
devtools::install_github("ncss-tech/sharpshootR", dependencies=FALSE, upgrade_dependencies=FALSE)Using this Tutorial
This tutorial can be used as a stand-alone demonstration of the functionality provided by fetchOSD and OSDquery, or as a template for your own project. Copy any paste blocks of code sequentially into your own script if you would like to follow-along.
First, you will need to load some packages. Be sure to install these if missing, and install the latest versions of soilDB and sharpshootR as much of this tutorial depends on recent updates.
library(soilDB)
library(raster)
library(maps)
library(sharpshootR)
library(latticeExtra)
library(ggplot2)
library(RColorBrewer)Climate Summaries
Estimates from a single point per map unit polygon
Annual and monthly climate summaries have been estimated from the SSR2 standard stack of 1981–2010 PRISM data. For now, percentiles are estimated from a single sample from within each map unit polygon, weighted by \(log(polygon\ area * component\ percentage)\).
soils <- c('argonaut', 'pierre', 'zook', 'cecil')
s <- fetchOSD(soils, extended = TRUE)| series | climate_var | minimum | q01 | q05 | q25 | q50 | q75 | q95 | q99 | maximum | n |
|---|---|---|---|---|---|---|---|---|---|---|---|
| ARGONAUT | Elevation (m) | 4.00 | 60.00 | 76.00 | 139.00 | 261.00 | 437.00 | 760.00 | 1034.00 | 1271.00 | 842 |
| ARGONAUT | Effective Precipitation (mm) | -368.26 | -284.42 | -268.45 | -132.91 | -16.34 | 77.01 | 290.86 | 460.90 | 747.55 | 842 |
| ARGONAUT | Frost-Free Days | 184.00 | 211.08 | 241.00 | 277.00 | 312.00 | 328.00 | 346.00 | 354.00 | 365.00 | 842 |
| ARGONAUT | Mean Annual Air Temperature (degrees C) | 12.33 | 13.76 | 14.82 | 15.52 | 16.15 | 16.48 | 16.66 | 16.88 | 17.00 | 842 |
| ARGONAUT | Mean Annual Precipitation (mm) | 498.00 | 568.00 | 580.00 | 730.00 | 830.00 | 905.00 | 1094.00 | 1238.00 | 1509.00 | 842 |
| ARGONAUT | Growing Degree Days (degrees C) | 1894.00 | 2148.45 | 2318.00 | 2431.00 | 2544.00 | 2599.00 | 2627.00 | 2674.00 | 2745.00 | 842 |
| ARGONAUT | Fraction of Annual PPT as Rain | 93.00 | 96.00 | 97.00 | 99.00 | 99.00 | 99.00 | 99.00 | 99.00 | 100.00 | 842 |
| ARGONAUT | Design Freeze Index (degrees C) | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 2.00 | 6.00 | 11.00 | 35.00 | 842 |
| series | climate_var | minimum | q01 | q05 | q25 | q50 | q75 | q95 | q99 | maximum | n | month | variable | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 9 | ARGONAUT | ppt1 | 94.00 | 101.00 | 104.00 | 126.00 | 144.00 | 153.00 | 182.00 | 223.00 | 250.00 | 842 | 1 | Precipitation (mm) |
| 10 | ARGONAUT | ppt2 | 90.00 | 98.00 | 101.00 | 126.00 | 144.00 | 156.00 | 181.00 | 210.00 | 253.00 | 842 | 2 | Precipitation (mm) |
| 11 | ARGONAUT | ppt3 | 74.00 | 88.00 | 89.00 | 107.00 | 126.00 | 140.00 | 165.00 | 190.00 | 224.00 | 842 | 3 | Precipitation (mm) |
| 21 | ARGONAUT | pet1 | 13.36 | 13.65 | 13.88 | 14.40 | 14.72 | 14.91 | 15.23 | 15.47 | 25.48 | 842 | 1 | Potential ET (mm) |
| 22 | ARGONAUT | pet2 | 14.51 | 16.53 | 18.30 | 20.65 | 21.27 | 21.77 | 22.36 | 22.83 | 26.54 | 842 | 2 | Potential ET (mm) |
| 23 | ARGONAUT | pet3 | 25.81 | 27.22 | 29.42 | 33.81 | 35.27 | 37.03 | 38.13 | 38.65 | 39.94 | 842 | 3 | Potential ET (mm) |
Visualization of the annual climate summaries using select percentiles.
# control color like this
trellis.par.set(plot.line=list(col='RoyalBlue'))
# control centers symbol and size here
res <- vizAnnualClimate(s$climate.annual, IQR.cex = 1.1, cex=1.1, pch=18)
print(res$fig)
One possible depiction of monthly PPT and PET.
# reasonable colors for a couple of groups
cols <- brewer.pal(9, 'Set1')
cols <- cols[c(1:5,7,9)]
ggplot(s$climate.monthly, aes(x = month, group=series)) +
geom_ribbon(aes(ymin = q25, ymax = q75, fill=series)) +
geom_line(aes(month, q25)) +
geom_line(aes(month, q75)) +
geom_abline(intercept=0, slope = 0, lty=2) +
xlab('') + ylab('mm') +
ggtitle('Monthly IQR') +
scale_fill_manual(values=alpha(cols, 0.75)) +
facet_wrap(vars(variable), scales = 'free_y') +
theme_bw()
Estimates from SoilWeb series extent boundary
# the result is a SpatialPolygonsDataFrame object
see.boundary <- seriesExtent('amador')
# collect as many samples as are present in the fetchOSD summary
see.samples <- spsample(see.boundary, n = 394, type='regular')
# quick map of samples
par(mar=c(0,0,0,0))
map(database='county', regions='california')
points(see.samples, col='royalblue', cex=0.25)
mtext('SEE boundary, 394 samples', side = 1, line = -1.25, at = -124.25, adj = 0)
box()
Sample raster stack and compute select percentiles
Estimates from sampling original map unit polygons
# get polygons associated with map units that contain "amador" as a major component
q <- "select G.MupolygonWktWgs84 as geom, mapunit.mukey, muname
FROM mapunit
CROSS APPLY SDA_Get_MupolygonWktWgs84_from_Mukey(mapunit.mukey) as G
WHERE mukey IN (SELECT DISTINCT mukey FROM component WHERE compname like 'amador%' AND majcompflag = 'Yes')"
# result is a data.frame, "MupolygonWktWgs84" contains WKT representation of geometry
res <- SDA_query(q)
# convert to SPDF
mu.poly <- processSDA_WKT(res)
# map
par(mar=c(0,0,0,0))
map(database='county', regions='california')
plot(mu.poly, add=TRUE, border='royalblue', col='royalblue')
mtext('MU polygons', side = 1, line = -1.25, at = -124.25, adj = 0)
box()
This document is based on aqp version 1.17.01, soilDB version 2.3.6, and sharpshootR version 1.4.01.