In phytomosaic/wildstew: Wilderness Stewardship Performance Reporting
knitr::opts_chunk$set(echo = FALSE)
options(knitr.kable.NA = '')
has_ttest <- FALSE
### set graphical parameters
`set_par` <- function (panels = 1, pty = 's', ...) {
mgp <- c(2, 0.4, 0)
mar <- c(4, 4, 0.5, 0.5)
oma <- c(0, 0, 0, 0)
auto_rowcol <- function(n = panels) {
if (n <= 3)
c(1, n)
else if (n <= 6)
c(2, (n + 1)%/%2)
else if (n <= 12)
c(3, (n + 2)%/%3)
else c(ceiling(n/(nr <- ceiling(sqrt(n)))), nr)
}
mfrow <- auto_rowcol()
if (panels > 4)
panels <- 4
switch(as.character(panels), `4` = par(mfrow = mfrow,
mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L',
las = 1, cex.lab = 1.2, tcl = -0.2, ...), `3` = par(mfrow = mfrow,
mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L',
las = 1, cex.lab = 1.4, cex.axis = 1.2, tcl = -0.2, ...),
`2` = par(mfrow = mfrow, mgp = mgp, mar = mar,
pty = pty, oma = oma, bty = 'L', las = 1, cex.axis = 0.85,
tcl = -0.2, ...), `1` = par(mfrow = mfrow,
mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L',
las = 1, cex.axis = 0.85, tcl = -0.2, ...))
}
### get years of sampling
`get_years` <- function (pick, ...) {
rng <- range(d[d$wa==pick, 'year'], na.rm=TRUE) ### <<--------------------
if (any(!is.finite(rng))) c(NA,NA) else rng
}
### get recent decade
`get_recent_decade` <- function (pick, ...) {
s <- sort(unique(a[a$warea==pick, 'tenyr']))
s <- s[length(s)]
if (any(!is.finite(s))) c(NA) else s
}
### get number of plots
`get_n` <- function (pick, ...) {
n <- unique(a[a$warea==pick, 'n_plots'])
if (any(!is.finite(n))) c(NA) else n
}
### table function with 5-y mean
`get_tab` <- function (pick, xvar, col.names, caption, format){
tab <- a[which(a$warea==pick), xvar]
if (missing(col.names)) col.names <- xvar
if (missing(caption)) caption <- NULL
if (missing(format)) format <- 'pandoc'
j <- xvar[xvar != 'tenyr'] # dont consider year column
tab <- tab[rowSums(is.na(tab[,j])) != ncol(tab[,j]),]
tab[tab %in% c('NA','***NA***','NaN','***NaN***')] <- NA
knitr::kable(x = tab,
format = format,
row.names = F,
align = c('l', rep('r',NCOL(tab)-1)),
longtable = T,
booktabs = T,
caption = caption,
col.names = col.names)
}
\vspace*{\fill}
\begin{wrapfigure}{l}{1.8cm}
\includegraphics{usfs.png}
\end{wrapfigure}
\
Air Resource Management Program
USDA Forest Service
1400 Independence Ave SW
Washington, DC 20250
\
\newpage
Overview
Motivation
Each Wilderness area managed by the USDA Forest Service must report annually on progress toward “wildernesses meeting baseline performance for preserving wilderness character”, otherwise known as Wilderness Stewardship Performance (WSP). Among other elements, Wilderness Air Quality Values (WAQVs) are a marker of natural quality of Wilderness character (WSP Guidebook 2020). This report streamlines annual WSP reporting for WAQVs from lichen-based Critical Loads (CLs).
Sensitive receptors
WAQVs are supported by monitoring changes in “sensitive receptors” like epiphytic macrolichens. Lichens are impacted by air pollutants and therefore indicate changing conditions. Workers long ago established links between lichen tissue content and the atmospheric deposition of pollutants like oxidized nitrogen (NO$_x$), reduced nitrogen (NH$_x$) and oxidized sulfur (SO$_x$). Each species also has a unique and individualistic pollution tolerance, so observing collective responses is far more robust than any single species alone.
Sensitive receptor indicators
For WSP reporting of air quality, “sensitive receptor indicators” can be estimated by the presence and abundance of lichen species. Nationally-consistent FHM/FIA surveys of lichen communities let us score and classify species by their relative pollution tolerance, and to calculate community metrics across all species at a site. Community metrics are more robust than tracking single species because the absence of a single species offers only circumstantial evidence of air quality (what if the species simply wasn't able to disperse to or tolerate microclimate at that site?). The joint absence of multiple pollution-intolerant species provides much stronger evidence of air-quality impacts.
Indicator: Air scores
Individual species' sensitivities are the value of N or S deposition (kg ha^-1^ y^-1^) at which the species was most frequently detected. Air scores for each site are the community-weighted mean of species' sensitivities, for all species present (kg ha^-1^ y^-1^). Therefore, site air scores represent the collective pollution tolerance of co-occurring species: "clean" sites have fairly low air scores driven by pollution-intolerant species, while "polluted" sites have higher air scores driven by pollution-tolerant ones. Air scores are reported as both raw scores (without regard to climate) and as climate-adjusted scores (which permit nationwide comparison).
Indicator: Critical loads exceedances
A critical load (CL) is the air pollutant deposition rate expected to result in harmful ecosystem changes. Exceedance is the degree to which locally observed deposition goes beyond an established CL. Lichen-based CLs were established in prior work (Geiser et al. 2019). Specifically, lichen-based CLs are the deposition rates at which we expect changes in 1) total species richness, 2) sensitive species richness, 3) forage lichen abundance, 4) cyanolichen abundance, and 5) community compositions (see critical loads we used).
Conservation status
Conservation rankings follow NatureServe (https://www.natureserve.org/). Species considered 'at risk' here are those ranked by NatureServe as Critically Imperiled, Imperiled, or Vulnerable, and/or those ranked by COSEWIC as Endangered, Special Concern, or Threatened.
Critical loads we used
We used critical loads from Table 3 in Geiser et al. (2019). These critical loads are invariant to geography and so apply nationwide. Critical loads in this table are deposition levels resulting in $\geq20\%$ decrease in the richness or abundance of the given lichen metric.
\
tab <- data.frame(Metric = c('Species richness',
'Sensitive species richness',
'Forage lichen abundance',
'Cyanolichen abundance',
'Comm composition shifts'),
Ndep = c(3.5, 3.1, 1.9, 1.3, 1.5),
Sdep = c(6.0, 2.5, 2.6, 2.3, 2.7))
knitr::kable(tab, format='pandoc', row.names = F, longtable = T,
booktabs = T,
caption = 'Critical loads used in this report (kg ha^-1^ y^-1^).',
escape = F,
col.names = c('Lichen metric', 'N deposition', 'S deposition'))
See the article:
Geiser, L., P. Nelson, S. Jovan, H. Root, and C. Clark. 2019. Assessing ecological risks from atmospheric deposition of nitrogen and sulfur to US forests using epiphytic macrolichens. Diversity 11 (87):1-29. doi:10.3390/d11060087{target="_blank"}
\
Why monitor lichens? Ecological functions
\begin{wrapfigure}{r}{5.0cm}
\includegraphics{bothlichens.png}
\end{wrapfigure}
Lichens are commonly divided into "functional groups" based on ecological functions: three common groups are forage lichens, cyanolichens and matrix lichens.
Forage lichens (including "old-man's beard" and "witches' hair" lichens) provide food, nesting materials and suitable habitat for many animals (from insects to amphibians to small mammals). Large ungulates like deer, elk or caribou may depend on forage lichens to eat during lean winter months. Loss of these pollution-sensitive lichens might have implications for the wildlife that depend on them.
Cyanolichens (including "lungwort" lichens) contain cyanobacterial partners that convert nitrogen from the air into a form which forest vegetation can absorb -- when cyanolichens fall to the forest floor and release this nitrogen, they serve as "forest fertilizer"! Cyanolichens are among the most intolerant of poor air quality, so their absence may indicate deteriorating conditions.
Matrix lichens are any not included as forage or cyanolichens. They are often more pollution-tolerant than those two groups, and perform an array of ecosystem functions. These functions include rainfall interception, water cycle buffering, animal habitat, and bioindication of forest health.
Paying attention to changes in lichen communities is informative -- losses of lichen functional groups can indicate declines in the ecosystem functions that they provide.
\newpage
Summary for r input$x
rng <- get_years(pick = rx())
n <- get_n(pick = rx())
dec <- get_recent_decade(pick = rx())
Number of surveyed plots: r n
Field sampling dates: r rng[1] to r rng[2]
Most recent ten-year sampling period: r dec
Nitrogen snapshot
### plot the valuebar for observed CMAQ deposition
`plot_valuebar` <- function(pick, brk=b, select='n', ...){
b <- brk
x <- a[which(a$warea==pick),
paste0('cmaq_',select,'_3yroll.mean')]
arrlab <- 'OBSERVED\nCMAQ deposition'
if(all(is.na(x))) { # try lichen airscores if CMAQ is missing
x <- a[which(a$warea==pick),
paste0(select,'_airscore.mean')]
arrlab <- 'OBSERVED\nlichen airscores'
}
if(all(is.na(x))) { # no more tries
x <- 999
}
obs <- x[length(x)] # take most recent one
par(mar=c(4,0.3,0.3,0.3), las=1, tcl=-0.2, mgp=c(1.8,0.4,0))
n <- 100
xseq <- seq(-0.5, 20.0, len=n)
plot(xseq, seq(1, 10, len=n), type='n', bty='n',
yaxt='n', ylab='', ylim=c(1,13.4), yaxs='i', xaxs='r', ...)
cc <- viridis::inferno(5,begin=0.1,end=0.9,alpha=seq(0.7,1,len=5))
rect(-9, -10, b[1], 20, col=cc[1], border=NA)
rect(b[1], -10, b[2], 20, col=cc[2], border=NA)
rect(b[2], -10, b[3], 20, col=cc[3], border=NA)
rect(b[3], -10, b[4], 20, col=cc[4], border=NA)
rect(b[4], -10, 99, 20, col=cc[5], border=NA)
text(x=b[1] - 0.5*(b[1]-(-1)), y=13.1, 'Acceptable', pos=1, cex=0.6)
text(x=b[2] - 0.5*(b[2]-b[1]), y=13.1, 'Marginal', pos=1, cex=0.6)
text(x=b[3] - 0.5*(b[3]-b[2]), y=13.1, 'Concerning', pos=1, cex=0.6)
text(x=b[4] - 0.5*(b[4]-b[3]), y=13.1, 'Poor', pos=1, cex=0.6)
text(x=21 - 0.5*(21-b[4]), y=13.1, 'Very poor', pos=1, cex=0.6)
if (obs >= 999) {
text(x=10, y=10, 'No data available...', pos=1, cex=0.7, font=2)
} else if (obs > 20 & obs < 999) {
arrows(20.2, 6.5, 20.2, 1, length=0.20, angle=15, lwd=3, col=1)
text(x=18.8, y=10, arrlab, pos=1, cex=0.7, font=2)
} else {
arrows(obs, 6.5, obs, 1, length=0.20, angle=15, lwd=3, col=1)
text(x=obs, y=10, arrlab, pos=1, cex=0.7, font=2)
}
}
b <- c(2.5, 6, 9.9, 14) # loosely based on 20/50/80
plot_valuebar(pick = rx(), brk=b, select='n',
xlab=expression(N~deposition~'('*kg~ha^-1~yr^-1*')'))
\
Sulfur snapshot
### plot the valuebar for observed CMAQ S deposition
b <- b + (b * 0.05) # new (loosely based on 20/50/80)
plot_valuebar(pick = rx(), brk=b, select='s',
xlab=expression(S~deposition~'('*kg~ha^-1~yr^-1*')'))
\newpage
x <- d[d$wa==rx(), ]
u <- sort(unique(na.omit(x[,'tenyr']))) # get unique time periods
len <- length(u)
if(len < 2) {
has_ttest <- FALSE
} else {
u <- u[(len-1):(len)] # take *most recent* two periods
x <- droplevels(x[x$tenyr %in% u,,drop=T])
x <- x[order(x$tenyr), ]
### either tenyr have *ALL* NA? (ttest still works if *SOME* NA)
`f` <- function(k) all(is.na(k))
### control logic to select CMAQ > AIRSCORE > NONE (nitrogen)
if (!any(tapply(x$cmaq_n_3yroll, list(x$tenyr), f))) {
fmla_n <- as.formula(cmaq_n_3yroll ~ tenyr)
} else if (!any(tapply(x$n_airscore, list(x$tenyr), f))) {
fmla_n <- as.formula(n_airscore ~ tenyr)
} else {
fmla_n <- NULL
}
### control logic to select CMAQ > AIRSCORE > NONE (sulfur)
if (!any(tapply(x$cmaq_s_3yroll, list(x$tenyr), f))) {
fmla_s <- as.formula(cmaq_s_3yroll ~ tenyr)
} else if (!any(tapply(x$s_airscore, list(x$tenyr), f))) {
fmla_s <- as.formula(s_airscore ~ tenyr)
} else {
fmla_s <- NULL
}
### ensure minimum sample size
`get_n_per_grp` <- function (x, na.rm = TRUE) {
n <- length(x) - sum(is.na(x))
return(c(n=n))
}
n_n <- tryCatch(aggregate(fmla_n, data=x, get_n_per_grp),
error=function(e) return(NULL))
n_s <- tryCatch(aggregate(fmla_s, data=x, get_n_per_grp),
error=function(e) return(NULL))
### control logic, if both are NULL or too-small sample-size,
### then dont do the ttests
if ((is.null(fmla_n) & is.null(fmla_s)) |
(any(n_n[,2] < 2) & any(n_s[,2] < 2))) {
has_ttest <- FALSE
} else {
### UNPAIRED 2-sample ttest (plots NOT paired...)
t1 <- tryCatch(t.test(fmla_n,data=x),
error=function(e)return(NULL))
t3 <- tryCatch(t.test(fmla_s,data=x),
error=function(e)return(NULL))
### descriptive summary table
`f` <- function (x, na.rm = TRUE) {
if (!is.numeric(x)) return(NULL)
m <- mean(x, na.rm = na.rm)
s <- stats::sd(x, na.rm = na.rm)
na <- sum(is.na(x))
n <- length(x) - na
se <- s/sqrt(n - 1)
ci <- 1.96 * se
return(c(mean=round(m,3), ci=round(ci,3), n=n))
}
tab1 <- tryCatch(do.call(data.frame,
aggregate(fmla_n, data=x, f)),
error=function(e)return(NULL))
tab3 <- tryCatch(do.call(data.frame,
aggregate(fmla_s, data=x, f)),
error=function(e)return(NULL))
has_ttest <- TRUE
### if ttests failed then cancel ttests
if (is.null(t1) & is.null(t3)) {
has_ttest <- FALSE
}
}
}
r if(has_ttest) {paste0('## Significant trends\n\nFrom the decade ', u[1], ' to ', u[2], ', a two-sample _t_-test revealed ', ifelse(t1$p.value < 0.05,'a ','NO '), 'SIGNIFICANT trend in N deposition, and ', ifelse(t3$p.value < 0.05,'a ','NO '), 'SIGNIFICANT trend in S deposition. ', 'See figure for test results.')}
### PLOT
`plot_ci` <- function (tab, mod, ylab, ...) {
nm <- tab[,1]
y <- tab[,2]
yci <- tab[,3]
n <- tab[,4]
x <- 1:2
ylm <- c(0, # min(c(y, y - yci), na.rm=T) - 0.4,
max(c(y, y + yci), na.rm=T) + 0.1)
if (all(is.infinite(ylm))) return(NULL)
if (anyNA(ylm)) { ylm[which(is.na(ylm))] <- NULL }
if (missing(ylab)) {
ylb <- expression(Air~score~(kg~ha^'-1'~y^'-1'))
} else {
ylb <- ylab
}
plot(x=x, y=tab[,2], pch=16, ylab=ylb, xlab='', xaxt='n',
xlim=c(0.4,2.6), ylim=ylm, ...)
xlb <- paste0(nm,'\n','n = ', n)
axis(1, at=x, labels=xlb, tick=F, line=1, cex.axis=0.9)
w <- 0.05
segments(x0=x, x1=x, y0=y-yci, y1=y+yci, lwd=1.75)
segments(x0=x-w, x1=x+w, y0=y+yci, y1=y+yci, lwd=1.75)
segments(x0=x-w, x1=x+w, y0=y-yci, y1=y-yci, lwd=1.75)
`add_text` <- function (x, y, labels, ...) {
text(x = graphics::grconvertX(x, from = "npc", to = "user"),
y = graphics::grconvertY(y, from = "npc", to = "user"),
labels = labels, adj = 0, ...)
}
add_text(0.05,0.19,
lab=paste0('t = ', formatC(round(mod$statistic,1),digits=3)),
cex=0.8)
add_text(0.05,0.12,
lab=ifelse(mod$p.value < 0.01, 'p < 0.01',
paste0('p = ',formatC(round(mod$p.value,2),digits=4))),
cex=0.8)
add_text(0.05,0.05,
lab=paste0('Est. diff. = ',
formatC(round(diff(mod$estimate),1),digits=2)),
cex=0.8)
}
### setup labeling
if (grepl('cmaq', names(tab1)[2])) {
ylb_n <- expression(CMAQ~value~(kg~N~ha^'-1'~y^'-1'))
} else if (grepl('airscore', names(tab1)[2])) {
ylb_n <- expression(N~air~score~(kg~N~ha^'-1'~y^'-1'))
} else {
ylb_n <- NA
}
if (grepl('cmaq', names(tab3)[2])) {
ylb_s <- expression(CMAQ~value~(kg~S~ha^'-1'~y^'-1'))
} else if (grepl('airscore', names(tab3)[2])) {
ylb_s <- expression(S~air~score~(kg~S~ha^'-1'~y^'-1'))
} else {
ylb_s <- NA
}
### plot
set_par(2)
par(mgp=c(1.3,0.4,0))
try(plot_ci(tab1, t1, ylab=ylb_n), silent=TRUE)
try(plot_ci(tab3, t3, ylab=ylb_s), silent=TRUE)
rm(get_n_per_grp, plot_ci, len, x, u, f, t1, t3, tab1, tab3, ylb_n, ylb_s, fmla_n, fmla_s)
r if(has_ttest) {"Power analysis: a _t_-test would require _n_ = 17 plots per time interval to detect a change of 1 kg ha^-1^ y^-1^ with 80\\% probability.\\newpage"}
Area map
### map wilderness boundary and sampling locations
`map_points` <- function (pick=pick, pts=d, shp=w, ras=dem,
lwd.cont=1, ...) {
aea_prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5
+lon_0=-110 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m"
trg_prj <- projection(shp)
shp <- shp[shp$NAME == pick,] # select your wilderness
pts <- pts[pts$wa == pick,] # select your Wilderness
if(all(is.na(pts$lon)) | all(is.na(pts$lon))) {
e <- extent(shp) * 1.25
has_points <- FALSE
} else {
has_points <- TRUE
`make_xy` <- function(xy, crs, ...){
xy <- as.data.frame(xy)
xy <- xy[!(is.na(xy$lon) | is.na(xy$lat)),] # rm NAs
coordinates(xy) <- ~lon+lat
proj4string(xy) <- CRS('+init=epsg:4269') # NAD83 dec degree
return(spTransform(xy, crs))
}
ppts <- make_xy(pts, crs=trg_prj)
near <- c(rgeos::gDistance(spTransform(ppts, CRS(aea_prj)),
spTransform(shp, CRS(aea_prj)),
byid=T) < 3000) # reject pts > 3 km
ppts <- ppts[near,]
e <- cbind(matrix(extent(shp)),matrix(extent(ppts)))
e <- c(min(e[1,]),max(e[2,]),min(e[3,]),max(e[4,])) # combined ext
`haversine` <- function(loc1, loc2) { # geographic distance
lat1 <- loc1[1] * pi / 180
lon1 <- loc1[2] * pi / 180
lat2 <- loc2[1] * pi / 180
lon2 <- loc2[2] * pi / 180
R <- 6371.0 # avg earth radius [km]
dlat <- lat2 - lat1
dlon <- lon2 - lon1
a <- sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2
c <- 2 * atan2(sqrt(a), sqrt(1-a))
return(R * c)
}
dy <- haversine(c(e[3],e[1]), c(e[4],e[1])) # kilometers y-axis
dx <- haversine(c(e[3],e[1]), c(e[3],e[2])) # kilometers x-axis
e <- extent(e)
if (dy > dx) {
e[1:2] <- matrix(e * dy/dx)[1:2,] # expand x
} else {
e[3:4] <- matrix(e * dx/dy)[3:4,] # expand y
}
dx <- as.numeric( # how large a scale bar to make?
as.character(
cut(dx,breaks=c(0,5,50,200,Inf),labels=c(1,5,10,50),incl=T)
)
)
e <- extent(e) * 1.25 # slight expansion
}
ras <- disaggregate(crop(ras, e), 10, method='bilinear')
### custom function for filled contours....
`fillcont` <- function (x, maxpixels = 1e+05, ...) {
`fc` <- function (
x = seq(0, 1, length.out = nrow(z)),
y = seq(0, 1, length.out = ncol(z)),
z,
xlim = range(x, finite = TRUE),
ylim = range(y, finite = TRUE),
zlim = range(z, finite = TRUE),
levels = pretty(zlim, nlevels),
nlevels = 20,
color.palette = function(n) viridis::viridis(n,begin=0,end=0.9),
col = color.palette(length(levels) - 1), ...) {
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
}
else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
}
else stop("no 'z' matrix specified")
}
else if (is.list(x)) {
y <- x$y
x <- x$x
}
if (any(diff(x) <= 0) || any(diff(y) <= 0))
stop("increasing 'x' and 'y' values expected")
plot.new()
plot.window(xlim, ylim, "", xaxs = "i", yaxs = "i", asp = NA)
.filled.contour(x, y, z, levels, col)
}
x <- sampleRegular(x, maxpixels, asRaster = T, useGDAL = T)
X <- xFromCol(x, 1:ncol(x))
Y <- yFromRow(x, nrow(x):1)
Z <- t(matrix(getValues(x), ncol = x@ncols, byrow=T)[nrow(x):1,])
lonlat <- couldBeLonLat(x, warnings = FALSE)
asp <- list(...)$asp
if (is.null(asp)) {
if (lonlat) {
ym <- mean(c(x@extent@ymax, x@extent@ymin))
asp <- 1/cos((ym * pi)/180)
}
else {
asp <- 1
}
fc(x = X, y = Y, z = Z, asp = asp, ...)
}
else {
fc(x = X, y = Y, z = Z, ...)
}
}
fillcont(ras, nlevels=20, lwd=lwd.cont)
plot(shp, col='transparent', border='red', add=T, lwd=2)
if (has_points) {
points(ppts, pch=16, col='#000000', cex=1)
scalebar(d=dx, xy=c(e[1] + abs(diff(e[1:2]))*0.05,
e[3] + abs(diff(e[3:4]))*0.08),
type='bar', divs=2, below=paste0('\n', dx, ' km'),
adj=c(0.5,-1), lab='', lwd=2, cex=0.7, font=2)
}
}
### some wilderness areas too large to plot (memory overrun on server)
too_big <- c(
'Oregon Islands Wilderness',
'Frank Church-River of No Return Wilderness',
'Kootznoowoo Wilderness',
'Boundary Waters Canoe Area Wilderness'
# 'Daniel J. Evans Wilderness' # <<--- AKA Olympic NP, should be ok
)
# error catch -- dont plot if shape too large to plot (server memory overrun)
if (rx() %in% too_big) {
plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE)
text(5,5,labels='Map not available',adj=c(0.5,0.5))
} else { # error catch -- dont plot if valid name not in shapefile
tryCatch(map_points(pick = rx(), pts=d, shp=w, ras=dem),
error = function(e) {
plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE)
text(5,5,labels='Map not available',adj=c(0.5,0.5))
})
}
# ### error catch -- dont plot if valid name not in shapefile
# `try_map` <- function() {
# tryCatch(map_points(pick = rx(), pts=d, shp=w, ras=dem),
# error = function(e) {
# plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE)
# text(5,5,labels='Map not available',adj=c(0.5,0.5))
# })
# }
# try_map()
# # map_points(pick = rx(), pts=d, shp=w, ras=dem)
\newpage
Site scores
You can think of site scores (or "air scores") as the "average" position of species along an air quality gradient. Technically, these are the community-weighted mean of species' sensitivities, where sensitivities are defined as the value of N or S deposition (kg ha^-1^ y^-1^) at which each species was most frequently detected. We report both raw and climate-adjusted scores, and their nationwide ranking as a percentile.^[See "Indicator: Air scores" section for more information.].
Nitrogen
N: mean deposition
get_tab(pick = rx(),
xvar=c('tenyr','n_airscore.mean','nadj.mean','cmaq_n_3yroll.mean'),
col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'),
caption= 'Site scores for estimated nitrogen deposition (N), as a *mean* (kg ha^-1^ y^-1^).')
N: percentiles
get_tab(pick = rx(),
xvar=c('tenyr','p_n_airscore.mean','p_n_adj.mean','p_n_cmaq.mean'),
col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'),
caption= 'Site scores for estimated nitrogen deposition (N), as a *percentile* (0-100).')
Sulfur
S: mean deposition
get_tab(pick = rx(),
xvar=c('tenyr','s_airscore.mean','sadj.mean','cmaq_s_3yroll.mean'),
col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'),
caption= 'Site scores for estimated sulfur deposition (S), as a *mean* (kg ha^-1^ y^-1^).')
S: percentiles
get_tab(pick = rx(),
xvar=c('tenyr','p_s_airscore.mean','p_s_adj.mean','p_s_cmaq.mean'),
col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'),
caption= 'Site scores for estimated sulfur deposition (S), as a *percentile* (0-100).')
\newpage
Critical loads exceedances
Exceedances are the departure of raw or climate-adjusted lichen airscores from established lichen-based critical loads^[See "Critical loads we used" section.], expressed as the mean of a ten-year period. Negative values indicate the "safety margin" below the critical load threshold, positive values indicate exceedance above the threshold.
Nitrogen
N exceedance: raw lichen airscores
get_tab(pick = rx(),
xvar=c('tenyr', 'n_raw_ex_sr.mean', 'n_raw_ex_sr_sens.mean',
'n_raw_ex_sr_forag.mean', 'n_raw_ex_sr_cyano.mean',
'n_raw_ex_compn.mean'),
col.names=c('Ten-year interval', 'Richness', 'Sensitive richness',
'Forage lichen abund', 'Cyanolichen abund',
'Comm comp'),
caption= 'CL exceedances for mean nitrogen deposition (kg N ha^-1^ y^-1^), based on raw lichen airscores.')
N exceedance: climate-adjusted lichen airscores
get_tab(pick = rx(),
xvar=c('tenyr', 'n_adj_ex_sr.mean',
'n_adj_ex_sr_sens.mean',
'n_adj_ex_sr_forag.mean',
'n_adj_ex_sr_cyano.mean',
'n_adj_ex_compn.mean'),
col.names=c('Ten-year interval', 'Richness', 'Sensitive richness',
'Forage lichen abund', 'Cyanolichen abund',
'Comm comp'),
caption= 'CL exceedances for mean nitrogen deposition (kg N ha^-1^ y^-1^), based on climate-adjusted lichen airscores.')
Sulfur
S exceedance: raw lichen airscores
get_tab(pick = rx(),
xvar=c('tenyr', 's_raw_ex_sr.mean', 's_raw_ex_sr_sens.mean',
's_raw_ex_sr_forag.mean', 's_raw_ex_sr_cyano.mean',
's_raw_ex_compn.mean'),
col.names=c('Ten-year interval', 'Richness', 'Sensitive richness',
'Forage lichen abund', 'Cyanolichen abund',
'Comm comp'),
caption= 'CL exceedances for mean sulfur deposition (kg S ha^-1^ y^-1^), based on raw lichen airscores.')
S exceedance: climate-adjusted lichen airscores
get_tab(pick = rx(),
xvar=c('tenyr', 's_adj_ex_sr.mean', 's_adj_ex_sr_sens.mean',
's_adj_ex_sr_forag.mean','s_adj_ex_sr_cyano.mean',
's_adj_ex_compn.mean'),
col.names=c('Ten-year interval', 'Richness', 'Sensitive richness',
'Forage lichen abund', 'Cyanolichen abund',
'Comm comp'),
caption= 'CL exceedances for mean sulfur deposition (kg S ha^-1^ y^-1^), based on climate-adjusted lichen airscores.')
\newpage
Lichen species observed
Sampling completeness
`get_richness` <- function (pick, ...) {
NROW(s[[which(names(s) == pick)]])
}
spp_obs <- get_richness(pick = rx())
spp_exp <- 4.01 + 25.81*log(n) # log-linear eqn
Number of surveyed plots: r n
Number of species observed: r spp_obs
Number of species expected: r round(spp_exp, 1)
\
x <- 1:2500 # n plots sampled
y <- 4.01 + 25.81*log(x) # n species expected
set_par(1)
plot(x, y, type='l', log='x', xaxt='n', xlim=c(1,2500),
xlab=bquote(Plots~sampled~(log[10]~scale)),
ylab='Species expected')
atx <- outer(1:9, 10^(0:3))
labx <- ifelse(log10(atx) %% 1 == 0, atx, NA)
axis(1, at=atx, labels=labx, las=1)
lines(rep(n,2), c(spp_exp,spp_obs))
points(n, spp_exp, pch=21, cex=1.2, bg='white')
points(n, spp_obs, pch=16, cex=1.2)
legend('topleft', pch=c(16,21), bty='n', cex=0.7, pt.cex=1.0,
legend=c(paste0('Observed = ', spp_obs),
paste0('Expected = ', round(spp_exp,1))))
Rare and pollution-sensitive species
`get_natureserve` <- function (pick, ...) {
tab <- s[[which(names(s) == pick)]]
rich <- NROW(tab)
nsens <- round(sum(tab$n_class == 'oligotroph') /
sum(tab$n_class != '') * 100, 1)
ssens <- round(sum(tab$s_class == 'sensitive') /
sum(tab$s_class != '') * 100, 1)
atrisk <- sum(1 * grepl('1|2|3|Endanger|Concern|Threaten',
tab$consrv_stat))
list(rich=rich, nsens=nsens, ssens=ssens, atrisk=atrisk)
}
x <- get_natureserve(pick = rx())
Number of species observed in total: r x$rich
Number of species with 'at risk' conservation status: r x$atrisk
Percentage of rated species with 'sensitive (oligotroph)' N tolerance: r x$nsens\%
Percentage of rated species with 'sensitive' S tolerance: r x$ssens\%
Species list
List of all macrolichens in plots at r input$x from r rng[1] to r rng[2]. True diversity is likely greater because plot footprints are just a fraction of total Wilderness area. Functional groups ('Fxl grp') are cyanolichens ('cyano') and forage lichens ('forage')(see descriptions). Sensitivity for each species ranges from low to high tolerance for nitrogen (oligo-, meso-, eutroph) and sulfur (sensitive, intermediate, tolerant). Critical loads ('CL') are deposition levels resulting in >20\% decrease in occurrence. Conservation from NatureServe (https://www.natureserve.org/) assigns each species both a global (G) and state (S) rank on a scale of 1 to 5 (1=Critically imperiled; 2=Imperiled; 3=Vulnerable; 4=Apparently secure; 5=Secure).
`get_tab_species` <- function (pick, ...) {
tab <- s[which(names(s) == pick)][[1]]
tab <- tab[order(tab$n_critload, tab$s_critload, tab$fxlgp),]
knitr::kable(x = tab,
format = 'markdown', ### format = 'pandoc',
row.names = FALSE,
longtable = TRUE,
booktabs = TRUE,
caption = NULL, ### 'Macrolichen species list.',
escape = FALSE,
col.names = c('Species','Fxl grp','N sensitivity','N CL',
'S sensitivity','S CL','Conservation')
)
}
get_tab_species(pick = rx())
phytomosaic/wildstew documentation built on Feb. 19, 2022, 12:40 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE) options(knitr.kable.NA = '') has_ttest <- FALSE
### set graphical parameters `set_par` <- function (panels = 1, pty = 's', ...) { mgp <- c(2, 0.4, 0) mar <- c(4, 4, 0.5, 0.5) oma <- c(0, 0, 0, 0) auto_rowcol <- function(n = panels) { if (n <= 3) c(1, n) else if (n <= 6) c(2, (n + 1)%/%2) else if (n <= 12) c(3, (n + 2)%/%3) else c(ceiling(n/(nr <- ceiling(sqrt(n)))), nr) } mfrow <- auto_rowcol() if (panels > 4) panels <- 4 switch(as.character(panels), `4` = par(mfrow = mfrow, mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L', las = 1, cex.lab = 1.2, tcl = -0.2, ...), `3` = par(mfrow = mfrow, mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L', las = 1, cex.lab = 1.4, cex.axis = 1.2, tcl = -0.2, ...), `2` = par(mfrow = mfrow, mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L', las = 1, cex.axis = 0.85, tcl = -0.2, ...), `1` = par(mfrow = mfrow, mgp = mgp, mar = mar, pty = pty, oma = oma, bty = 'L', las = 1, cex.axis = 0.85, tcl = -0.2, ...)) } ### get years of sampling `get_years` <- function (pick, ...) { rng <- range(d[d$wa==pick, 'year'], na.rm=TRUE) ### <<-------------------- if (any(!is.finite(rng))) c(NA,NA) else rng } ### get recent decade `get_recent_decade` <- function (pick, ...) { s <- sort(unique(a[a$warea==pick, 'tenyr'])) s <- s[length(s)] if (any(!is.finite(s))) c(NA) else s } ### get number of plots `get_n` <- function (pick, ...) { n <- unique(a[a$warea==pick, 'n_plots']) if (any(!is.finite(n))) c(NA) else n } ### table function with 5-y mean `get_tab` <- function (pick, xvar, col.names, caption, format){ tab <- a[which(a$warea==pick), xvar] if (missing(col.names)) col.names <- xvar if (missing(caption)) caption <- NULL if (missing(format)) format <- 'pandoc' j <- xvar[xvar != 'tenyr'] # dont consider year column tab <- tab[rowSums(is.na(tab[,j])) != ncol(tab[,j]),] tab[tab %in% c('NA','***NA***','NaN','***NaN***')] <- NA knitr::kable(x = tab, format = format, row.names = F, align = c('l', rep('r',NCOL(tab)-1)), longtable = T, booktabs = T, caption = caption, col.names = col.names) }
\vspace*{\fill} \begin{wrapfigure}{l}{1.8cm} \includegraphics{usfs.png} \end{wrapfigure}
\
Air Resource Management Program
USDA Forest Service
1400 Independence Ave SW
Washington, DC 20250
\
\newpage
Overview
Motivation
Each Wilderness area managed by the USDA Forest Service must report annually on progress toward “wildernesses meeting baseline performance for preserving wilderness character”, otherwise known as Wilderness Stewardship Performance (WSP). Among other elements, Wilderness Air Quality Values (WAQVs) are a marker of natural quality of Wilderness character (WSP Guidebook 2020). This report streamlines annual WSP reporting for WAQVs from lichen-based Critical Loads (CLs).
Sensitive receptors
WAQVs are supported by monitoring changes in “sensitive receptors” like epiphytic macrolichens. Lichens are impacted by air pollutants and therefore indicate changing conditions. Workers long ago established links between lichen tissue content and the atmospheric deposition of pollutants like oxidized nitrogen (NO$_x$), reduced nitrogen (NH$_x$) and oxidized sulfur (SO$_x$). Each species also has a unique and individualistic pollution tolerance, so observing collective responses is far more robust than any single species alone.
Sensitive receptor indicators
For WSP reporting of air quality, “sensitive receptor indicators” can be estimated by the presence and abundance of lichen species. Nationally-consistent FHM/FIA surveys of lichen communities let us score and classify species by their relative pollution tolerance, and to calculate community metrics across all species at a site. Community metrics are more robust than tracking single species because the absence of a single species offers only circumstantial evidence of air quality (what if the species simply wasn't able to disperse to or tolerate microclimate at that site?). The joint absence of multiple pollution-intolerant species provides much stronger evidence of air-quality impacts.
Indicator: Air scores
Individual species' sensitivities are the value of N or S deposition (kg ha^-1^ y^-1^) at which the species was most frequently detected. Air scores for each site are the community-weighted mean of species' sensitivities, for all species present (kg ha^-1^ y^-1^). Therefore, site air scores represent the collective pollution tolerance of co-occurring species: "clean" sites have fairly low air scores driven by pollution-intolerant species, while "polluted" sites have higher air scores driven by pollution-tolerant ones. Air scores are reported as both raw scores (without regard to climate) and as climate-adjusted scores (which permit nationwide comparison).
Indicator: Critical loads exceedances
A critical load (CL) is the air pollutant deposition rate expected to result in harmful ecosystem changes. Exceedance is the degree to which locally observed deposition goes beyond an established CL. Lichen-based CLs were established in prior work (Geiser et al. 2019). Specifically, lichen-based CLs are the deposition rates at which we expect changes in 1) total species richness, 2) sensitive species richness, 3) forage lichen abundance, 4) cyanolichen abundance, and 5) community compositions (see critical loads we used).
Conservation status
Conservation rankings follow NatureServe (https://www.natureserve.org/). Species considered 'at risk' here are those ranked by NatureServe as Critically Imperiled, Imperiled, or Vulnerable, and/or those ranked by COSEWIC as Endangered, Special Concern, or Threatened.
Critical loads we used
We used critical loads from Table 3 in Geiser et al. (2019). These critical loads are invariant to geography and so apply nationwide. Critical loads in this table are deposition levels resulting in $\geq20\%$ decrease in the richness or abundance of the given lichen metric.
\
tab <- data.frame(Metric = c('Species richness', 'Sensitive species richness', 'Forage lichen abundance', 'Cyanolichen abundance', 'Comm composition shifts'), Ndep = c(3.5, 3.1, 1.9, 1.3, 1.5), Sdep = c(6.0, 2.5, 2.6, 2.3, 2.7)) knitr::kable(tab, format='pandoc', row.names = F, longtable = T, booktabs = T, caption = 'Critical loads used in this report (kg ha^-1^ y^-1^).', escape = F, col.names = c('Lichen metric', 'N deposition', 'S deposition'))
See the article:
Geiser, L., P. Nelson, S. Jovan, H. Root, and C. Clark. 2019. Assessing ecological risks from atmospheric deposition of nitrogen and sulfur to US forests using epiphytic macrolichens. Diversity 11 (87):1-29. doi:10.3390/d11060087{target="_blank"}
\
Why monitor lichens? Ecological functions
\begin{wrapfigure}{r}{5.0cm} \includegraphics{bothlichens.png} \end{wrapfigure}
Lichens are commonly divided into "functional groups" based on ecological functions: three common groups are forage lichens, cyanolichens and matrix lichens.
Forage lichens (including "old-man's beard" and "witches' hair" lichens) provide food, nesting materials and suitable habitat for many animals (from insects to amphibians to small mammals). Large ungulates like deer, elk or caribou may depend on forage lichens to eat during lean winter months. Loss of these pollution-sensitive lichens might have implications for the wildlife that depend on them.
Cyanolichens (including "lungwort" lichens) contain cyanobacterial partners that convert nitrogen from the air into a form which forest vegetation can absorb -- when cyanolichens fall to the forest floor and release this nitrogen, they serve as "forest fertilizer"! Cyanolichens are among the most intolerant of poor air quality, so their absence may indicate deteriorating conditions.
Matrix lichens are any not included as forage or cyanolichens. They are often more pollution-tolerant than those two groups, and perform an array of ecosystem functions. These functions include rainfall interception, water cycle buffering, animal habitat, and bioindication of forest health.
Paying attention to changes in lichen communities is informative -- losses of lichen functional groups can indicate declines in the ecosystem functions that they provide.
\newpage
Summary for r input$x
rng <- get_years(pick = rx()) n <- get_n(pick = rx()) dec <- get_recent_decade(pick = rx())
Number of surveyed plots: r n
Field sampling dates: r rng[1] to r rng[2]
Most recent ten-year sampling period: r dec
Nitrogen snapshot
### plot the valuebar for observed CMAQ deposition `plot_valuebar` <- function(pick, brk=b, select='n', ...){ b <- brk x <- a[which(a$warea==pick), paste0('cmaq_',select,'_3yroll.mean')] arrlab <- 'OBSERVED\nCMAQ deposition' if(all(is.na(x))) { # try lichen airscores if CMAQ is missing x <- a[which(a$warea==pick), paste0(select,'_airscore.mean')] arrlab <- 'OBSERVED\nlichen airscores' } if(all(is.na(x))) { # no more tries x <- 999 } obs <- x[length(x)] # take most recent one par(mar=c(4,0.3,0.3,0.3), las=1, tcl=-0.2, mgp=c(1.8,0.4,0)) n <- 100 xseq <- seq(-0.5, 20.0, len=n) plot(xseq, seq(1, 10, len=n), type='n', bty='n', yaxt='n', ylab='', ylim=c(1,13.4), yaxs='i', xaxs='r', ...) cc <- viridis::inferno(5,begin=0.1,end=0.9,alpha=seq(0.7,1,len=5)) rect(-9, -10, b[1], 20, col=cc[1], border=NA) rect(b[1], -10, b[2], 20, col=cc[2], border=NA) rect(b[2], -10, b[3], 20, col=cc[3], border=NA) rect(b[3], -10, b[4], 20, col=cc[4], border=NA) rect(b[4], -10, 99, 20, col=cc[5], border=NA) text(x=b[1] - 0.5*(b[1]-(-1)), y=13.1, 'Acceptable', pos=1, cex=0.6) text(x=b[2] - 0.5*(b[2]-b[1]), y=13.1, 'Marginal', pos=1, cex=0.6) text(x=b[3] - 0.5*(b[3]-b[2]), y=13.1, 'Concerning', pos=1, cex=0.6) text(x=b[4] - 0.5*(b[4]-b[3]), y=13.1, 'Poor', pos=1, cex=0.6) text(x=21 - 0.5*(21-b[4]), y=13.1, 'Very poor', pos=1, cex=0.6) if (obs >= 999) { text(x=10, y=10, 'No data available...', pos=1, cex=0.7, font=2) } else if (obs > 20 & obs < 999) { arrows(20.2, 6.5, 20.2, 1, length=0.20, angle=15, lwd=3, col=1) text(x=18.8, y=10, arrlab, pos=1, cex=0.7, font=2) } else { arrows(obs, 6.5, obs, 1, length=0.20, angle=15, lwd=3, col=1) text(x=obs, y=10, arrlab, pos=1, cex=0.7, font=2) } } b <- c(2.5, 6, 9.9, 14) # loosely based on 20/50/80 plot_valuebar(pick = rx(), brk=b, select='n', xlab=expression(N~deposition~'('*kg~ha^-1~yr^-1*')'))
\
Sulfur snapshot
### plot the valuebar for observed CMAQ S deposition b <- b + (b * 0.05) # new (loosely based on 20/50/80) plot_valuebar(pick = rx(), brk=b, select='s', xlab=expression(S~deposition~'('*kg~ha^-1~yr^-1*')'))
\newpage
x <- d[d$wa==rx(), ] u <- sort(unique(na.omit(x[,'tenyr']))) # get unique time periods len <- length(u) if(len < 2) { has_ttest <- FALSE } else { u <- u[(len-1):(len)] # take *most recent* two periods x <- droplevels(x[x$tenyr %in% u,,drop=T]) x <- x[order(x$tenyr), ] ### either tenyr have *ALL* NA? (ttest still works if *SOME* NA) `f` <- function(k) all(is.na(k)) ### control logic to select CMAQ > AIRSCORE > NONE (nitrogen) if (!any(tapply(x$cmaq_n_3yroll, list(x$tenyr), f))) { fmla_n <- as.formula(cmaq_n_3yroll ~ tenyr) } else if (!any(tapply(x$n_airscore, list(x$tenyr), f))) { fmla_n <- as.formula(n_airscore ~ tenyr) } else { fmla_n <- NULL } ### control logic to select CMAQ > AIRSCORE > NONE (sulfur) if (!any(tapply(x$cmaq_s_3yroll, list(x$tenyr), f))) { fmla_s <- as.formula(cmaq_s_3yroll ~ tenyr) } else if (!any(tapply(x$s_airscore, list(x$tenyr), f))) { fmla_s <- as.formula(s_airscore ~ tenyr) } else { fmla_s <- NULL } ### ensure minimum sample size `get_n_per_grp` <- function (x, na.rm = TRUE) { n <- length(x) - sum(is.na(x)) return(c(n=n)) } n_n <- tryCatch(aggregate(fmla_n, data=x, get_n_per_grp), error=function(e) return(NULL)) n_s <- tryCatch(aggregate(fmla_s, data=x, get_n_per_grp), error=function(e) return(NULL)) ### control logic, if both are NULL or too-small sample-size, ### then dont do the ttests if ((is.null(fmla_n) & is.null(fmla_s)) | (any(n_n[,2] < 2) & any(n_s[,2] < 2))) { has_ttest <- FALSE } else { ### UNPAIRED 2-sample ttest (plots NOT paired...) t1 <- tryCatch(t.test(fmla_n,data=x), error=function(e)return(NULL)) t3 <- tryCatch(t.test(fmla_s,data=x), error=function(e)return(NULL)) ### descriptive summary table `f` <- function (x, na.rm = TRUE) { if (!is.numeric(x)) return(NULL) m <- mean(x, na.rm = na.rm) s <- stats::sd(x, na.rm = na.rm) na <- sum(is.na(x)) n <- length(x) - na se <- s/sqrt(n - 1) ci <- 1.96 * se return(c(mean=round(m,3), ci=round(ci,3), n=n)) } tab1 <- tryCatch(do.call(data.frame, aggregate(fmla_n, data=x, f)), error=function(e)return(NULL)) tab3 <- tryCatch(do.call(data.frame, aggregate(fmla_s, data=x, f)), error=function(e)return(NULL)) has_ttest <- TRUE ### if ttests failed then cancel ttests if (is.null(t1) & is.null(t3)) { has_ttest <- FALSE } } }
r if(has_ttest) {paste0('## Significant trends\n\nFrom the decade ', u[1], ' to ', u[2], ', a two-sample _t_-test revealed ', ifelse(t1$p.value < 0.05,'a ','NO '), 'SIGNIFICANT trend in N deposition, and ', ifelse(t3$p.value < 0.05,'a ','NO '), 'SIGNIFICANT trend in S deposition. ', 'See figure for test results.')}
### PLOT `plot_ci` <- function (tab, mod, ylab, ...) { nm <- tab[,1] y <- tab[,2] yci <- tab[,3] n <- tab[,4] x <- 1:2 ylm <- c(0, # min(c(y, y - yci), na.rm=T) - 0.4, max(c(y, y + yci), na.rm=T) + 0.1) if (all(is.infinite(ylm))) return(NULL) if (anyNA(ylm)) { ylm[which(is.na(ylm))] <- NULL } if (missing(ylab)) { ylb <- expression(Air~score~(kg~ha^'-1'~y^'-1')) } else { ylb <- ylab } plot(x=x, y=tab[,2], pch=16, ylab=ylb, xlab='', xaxt='n', xlim=c(0.4,2.6), ylim=ylm, ...) xlb <- paste0(nm,'\n','n = ', n) axis(1, at=x, labels=xlb, tick=F, line=1, cex.axis=0.9) w <- 0.05 segments(x0=x, x1=x, y0=y-yci, y1=y+yci, lwd=1.75) segments(x0=x-w, x1=x+w, y0=y+yci, y1=y+yci, lwd=1.75) segments(x0=x-w, x1=x+w, y0=y-yci, y1=y-yci, lwd=1.75) `add_text` <- function (x, y, labels, ...) { text(x = graphics::grconvertX(x, from = "npc", to = "user"), y = graphics::grconvertY(y, from = "npc", to = "user"), labels = labels, adj = 0, ...) } add_text(0.05,0.19, lab=paste0('t = ', formatC(round(mod$statistic,1),digits=3)), cex=0.8) add_text(0.05,0.12, lab=ifelse(mod$p.value < 0.01, 'p < 0.01', paste0('p = ',formatC(round(mod$p.value,2),digits=4))), cex=0.8) add_text(0.05,0.05, lab=paste0('Est. diff. = ', formatC(round(diff(mod$estimate),1),digits=2)), cex=0.8) } ### setup labeling if (grepl('cmaq', names(tab1)[2])) { ylb_n <- expression(CMAQ~value~(kg~N~ha^'-1'~y^'-1')) } else if (grepl('airscore', names(tab1)[2])) { ylb_n <- expression(N~air~score~(kg~N~ha^'-1'~y^'-1')) } else { ylb_n <- NA } if (grepl('cmaq', names(tab3)[2])) { ylb_s <- expression(CMAQ~value~(kg~S~ha^'-1'~y^'-1')) } else if (grepl('airscore', names(tab3)[2])) { ylb_s <- expression(S~air~score~(kg~S~ha^'-1'~y^'-1')) } else { ylb_s <- NA } ### plot set_par(2) par(mgp=c(1.3,0.4,0)) try(plot_ci(tab1, t1, ylab=ylb_n), silent=TRUE) try(plot_ci(tab3, t3, ylab=ylb_s), silent=TRUE) rm(get_n_per_grp, plot_ci, len, x, u, f, t1, t3, tab1, tab3, ylb_n, ylb_s, fmla_n, fmla_s)
r if(has_ttest) {"Power analysis: a _t_-test would require _n_ = 17 plots per time interval to detect a change of 1 kg ha^-1^ y^-1^ with 80\\% probability.\\newpage"}
Area map
### map wilderness boundary and sampling locations `map_points` <- function (pick=pick, pts=d, shp=w, ras=dem, lwd.cont=1, ...) { aea_prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-110 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m" trg_prj <- projection(shp) shp <- shp[shp$NAME == pick,] # select your wilderness pts <- pts[pts$wa == pick,] # select your Wilderness if(all(is.na(pts$lon)) | all(is.na(pts$lon))) { e <- extent(shp) * 1.25 has_points <- FALSE } else { has_points <- TRUE `make_xy` <- function(xy, crs, ...){ xy <- as.data.frame(xy) xy <- xy[!(is.na(xy$lon) | is.na(xy$lat)),] # rm NAs coordinates(xy) <- ~lon+lat proj4string(xy) <- CRS('+init=epsg:4269') # NAD83 dec degree return(spTransform(xy, crs)) } ppts <- make_xy(pts, crs=trg_prj) near <- c(rgeos::gDistance(spTransform(ppts, CRS(aea_prj)), spTransform(shp, CRS(aea_prj)), byid=T) < 3000) # reject pts > 3 km ppts <- ppts[near,] e <- cbind(matrix(extent(shp)),matrix(extent(ppts))) e <- c(min(e[1,]),max(e[2,]),min(e[3,]),max(e[4,])) # combined ext `haversine` <- function(loc1, loc2) { # geographic distance lat1 <- loc1[1] * pi / 180 lon1 <- loc1[2] * pi / 180 lat2 <- loc2[1] * pi / 180 lon2 <- loc2[2] * pi / 180 R <- 6371.0 # avg earth radius [km] dlat <- lat2 - lat1 dlon <- lon2 - lon1 a <- sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2 c <- 2 * atan2(sqrt(a), sqrt(1-a)) return(R * c) } dy <- haversine(c(e[3],e[1]), c(e[4],e[1])) # kilometers y-axis dx <- haversine(c(e[3],e[1]), c(e[3],e[2])) # kilometers x-axis e <- extent(e) if (dy > dx) { e[1:2] <- matrix(e * dy/dx)[1:2,] # expand x } else { e[3:4] <- matrix(e * dx/dy)[3:4,] # expand y } dx <- as.numeric( # how large a scale bar to make? as.character( cut(dx,breaks=c(0,5,50,200,Inf),labels=c(1,5,10,50),incl=T) ) ) e <- extent(e) * 1.25 # slight expansion } ras <- disaggregate(crop(ras, e), 10, method='bilinear') ### custom function for filled contours.... `fillcont` <- function (x, maxpixels = 1e+05, ...) { `fc` <- function ( x = seq(0, 1, length.out = nrow(z)), y = seq(0, 1, length.out = ncol(z)), z, xlim = range(x, finite = TRUE), ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE), levels = pretty(zlim, nlevels), nlevels = 20, color.palette = function(n) viridis::viridis(n,begin=0,end=0.9), col = color.palette(length(levels) - 1), ...) { if (missing(z)) { if (!missing(x)) { if (is.list(x)) { z <- x$z y <- x$y x <- x$x } else { z <- x x <- seq.int(0, 1, length.out = nrow(z)) } } else stop("no 'z' matrix specified") } else if (is.list(x)) { y <- x$y x <- x$x } if (any(diff(x) <= 0) || any(diff(y) <= 0)) stop("increasing 'x' and 'y' values expected") plot.new() plot.window(xlim, ylim, "", xaxs = "i", yaxs = "i", asp = NA) .filled.contour(x, y, z, levels, col) } x <- sampleRegular(x, maxpixels, asRaster = T, useGDAL = T) X <- xFromCol(x, 1:ncol(x)) Y <- yFromRow(x, nrow(x):1) Z <- t(matrix(getValues(x), ncol = x@ncols, byrow=T)[nrow(x):1,]) lonlat <- couldBeLonLat(x, warnings = FALSE) asp <- list(...)$asp if (is.null(asp)) { if (lonlat) { ym <- mean(c(x@extent@ymax, x@extent@ymin)) asp <- 1/cos((ym * pi)/180) } else { asp <- 1 } fc(x = X, y = Y, z = Z, asp = asp, ...) } else { fc(x = X, y = Y, z = Z, ...) } } fillcont(ras, nlevels=20, lwd=lwd.cont) plot(shp, col='transparent', border='red', add=T, lwd=2) if (has_points) { points(ppts, pch=16, col='#000000', cex=1) scalebar(d=dx, xy=c(e[1] + abs(diff(e[1:2]))*0.05, e[3] + abs(diff(e[3:4]))*0.08), type='bar', divs=2, below=paste0('\n', dx, ' km'), adj=c(0.5,-1), lab='', lwd=2, cex=0.7, font=2) } } ### some wilderness areas too large to plot (memory overrun on server) too_big <- c( 'Oregon Islands Wilderness', 'Frank Church-River of No Return Wilderness', 'Kootznoowoo Wilderness', 'Boundary Waters Canoe Area Wilderness' # 'Daniel J. Evans Wilderness' # <<--- AKA Olympic NP, should be ok ) # error catch -- dont plot if shape too large to plot (server memory overrun) if (rx() %in% too_big) { plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE) text(5,5,labels='Map not available',adj=c(0.5,0.5)) } else { # error catch -- dont plot if valid name not in shapefile tryCatch(map_points(pick = rx(), pts=d, shp=w, ras=dem), error = function(e) { plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE) text(5,5,labels='Map not available',adj=c(0.5,0.5)) }) } # ### error catch -- dont plot if valid name not in shapefile # `try_map` <- function() { # tryCatch(map_points(pick = rx(), pts=d, shp=w, ras=dem), # error = function(e) { # plot(1:10, 1:10, type='n', axes=FALSE, ann=FALSE) # text(5,5,labels='Map not available',adj=c(0.5,0.5)) # }) # } # try_map() # # map_points(pick = rx(), pts=d, shp=w, ras=dem)
\newpage
Site scores
You can think of site scores (or "air scores") as the "average" position of species along an air quality gradient. Technically, these are the community-weighted mean of species' sensitivities, where sensitivities are defined as the value of N or S deposition (kg ha^-1^ y^-1^) at which each species was most frequently detected. We report both raw and climate-adjusted scores, and their nationwide ranking as a percentile.^[See "Indicator: Air scores" section for more information.].
Nitrogen
N: mean deposition
get_tab(pick = rx(), xvar=c('tenyr','n_airscore.mean','nadj.mean','cmaq_n_3yroll.mean'), col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'), caption= 'Site scores for estimated nitrogen deposition (N), as a *mean* (kg ha^-1^ y^-1^).')
N: percentiles
get_tab(pick = rx(), xvar=c('tenyr','p_n_airscore.mean','p_n_adj.mean','p_n_cmaq.mean'), col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'), caption= 'Site scores for estimated nitrogen deposition (N), as a *percentile* (0-100).')
Sulfur
S: mean deposition
get_tab(pick = rx(), xvar=c('tenyr','s_airscore.mean','sadj.mean','cmaq_s_3yroll.mean'), col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'), caption= 'Site scores for estimated sulfur deposition (S), as a *mean* (kg ha^-1^ y^-1^).')
S: percentiles
get_tab(pick = rx(), xvar=c('tenyr','p_s_airscore.mean','p_s_adj.mean','p_s_cmaq.mean'), col.names=c('Ten-year interval', 'Lichen airscore (raw)', 'Lichen airscore (adj)', 'CMAQ model'), caption= 'Site scores for estimated sulfur deposition (S), as a *percentile* (0-100).')
\newpage
Critical loads exceedances
Exceedances are the departure of raw or climate-adjusted lichen airscores from established lichen-based critical loads^[See "Critical loads we used" section.], expressed as the mean of a ten-year period. Negative values indicate the "safety margin" below the critical load threshold, positive values indicate exceedance above the threshold.
Nitrogen
N exceedance: raw lichen airscores
get_tab(pick = rx(), xvar=c('tenyr', 'n_raw_ex_sr.mean', 'n_raw_ex_sr_sens.mean', 'n_raw_ex_sr_forag.mean', 'n_raw_ex_sr_cyano.mean', 'n_raw_ex_compn.mean'), col.names=c('Ten-year interval', 'Richness', 'Sensitive richness', 'Forage lichen abund', 'Cyanolichen abund', 'Comm comp'), caption= 'CL exceedances for mean nitrogen deposition (kg N ha^-1^ y^-1^), based on raw lichen airscores.')
N exceedance: climate-adjusted lichen airscores
get_tab(pick = rx(), xvar=c('tenyr', 'n_adj_ex_sr.mean', 'n_adj_ex_sr_sens.mean', 'n_adj_ex_sr_forag.mean', 'n_adj_ex_sr_cyano.mean', 'n_adj_ex_compn.mean'), col.names=c('Ten-year interval', 'Richness', 'Sensitive richness', 'Forage lichen abund', 'Cyanolichen abund', 'Comm comp'), caption= 'CL exceedances for mean nitrogen deposition (kg N ha^-1^ y^-1^), based on climate-adjusted lichen airscores.')
Sulfur
S exceedance: raw lichen airscores
get_tab(pick = rx(), xvar=c('tenyr', 's_raw_ex_sr.mean', 's_raw_ex_sr_sens.mean', 's_raw_ex_sr_forag.mean', 's_raw_ex_sr_cyano.mean', 's_raw_ex_compn.mean'), col.names=c('Ten-year interval', 'Richness', 'Sensitive richness', 'Forage lichen abund', 'Cyanolichen abund', 'Comm comp'), caption= 'CL exceedances for mean sulfur deposition (kg S ha^-1^ y^-1^), based on raw lichen airscores.')
S exceedance: climate-adjusted lichen airscores
get_tab(pick = rx(), xvar=c('tenyr', 's_adj_ex_sr.mean', 's_adj_ex_sr_sens.mean', 's_adj_ex_sr_forag.mean','s_adj_ex_sr_cyano.mean', 's_adj_ex_compn.mean'), col.names=c('Ten-year interval', 'Richness', 'Sensitive richness', 'Forage lichen abund', 'Cyanolichen abund', 'Comm comp'), caption= 'CL exceedances for mean sulfur deposition (kg S ha^-1^ y^-1^), based on climate-adjusted lichen airscores.')
\newpage
Lichen species observed
Sampling completeness
`get_richness` <- function (pick, ...) { NROW(s[[which(names(s) == pick)]]) } spp_obs <- get_richness(pick = rx()) spp_exp <- 4.01 + 25.81*log(n) # log-linear eqn
Number of surveyed plots: r n
Number of species observed: r spp_obs
Number of species expected: r round(spp_exp, 1)
\
x <- 1:2500 # n plots sampled y <- 4.01 + 25.81*log(x) # n species expected set_par(1) plot(x, y, type='l', log='x', xaxt='n', xlim=c(1,2500), xlab=bquote(Plots~sampled~(log[10]~scale)), ylab='Species expected') atx <- outer(1:9, 10^(0:3)) labx <- ifelse(log10(atx) %% 1 == 0, atx, NA) axis(1, at=atx, labels=labx, las=1) lines(rep(n,2), c(spp_exp,spp_obs)) points(n, spp_exp, pch=21, cex=1.2, bg='white') points(n, spp_obs, pch=16, cex=1.2) legend('topleft', pch=c(16,21), bty='n', cex=0.7, pt.cex=1.0, legend=c(paste0('Observed = ', spp_obs), paste0('Expected = ', round(spp_exp,1))))
Rare and pollution-sensitive species
`get_natureserve` <- function (pick, ...) { tab <- s[[which(names(s) == pick)]] rich <- NROW(tab) nsens <- round(sum(tab$n_class == 'oligotroph') / sum(tab$n_class != '') * 100, 1) ssens <- round(sum(tab$s_class == 'sensitive') / sum(tab$s_class != '') * 100, 1) atrisk <- sum(1 * grepl('1|2|3|Endanger|Concern|Threaten', tab$consrv_stat)) list(rich=rich, nsens=nsens, ssens=ssens, atrisk=atrisk) } x <- get_natureserve(pick = rx())
Number of species observed in total: r x$rich
Number of species with 'at risk' conservation status: r x$atrisk
Percentage of rated species with 'sensitive (oligotroph)' N tolerance: r x$nsens\%
Percentage of rated species with 'sensitive' S tolerance: r x$ssens\%
Species list
List of all macrolichens in plots at r input$x from r rng[1] to r rng[2]. True diversity is likely greater because plot footprints are just a fraction of total Wilderness area. Functional groups ('Fxl grp') are cyanolichens ('cyano') and forage lichens ('forage')(see descriptions). Sensitivity for each species ranges from low to high tolerance for nitrogen (oligo-, meso-, eutroph) and sulfur (sensitive, intermediate, tolerant). Critical loads ('CL') are deposition levels resulting in >20\% decrease in occurrence. Conservation from NatureServe (https://www.natureserve.org/) assigns each species both a global (G) and state (S) rank on a scale of 1 to 5 (1=Critically imperiled; 2=Imperiled; 3=Vulnerable; 4=Apparently secure; 5=Secure).
`get_tab_species` <- function (pick, ...) { tab <- s[which(names(s) == pick)][[1]] tab <- tab[order(tab$n_critload, tab$s_critload, tab$fxlgp),] knitr::kable(x = tab, format = 'markdown', ### format = 'pandoc', row.names = FALSE, longtable = TRUE, booktabs = TRUE, caption = NULL, ### 'Macrolichen species list.', escape = FALSE, col.names = c('Species','Fxl grp','N sensitivity','N CL', 'S sensitivity','S CL','Conservation') ) } get_tab_species(pick = rx())
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.