## SAVING INITIAL PARAMETERS FOR VARIOGRAMS

save(am.eye, cn.eye, cnh.eye, cnl.eye, hc.eye, hn.eye, jm.eye, lc.eye, ln.eye, ph.eye, tc.eye, tn.eye, file = "/Users/anthony/sci/sws/code/sws_pools_eyefit.Rdata")

## DATA ACQUISITION

require(RODBC)
require(lattice)
require(grid)
require(Hmisc)
require(geoR)
source("http://anthony.darrouzet-nardi.net/works/adnlattice.R")
sws <- odbcConnect("swsartifex") # Tunnel to artifex port 3306

panel.sws = function(x,y,...) {
	panel.xyplot(x,y,...)
	n <- data.frame(x = seq(min(as.numeric(x), na.rm = T), 
		max(as.numeric(x), na.rm = T), .01))
	yhat <- predict(lm(y ~ as.numeric(x), weights = pools$bin_area), n, interval = "confidence")
	llines(n$x, yhat[,1], col = "black")
	llines(n$x, yhat[,2], col = "gray")
	llines(n$x, yhat[,3], col = "gray")
	#panel.loess(x,y, col = "red", lwd = 1.5)
	grid.text(as.character(round(summary(lm(y ~ as.numeric(x), weights = pools$bin_area))$r.squared, 2)), 0.2,0.9) }


pools <- sqlQuery(sws, "SELECT bin, cover, northing_mean as northing, easting_mean as easting, bin_area, light_pc_core, total_percentc, heavy_pc_core, light_pn_core, total_percentn, heavy_pn_core, light_C_to_N, total_C_to_N, heavy_C_to_N, percent_light, percent_heavysmall, percent_heavylarge, heavy_percentn, heavy_percentc, light_percentn, light_percentc FROM soilbins2006 WHERE cover NOT IN ('stream', 'road') GROUP BY bin ORDER BY bin")


## SEMIVARIOGRAMS (Figure 2)

quartz()

load(url("http://anthony.darrouzet-nardi.net/data/sws_pools_eyefit.Rdata"))

layout(matrix(c(rep(1:9)), 3,3))
lc.v <- variog(as.geodata(pools, 3:4, 6), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(lc.v, lwd = 5, main = "Light carbon")
	lc.vfit <- variofit(lc.v, ini = lc.eye, cov.model = "exponential")
	lines.variomodel(lc.vfit)
ln.v <- variog(as.geodata(pools, 3:4, 9), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(ln.v, lwd = 5, main = "Light nitrogen")
	ln.vfit <- variofit(ln.v, ini = ln.eye, cov.model = "exponential")
	lines.variomodel(ln.vfit)
cnl.v <- variog(as.geodata(pools, 3:4, 12), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(cnl.v, lwd = 5, main = "C:N light")
	cnl.vfit <- variofit(cnl.v, ini = cnl.eye, cov.model = "exponential")
	lines.variomodel(cnl.vfit)
tc.v <- variog(as.geodata(pools, 3:4, 7), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(tc.v, lwd = 5, main = "Total carbon")
	tc.vfit <- variofit(tc.v, ini = tc.eye, cov.model = "exponential")
	lines.variomodel(tc.vfit)
tn.v <- variog(as.geodata(pools, 3:4, 10), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(tn.v, lwd = 5, main = "Total nitrogen")
	tn.vfit <- variofit(tn.v, ini = tn.eye, cov.model = "exponential")
	lines.variomodel(tn.vfit)
cn.v <- variog(as.geodata(pools, 3:4, 13), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(cn.v, lwd = 5, main = "C:N")
	cn.vfit <- variofit(cn.v, ini = cn.eye, cov.model = "exponential")
	lines.variomodel(cn.vfit)
hc.v <- variog(as.geodata(pools, 3:4, 8), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(hc.v, lwd = 5, main = "Heavy carbon")
	hc.vfit <- variofit(hc.v, ini = hc.eye, cov.model = "exponential")
	lines.variomodel(hc.vfit)
hn.v <- variog(as.geodata(pools, 3:4, 11), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(hn.v, lwd = 5, main = "Heavy nitrogen")
	hn.vfit <- variofit(hn.v, , ini = hn.eye, cov.model = "exponential")
	lines.variomodel(hn.vfit)
cnh.v <- variog(as.geodata(pools, 3:4, 14), option = "bin", max.dist = 500, 
		uvec = 7, estimator.type = "modulus")
	plot(cnh.v, lwd = 5, main = "C:N heavy")
	cnh.vfit <- variofit(cnh.v, ini = cnh.eye, cov.model = "exponential")
	lines.variomodel(cnh.vfit)


## BASE STATS (Table 1)

bs.mat <- matrix(nrow = 9, ncol = 8, 
	dimnames = list(c("% Soil C, light fraction", "% Soil C, heavy fraction", "% Soil C, total",
		"% Soil N, light fraction", "% Soil N, heavy fraction", "% Soil N, total",
		"C:N ratio, light","C:N ratio, total","C:N ratio, heavy"),
	c("Mean", "SD", "CV%", "Order", "Skew", "min","max", "n")))

#mean
bs.mat[1,1] <- wtd.mean(pools$light_pc_core, w = pools$bin_area)
bs.mat[2,1] <- wtd.mean(pools$heavy_pc_core, w = pools$bin_area)
bs.mat[3,1] <- wtd.mean(pools$total_percentc, w = pools$bin_area)
bs.mat[4,1] <- wtd.mean(pools$light_pn_core, w = pools$bin_area)
bs.mat[5,1] <- wtd.mean(pools$heavy_pn_core, w = pools$bin_area)
bs.mat[6,1] <- wtd.mean(pools$total_percentn, w = pools$bin_area)
bs.mat[7,1] <- wtd.mean(pools$light_C_to_N, w = pools$bin_area)
bs.mat[8,1] <- wtd.mean(pools$total_C_to_N, w = pools$bin_area)
bs.mat[9,1] <- wtd.mean(pools$heavy_C_to_N, w = pools$bin_area)

#standard deviation
bs.mat[1,2] <- sqrt(wtd.var(pools$light_pc_core, w = pools$bin_area, normwt = T))
bs.mat[2,2] <- sqrt(wtd.var(pools$heavy_pc_core, w = pools$bin_area, normwt = T))
bs.mat[3,2] <- sqrt(wtd.var(pools$total_percentc, w = pools$bin_area, normwt = T))
bs.mat[4,2] <- sqrt(wtd.var(pools$light_pn_core, w = pools$bin_area, normwt = T))
bs.mat[5,2] <- sqrt(wtd.var(pools$heavy_pn_core, w = pools$bin_area, normwt = T))
bs.mat[6,2] <- sqrt(wtd.var(pools$total_percentn, w = pools$bin_area, normwt = T))
bs.mat[7,2] <- sqrt(wtd.var(pools$light_C_to_N, w = pools$bin_area, normwt = T))
bs.mat[8,2] <- sqrt(wtd.var(pools$total_C_to_N, w = pools$bin_area, normwt = T))
bs.mat[9,2] <- sqrt(wtd.var(pools$heavy_C_to_N, w = pools$bin_area, normwt = T))

#coefficient of variation as percent
bs.mat[1,3] <- sqrt(wtd.var(pools$light_pc_core, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$light_pc_core, w = pools$bin_area)*100
bs.mat[2,3] <- sqrt(wtd.var(pools$heavy_pc_core, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$heavy_pc_core, w = pools$bin_area)*100
bs.mat[3,3] <- sqrt(wtd.var(pools$total_percentc, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$total_percentc, w = pools$bin_area)*100
bs.mat[4,3] <- sqrt(wtd.var(pools$light_pn_core, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$light_pn_core, w = pools$bin_area)*100
bs.mat[5,3] <- sqrt(wtd.var(pools$heavy_pn_core, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$heavy_pn_core, w = pools$bin_area)*100
bs.mat[6,3] <- sqrt(wtd.var(pools$total_percentn, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$total_percentn, w = pools$bin_area)*100
bs.mat[7,3] <- sqrt(wtd.var(pools$light_C_to_N, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$light_C_to_N, w = pools$bin_area)*100
bs.mat[8,3] <- sqrt(wtd.var(pools$total_C_to_N, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$total_C_to_N, w = pools$bin_area)*100
bs.mat[9,3] <- sqrt(wtd.var(pools$heavy_C_to_N, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$heavy_C_to_N, w = pools$bin_area)*100

#variation by orders of magnitude
orders <- function(x) {log10(max(x, na.rm = T)/min(x, na.rm = T))}
bs.mat[1,4] <- orders(pools$light_pc_core)
bs.mat[2,4] <- orders(pools$heavy_pc_core)
bs.mat[3,4] <- orders(pools$total_percentc)
bs.mat[4,4] <- orders(pools$light_pn_core)
bs.mat[5,4] <- orders(pools$heavy_pn_core)
bs.mat[6,4] <- orders(pools$total_percentn)
bs.mat[7,4] <- orders(pools$light_C_to_N)
bs.mat[8,4] <- orders(pools$total_C_to_N)
bs.mat[9,4] <- orders(pools$heavy_C_to_N)

#skewness

wtd.skew <- function(x, w = NULL) { 
	m3 <- wtd.mean((x-wtd.mean(x, weights = w))^3, w = w) 
	skew <- m3/(sqrt(wtd.var(x, weights = w, normwt = T)^3))
	skew }
	
bs.mat[1,5] <- wtd.skew(pools$light_pc_core, w = pools$bin_area)
bs.mat[2,5] <- wtd.skew(pools$heavy_pc_core, w = pools$bin_area)
bs.mat[3,5] <- wtd.skew(pools$total_percentc, w = pools$bin_area)
bs.mat[4,5] <- wtd.skew(pools$light_pn_core, w = pools$bin_area)
bs.mat[5,5] <- wtd.skew(pools$heavy_pn_core, w = pools$bin_area)
bs.mat[6,5] <- wtd.skew(pools$total_percentn, w = pools$bin_area)
bs.mat[7,5] <- wtd.skew(pools$light_C_to_N, w = pools$bin_area)
bs.mat[8,5] <- wtd.skew(pools$total_C_to_N, w = pools$bin_area)
bs.mat[9,5] <- wtd.skew(pools$heavy_C_to_N, w = pools$bin_area)

# min/max

bs.mat[1,6] <- min(pools$light_pc_core, na.rm = T)
bs.mat[2,6] <- min(pools$heavy_pc_core, na.rm = T)
bs.mat[3,6] <- min(pools$total_percentc, na.rm = T)
bs.mat[4,6] <- min(pools$light_pn_core, na.rm = T)
bs.mat[5,6] <- min(pools$heavy_pn_core, na.rm = T)
bs.mat[6,6] <- min(pools$total_percentn, na.rm = T)
bs.mat[7,6] <- min(pools$light_C_to_N, na.rm = T)
bs.mat[8,6] <- min(pools$total_C_to_N, na.rm = T)
bs.mat[9,6] <- min(pools$heavy_C_to_N, na.rm = T)

bs.mat[1,7] <- max(pools$light_pc_core, na.rm = T)
bs.mat[2,7] <- max(pools$heavy_pc_core, na.rm = T)
bs.mat[3,7] <- max(pools$total_percentc, na.rm = T)
bs.mat[4,7] <- max(pools$light_pn_core, na.rm = T)
bs.mat[5,7] <- max(pools$heavy_pn_core, na.rm = T)
bs.mat[6,7] <- max(pools$total_percentn, na.rm = T)
bs.mat[7,7] <- max(pools$light_C_to_N, na.rm = T)
bs.mat[8,7] <- max(pools$total_C_to_N, na.rm = T)
bs.mat[9,7] <- max(pools$heavy_C_to_N, na.rm = T)

#n
bs.mat[1,8] <- length(pools$light_pc_core[is.na(pools$light_pc_core) == F])
bs.mat[2,8] <- length(pools$heavy_pc_core[is.na(pools$heavy_pc_core) == F])
bs.mat[3,8] <- length(pools$total_percentc[is.na(pools$total_percentc) == F])
bs.mat[4,8] <- length(pools$light_pn_core[is.na(pools$light_pn_core) == F])
bs.mat[5,8] <- length(pools$heavy_pn_core[is.na(pools$heavy_pn_core) == F])
bs.mat[6,8] <- length(pools$total_percentn[is.na(pools$total_percentn) == F])
bs.mat[7,8] <- length(pools$light_C_to_N[is.na(pools$total_C_to_N) == F])
bs.mat[8,8] <- length(pools$total_C_to_N[is.na(pools$total_C_to_N) == F])
bs.mat[9,8] <- length(pools$heavy_C_to_N[is.na(pools$total_C_to_N) == F])


## SPATIAL STATS (Table 2)

ss.mat <- matrix(nrow = 9, ncol = 7,
	dimnames = list(c("% Soil C, light fraction", "% Soil C, heavy fraction", "% Soil C, total",
		"% Soil N, light fraction", "% Soil N, heavy fraction", "% Soil N, total",
		"C:N ratio, light","C:N ratio, total","C:N ratio, heavy"),
	c("Nugget", "Sill", "Nugget % of sill", "phi", "Kappa", "Range", "H")))


# helps finds what phi is multiplied by to get range
f <- function(x,kappa) cov.spatial(x,cov.pars=c(1,1),kappa=kappa)-.05 

# tausq,nugget
ss.mat[1,1] <- lc.vfit$nugget
ss.mat[2,1] <- hc.vfit$nugget
ss.mat[3,1] <- tc.vfit$nugget
ss.mat[4,1] <- ln.vfit$nugget
ss.mat[5,1] <- hn.vfit$nugget
ss.mat[6,1] <- tn.vfit$nugget
ss.mat[7,1] <- cnl.vfit$nugget
ss.mat[8,1] <- cn.vfit$nugget
ss.mat[9,1] <- cnh.vfit$nugget

# sill
ss.mat[1,2] <- lc.vfit$nugget + lc.vfit$cov.pars[1]
ss.mat[2,2] <- hc.vfit$nugget + hc.vfit$cov.pars[1]
ss.mat[3,2] <- tc.vfit$nugget + tc.vfit$cov.pars[1]
ss.mat[4,2] <- ln.vfit$nugget + ln.vfit$cov.pars[1]
ss.mat[5,2] <- hn.vfit$nugget + hn.vfit$cov.pars[1]
ss.mat[6,2] <- tn.vfit$nugget + tn.vfit$cov.pars[1]
ss.mat[7,2] <- cnl.vfit$nugget + cnl.vfit$cov.pars[1]
ss.mat[8,2] <- cn.vfit$nugget + cn.vfit$cov.pars[1]
ss.mat[9,2] <- cnh.vfit$nugget + cnh.vfit$cov.pars[1]

# nugget % of sill
ss.mat[1,3] <-lc.vfit$nugget / (lc.vfit$nugget + lc.vfit$cov.pars[1])*100
ss.mat[2,3] <-hc.vfit$nugget / (hc.vfit$nugget + hc.vfit$cov.pars[1])*100
ss.mat[3,3] <-tc.vfit$nugget / (tc.vfit$nugget + tc.vfit$cov.pars[1])*100
ss.mat[4,3] <-ln.vfit$nugget / (ln.vfit$nugget + ln.vfit$cov.pars[1])*100
ss.mat[5,3] <-hn.vfit$nugget / (hn.vfit$nugget + hn.vfit$cov.pars[1])*100
ss.mat[6,3] <-tn.vfit$nugget / (tn.vfit$nugget + tn.vfit$cov.pars[1])*100
ss.mat[7,3] <-cnl.vfit$nugget / (cnl.vfit$nugget + cnl.vfit$cov.pars[1])*100
ss.mat[8,3] <-cn.vfit$nugget / (cn.vfit$nugget + cnl.vfit$cov.pars[1])*100
ss.mat[9,3] <-cnh.vfit$nugget / (cnh.vfit$nugget + cnh.vfit$cov.pars[1])*100

#phi
ss.mat[1,4] <- lc.vfit$cov.pars[2]
ss.mat[2,4] <- hc.vfit$cov.pars[2]
ss.mat[3,4] <- tc.vfit$cov.pars[2]
ss.mat[4,4] <- ln.vfit$cov.pars[2]
ss.mat[5,4] <- hn.vfit$cov.pars[2]
ss.mat[6,4] <- tn.vfit$cov.pars[2]
ss.mat[7,4] <- cnl.vfit$cov.pars[2]
ss.mat[8,4] <- cn.vfit$cov.pars[2]
ss.mat[9,4] <- cnh.vfit$cov.pars[2]

# kappa
ss.mat[1,5] <- lc.vfit$kappa
ss.mat[2,5] <- hc.vfit$kappa
ss.mat[3,5] <- tc.vfit$kappa
ss.mat[4,5] <- ln.vfit$kappa
ss.mat[5,5] <- hn.vfit$kappa
ss.mat[6,5] <- tn.vfit$kappa
ss.mat[7,5] <- cnl.vfit$kappa
ss.mat[8,5] <- cn.vfit$kappa
ss.mat[9,5] <- cnh.vfit$kappa

# practical range
ss.mat[1,6] <- lc.vfit$cov.pars[2]*3
ss.mat[2,6] <- hc.vfit$cov.pars[2]*3
ss.mat[3,6] <- tc.vfit$cov.pars[2]*3
ss.mat[4,6] <- ln.vfit$cov.pars[2]*3
ss.mat[5,6] <- hn.vfit$cov.pars[2]*3
ss.mat[6,6] <- tn.vfit$cov.pars[2]*3
ss.mat[7,6] <- cnl.vfit$cov.pars[2]*3
ss.mat[8,6] <- cn.vfit$cov.pars[2]*3
ss.mat[9,6] <- cnh.vfit$cov.pars[2]*3

# H (heterogeneity index)
ss.mat[,7] <- (bs.mat[,3]/100) / (((100-ss.mat[,3])/100)*ss.mat[,6]/500)


signif(bs.mat, 2)
signif(ss.mat, 2)

## SOIL HEAVY LIGHT C:N SCATTERPLOTS AND DENSITY PLOT (Figure 3)

quartz()

print(xyplot(light_percentn ~ light_percentc, data = pools, panel = panel.sws, 
	groups = pools$cover, pch = c(1,2,119), cex = .7, col = "black"),
	split = c(1,1,3,1), more = T)

print(xyplot(heavy_percentn ~ heavy_percentc, data = pools, panel = panel.sws,
	groups = pools$cover, pch = c(1,2,119), cex = .7, col = "black"),
	split = c(2,1,3,1), more = T)

lcn <- pools$light_percentc/pools$light_percentn
hcn <- pools$heavy_percentc/pools$heavy_percentn
print(densityplot(~lcn + hcn, pch = c(3,1), col = c("black", "gray"), 
subset = hcn < 25, lty = 1), split = c(3,1,3,1))


## SPLOM OF TOTAL TOTAL C, N, PH, JULY WATER, AUGUST WATER (Figure 4)

quartz()


splom(pools[c(7,10,11,14,15)], panel = panel.sws, groups = pools$cover,
	pch = c(1,2,119), cex = .7, col = "black")


## SPLOM OF C, N, C:N; light, total, heavy

quartz()


splom(pools[c(5:10,15,18,19)], panel = panel.sws, groups = pools$cover,
	pch = c(1,2,119), cex = .7, col = "black")


## DENSITY PLOTS (Figures 3,4)

quartz()


print(densityplot(pools$light_percentc, pch = 124), split = c(1,3,3,3), more = T)
print(densityplot(pools$heavy_percentc, pch = 124), split = c(2,3,3,3), more = T)
print(densityplot(pools$total_percentc, pch = 124), split = c(3,3,3,3), more = T)
print(densityplot(pools$light_percentn, pch = 124), split = c(1,2,3,3), more = T)
print(densityplot(pools$heavy_percentn, pch = 124), split = c(2,2,3,3), more = T)
print(densityplot(pools$total_percentn, pch = 124), split = c(3,2,3,3), more = T)
print(densityplot(pools$jul_moisture, pch = 124), split = c(1,1,3,3), more = T)
print(densityplot(pools$aug_moisture, pch = 124), split = c(2,1,3,3), more = T)
print(densityplot(pools$pH, pch = 124), split = c(3,1,3,3))


## HETEROGENEITY DEMONSTRATION

layout(matrix(c(rep(1:9)), 3,3))
heterogeneity
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,0), nugget = 0),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
#CV
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = 0),
	col = gray(seq(.8,.2,l = 20)), xlab = "", ylab = "")
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = 0),
	col = gray(.5), xlab = "", ylab = "")
#partial sill
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = 0),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = .1),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = .5),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
#range
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.5), nugget = 0),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.3), nugget = 0),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")
set.seed(500)
image(grf(2500, grid = "reg", cov.model = "exponential", 
	cov.pars = c(1,.03), nugget = 0),
	col = gray(seq(1,0,l = 20)), xlab = "", ylab = "")




## SUPPLEMENTARY STATS FOR PAPER

wtd.mean(pools$light_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"])
wtd.mean(pools$total_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"])
wtd.mean(pools$heavy_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"])

wtd.mean(pools$light_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"])
wtd.mean(pools$total_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"])
wtd.mean(pools$heavy_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"])

sqrt(wtd.var(pools$light_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"]))
sqrt(wtd.var(pools$total_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"]))
sqrt(wtd.var(pools$heavy_C_to_N[pools$cover == "open"], w = pools$bin_area[pools$cover == "open"]))

sqrt(wtd.var(pools$light_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"]))
sqrt(wtd.var(pools$total_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"]))
sqrt(wtd.var(pools$heavy_C_to_N[pools$cover == "trees"], w = pools$bin_area[pools$cover == "trees"]))

cndif <- lm(formula = total_C_to_N ~ cover, data = soilbins, 
	subset = cover == "open" | cover == "trees", weights = bin_area)
coef(cndif)
confint(cndif)

cdif <- lm(formula = total_percentc ~ cover, data = soilbins, 
	subset = cover == "open" | cover == "trees", weights = bin_area)
coef(cdif)
confint(cdif)

ndif <- lm(formula = total_percentn ~ cover, data = soilbins, 
	subset = cover == "open" | cover == "trees", weights = bin_area)
coef(ndif)
confint(ndif)


sqrt(wtd.var(pools$light_percentc, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$light_percentc, w = pools$bin_area)*100
sqrt(wtd.var(pools$heavy_percentc, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$heavy_percentc, w = pools$bin_area)*100
sqrt(wtd.var(pools$light_percentn, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$light_percentn, w = pools$bin_area)*100
sqrt(wtd.var(pools$heavy_percentn, w = pools$bin_area, normwt = T))/
	wtd.mean(pools$heavy_percentn, w = pools$bin_area)*100


wtd.mean(pools$percent_heavysmall[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"])
sqrt(wtd.var(pools$percent_heavysmall[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"], normwt = T))
	
wtd.mean(pools$heavy_pc_core[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"])
sqrt(wtd.var(pools$heavy_pc_core[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"], normwt = T))

wtd.mean(pools$heavy_pn_core[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"])
sqrt(wtd.var(pools$heavy_pn_core[pools$cover == "trees"], 
	w = pools$bin_area[pools$cover == "trees"], normwt = T))


wtd.mean(pools$percent_heavysmall[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"])
sqrt(wtd.var(pools$percent_heavysmall[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"], normwt = T))
	
wtd.mean(pools$heavy_pc_core[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"])
sqrt(wtd.var(pools$heavy_pc_core[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"], normwt = T))

wtd.mean(pools$heavy_pn_core[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"])
sqrt(wtd.var(pools$heavy_pn_core[pools$cover == "open"], 
	w = pools$bin_area[pools$cover == "open"], normwt = T))