## LANDSCAPE BIOGEOCHEMISTRY OF A MOUNTAIN ECOSYSTEM
## SODDIE WATERSHED (SWS)

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

## ELEMENTAL ABUNDANCE AND ISOTOPES IN DRY MEADOW AND WET MEADOW

p15N <- sqlQuery(sws, "SELECT * FROM prelimsoils2006")

print(xyplot(percentn ~ site | fraction, data = p15N, layout = c(4,1)),
	split = c(1,1,2,2), more = T)
print(xyplot(percentc ~ site | fraction, data = p15N, layout = c(4,1)),
	split = c(2,1,2,2), more = T)
print(xyplot(d15n ~ site | fraction, data = p15N, layout = c(4,1),
		subset = p15N$d15n > 0),
	split = c(1,2,2,2), more = T)
print(xyplot(d13c ~ site | fraction, data = p15N, layout = c(4,1)),
	split = c(2,2,2,2))


#interrelations of isotope lab data
panel.lmci = function(x,y,...) {
	panel.xyplot(x,y,...)
	panel.lmline(x,y)
	n <- data.frame(x = seq(min(as.numeric(x), na.rm = T), 
		max(as.numeric(x), na.rm = T), .1))
	yhat <- predict(lm(y ~ as.numeric(x)), n, interval = "confidence")
	llines(n$x, yhat[,2], col = "gray49")
	llines(n$x, yhat[,3], col = "gray49")
	ltext(x = max(x)*.8, y = max(y)*0.8, round(summary(lm(y ~ as.numeric(x)))$r.squared, 2)) }

soilbins <- sqlQuery(sws, "SELECT * FROM soilbins2006")

splom(soilbins[c(2,20:22,3,7,11,5,9,6,10)], panel = panel.lmci, groups = soilbins$cover,
	subset = soilbins$cover != "road" & soilbins$cover != "stream",
	col = c("gray2", "gray3", "blue", "green4", "darkolivegreen2"), 
	pch = 16, alpha = .35, cex = .6)

splom(soilbins[c(2,20,22,3,7,11,5,9,6,10)], panel = panel.lmci, groups = soilbins$cover,
	subset = soilbins$cover != "road" & soilbins$cover != "stream",
	col = c("gray2", "gray3", "blue", "green4", "darkolivegreen2"), 
	pch = 16, alpha = .35, cex = .6)

sbtree <- sqlQuery(sws, "SELECT * FROM soilbins2006 WHERE cover = 'trees'")
sbopen <- sqlQuery(sws, "SELECT * FROM soilbins2006 WHERE cover = 'open'")

splom(sbtree[c(25,23,3,20,22)], panel = panel.lmci,
	pch = 16)

splom(sbopen[c(25,23,3,20,22)], panel = panel.lmci,
	pch = 16)



## BURIED BAGS

bb <- sqlQuery(sws, "SELECT * FROM buriedbag2007 NATURAL JOIN soilbins2006")
bbno <- sqlQuery(sws, "SELECT * FROM buriedbag2007 NATURAL JOIN soilbins2006 
	WHERE net_min <5 AND net_min >-5 AND net_nit<2")

splom(bb[c(2,3,15,16,20,36:40,33,34)], panel = panel.lmci, groups = bb$cover,
	subset = bb$cover != "road" & bb$cover != "stream",
	col = c("gray2", "gray3", "blue", "green4", "darkolivegreen2"), 
	pch = 16, alpha = .35, cex = .6)


## SIMPLE TREE VS. OPEN REGRESSIONS

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

soilbins <- sqlQuery(sws, "SELECT * FROM soilbins2006")

lm1 <- lm(total_del15n ~ slope + late_moisture, data = soilbins, 
	subset = cover == "open")
summary(lm1)$r.squared # r2 = 0.05
xyplot(predict(lm1) ~ model.extract(model.frame(lm1), "response"))

lm2 <- lm(total_del15n ~ slope + late_moisture, data = soilbins, 
	subset = cover == "trees")
summary(lm2)$r.squared # r2 = 0.33
xyplot(predict(lm2) ~ model.extract(model.frame(lm2), "response"))

summary(lm(total_del15n ~ total_percentn, data = soilbins, 
	subset = cover == "open"))$r.squared
summary(lm(total_del15n ~ pH , data = soilbins, 
	subset = cover == "open"))$r.squared
summary(lm(total_del15n ~ slope , data = soilbins, 
	subset = cover == "open"))$r.squared
summary(lm(total_del15n ~ late_moisture, data = soilbins, 
	subset = cover == "open"))$r.squared
summary(lm(total_del15n ~ total_percentn, data = soilbins, 
	subset = cover == "trees"))$r.squared
summary(lm(total_del15n ~ pH, data = soilbins, 
	subset = cover == "trees"))$r.squared
summary(lm(total_del15n ~ slope, data = soilbins, 
	subset = cover == "trees"))$r.squared
summary(lm(total_del15n ~ late_moisture, data = soilbins, 
	subset = cover == "trees"))$r.squared

lm2 <- lm(total_del15n ~ slope + late_moisture, data = soilbins, 
	subset = cover == "trees")
summary(lm2)$r.squared # r2 = 0.33



summary(lm(net_min ~ slope + late_moisture, data = bbno, 
	subset = cover == "open"))$r.squared # r2 = 0.07
summary(lm(net_min ~ slope + late_moisture, data = bbno, 
	subset = cover == "trees"))$r.squared # r2 = 0.01
summary(lm(net_nit ~ slope + late_moisture, data = bbno, 
	subset = cover == "open"))$r.squared # r2 = 0.04
summary(lm(net_nit ~ slope + late_moisture, data = bbno, 
	subset = cover == "trees"))$r.squared # r2 = 0.05





summary(lm(total_del15n ~ total_percentn, data = soilbins, 
	subset = cover == "open"))
summary(lm(total_del15n ~ pH, data = soilbins, 
	subset = cover == "open"))
summary(lm(total_del15n ~ slope, data = soilbins, 
	subset = cover == "open"))
summary(lm(total_del15n ~ late_moisture, data = soilbins, 
	subset = cover == "open"))
summary(lm(total_del15n ~ total_percentn, data = soilbins, 
	subset = cover == "trees"))
summary(lm(total_del15n ~ pH, data = soilbins, 
	subset = cover == "trees"))
summary(lm(total_del15n ~ slope, data = soilbins, 
	subset = cover == "trees"))
summary(lm(total_del15n ~ late_moisture, data = soilbins, 
	subset = cover == "trees"))



## MEAN WIND DIRECTION

soddiewind <- sqlQuery(sws, "SELECT wind_dir FROM soddiemet2006")
soddiewind <- as.numeric(soddiewind[[1]])
density(soddiewind)$x[density(soddiewind)$y == max(density(soddiewind)$y)]

saddlewind <- sqlQuery(sws, "SELECT wind_dir FROM saddlemet2006")
saddlewind <- as.numeric(saddlewind[[1]])
density(saddlewind)$x[density(saddlewind)$y == max(density(saddlewind)$y)]



## RESINS

resins <- sqlQuery(sws, "SELECT * FROM resins2007")

require(grid)
panel.lmci = 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)), 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)))$r.squared, 2)), 0.2,0.9) }

print(xyplot(log10(resins$NO3a+1) ~ log10(resins$NO3b+1), panel = panel.lmci,
	layout = c(1,1)), split = c(1,1,3,2), more = T)
print(xyplot(log10(resins$NH4a+1) ~ log10(resins$NH4b+1), panel = panel.lmci,
	layout = c(1,1)), split = c(2,1,3,2), more = T)
print(xyplot(log10(resins$mean_NH4+1) ~ log10(resins$mean_NO3+1), panel = panel.lmci,
	layout = c(1,1)), split = c(1,2,3,2), more = T)
print(xyplot(log10(resins$NH4a+1) ~ log10(resins$NO3b+1), panel = panel.lmci, 
	layout = c(1,1)), split = c(2,2,3,2), more = T)
print(xyplot(log10(resins$NH4b+1) ~ log10(resins$NO3a+1), panel = panel.lmci, 
	layout = c(1,1)), split = c(3,2,3,2), more = F)



## BURIED BAG TIME SERIES

bb08 <- sqlQuery(sws, "SELECT * FROM bb08")
bb08$location[bb08$bin == '376 to 385'] <- 'wet'
bb08$location[bb08$bin == '856 to 865'] <- 'forest'
bb08$location[bb08$bin == '236 to 244'] <- 'wet'
bb08$location[bb08$bin == '674 to 677'] <- 'swamp'
bb08$location[bb08$bin == '582 to 588'] <- 'slope'
bb08$location[bb08$bin == '745 to 753'] <- 'moist'
bb08$location[bb08$bin == '1920 to 1926'] <- 'willow'
bb08$location[bb08$bin == '1808 to 1811'] <- 'krummholz'
bb08$location[bb08$bin == '1659 to 1668'] <- 'moist'
bb08$location[bb08$bin == '1290 to 1299'] <- 'dry'

print(xyplot(log(NO3i+1) ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(1,1,2,2), more = T)
print(xyplot(NO3i ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(1,2,2,2), more = T)
print(xyplot(log(NH4i+1) ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(2,1,2,2), more = T)
print(xyplot(NH4i ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(2,2,2,2), )

print(xyplot(log(net_min+1) ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(1,1,2,2), more = T)
print(xyplot(net_min ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(1,2,2,2), more = T)
print(xyplot(log(net_nit+1) ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(2,1,2,2), more = T)
print(xyplot(net_nit ~ date_i, g = bin, data = bb08, panel = function(...) {
       panel.xyplot(..., type = "l")
       panel.text(..., labels=bb08$location) }), split = c(2,2,2,2), )


xyplot(log(IN+1)+log(NH4i+1)+log(NO3i+1) ~ date_i | 
	reorder(interaction(bb08$location, bb08$bin), bb08$IN, FUN = max),
	data = bb08, type = "l", scales = list(at = unique(bb08$date_i)), aspect = 1)


## SNOW DATA

snow <- 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 FROM soilbins2006 WHERE cover NOT IN ('stream', 'road') GROUP BY bin ORDER BY bin")


## CATIONS

cations <- sqlQuery(sws, "SELECT * FROM cationsview")

splom(cations[cations$Fe < 14,4:10])



## INTERMITTENT STREAMS 2008


## RESHAPING VEGETATION DATA

v <- read.delim("/Users/anthony/sci/sws/rawdata/2009 Vegetation/veg09p.txt")
	# should contain only bin, species, frequency
reshape(v, timevar = "species", idvar = "bin", direction = "wide") -> v2
write.table(v2, "/Users/anthony/sci/sws/rawdata/2009 Vegetation/veg09p.csv",
	sep = ",")
	# go into the csv and rename the columns to just the species names
	# and replace all NA's with 0's
	# and move the column names one column over so they line up right
	# and delete the first column that R adds

veg <- sqlQuery(sws, "SELECT * FROM vegetationview")
k <- names(veg)[6:length(veg)]
j <- NA
for (x in 6:length(veg)) j[x-5] <- sum(veg[,x])
data.frame(species = k, abundance = j)

## COMPARING PERCENT N MEASUREMENTS BETWEEN BINS AND INDIVIDUAL CORES

nc <- sqlQuery(sws, "SELECT soilbins2006.bin, pooldilutions2008.percentn,
	soilbins2006.total_percentn FROM soilbins2006 INNER JOIN pooldilutions2008
	ON soilbins2006.bin = pooldilutions2008.bin")
plot(nc$percentn ~ nc$total_percentn)
summary(lm(nc$percentn ~ nc$total_percentn))

## CONVERT DEL15N TO ATOM PERCENT

to.atom <- function(x, ar = 0.0036764) {
	(100 * ar * (x/1000 + 1)) /
	(1 + ar + (x/1000 + 1)) }
		#C13 VPDB 0.0111796
		#N15 AIR 0.0036764
		#D VSMOW 0.0001557
		#O18 VSMOW 0.0020051



## HOTSPOTS

bb <- sqlQuery(sws, "SELECT * FROM buriedbag2007 NATURAL JOIN soilbins2006")
resins <- sqlQuery(sws, "SELECT * FROM resins2007")
soilbins <- sqlQuery(sws, "SELECT * FROM soilbins2006")
pd <- sqlQuery(sws, "SELECT * FROM pooldilutions2008view")
cations <- sqlQuery(sws, "SELECT * FROM cationsview")
n2 <- read.delim("/Users/anthony/sci/sws/rawdata/2006 Soil samples/del15n_total_gis_polygons.txt")
n2[[1]] -> k
k-mean(k) -> j

# % N
plot(hotspots(soilbins$total_percentn)) # this should probably be done as the del15n was done too.
plot(hotspots(j))

# bb and resins
plot(hotspots(bb$net_nit, classes = 3), pch = 16, numbers = F)
plot(hotspots(bb$net_min, classes = 3), pch = 16, numbers = F)
plot(hotspots(resins$mean_NO3, classes = 3, tail = "hot"), pch = 16, numbers = F)
plot(hotspots(resins$mean_NH4, classes = 3, tail = "hot"), pch = 16, numbers = F)
summary(hotspots(resins$mean_NO3, classes = 3))
summary(hotspots(resins$mean_NH4, classes = 3))
wealth.gap(hotspots(resins$mean_NO3, classes = 3))
wealth.gap(hotspots(resins$mean_NH4, classes = 3))


#NH4 AND NO4
plot(hotspots(pd$NO3i[as.Date(pd$datetime_i) == "2008-07-09"], classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(pd$NO3i[as.Date(pd$datetime_i) == "2008-07-30"], classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(pd$NH4i[as.Date(pd$datetime_i) == "2008-07-09"], classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(pd$NH4i[as.Date(pd$datetime_i) == "2008-07-30"], classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(bb$soil_NO3_i, classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(bb$soil_NH4_i, classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(cations$NO3, classes = 4), disp = 2, num = F, pch = 16)
plot(hotspots(cations$NH4, classes = 4), disp = 2, num = F, pch = 16) # this one is pretty whack



## PH DATA

aa <- sqlQuery(sws, "SELECT * FROM aa")
#aa$bin <- reorder(aa$bin, aa$pH_5week, max)
#xyplot(pH_24h ~ acid | bin, data = aa)
#xyplot(pH_5week ~ acid | bin, data = aa)

aa$pH_5week_st <- NA
for (v in 1:length(levels(aa$bin))) {
	df <- aa[aa$bin == levels(aa$bin)[v],]
	aa$pH_5week_st[aa$bin == levels(aa$bin)[v]] <- spline(df$pH_5week ~ df$acid, 
		method = "natural", 
		n = 6, 
		xmax = median(aa$acid[aa$HNO3_ml == 5], na.rm = T))$y }

aa$pH_24h_st <- NA
for (v in 1:length(levels(aa$bin))) {
	df <- aa[aa$bin == levels(aa$bin)[v],]
	aa$pH_24h_st[aa$bin == levels(aa$bin)[v]] <- spline(df$pH_24h ~ df$acid, 
		method = "natural", 
		n = 6, 
		xmax = median(aa$acid[aa$HNO3_ml == 5], na.rm = T))$y }

par(mfrow = c(10,11))
for (v in 1:20) {
	df <- aa[aa$bin == levels(aa$bin)[v],]
	plot(df$pH_5week ~ df$acid)
	lines(spline(df$pH_5week ~ df$acid, method = "natural", n = 6, 
		xmax = median(aa$acid[aa$HNO3_ml == 5], na.rm = T))$y) }
par(mfrow = c(1,1))

xyplot(pH_5week_st ~ HNO3_ml | bin, data = aa)
xyplot(pH_24h_st ~ HNO3_ml | bin, data = aa)

aar <- reshape(aa[,c(1,8:10)], direction = "wide", v.names = c("pH_24h_st", "pH_5week_st"),
	timevar = "HNO3_ml", idvar = "bin")

aar$pH_5week.dif_1_0 <- aar$pH_5week_st.0 - aar$pH_5week_st.1
aar$pH_5week.dif_2_0 <- aar$pH_5week_st.0 - aar$pH_5week_st.2
aar$pH_5week.dif_3_0 <- aar$pH_5week_st.0 - aar$pH_5week_st.3
aar$pH_5week.dif_4_0 <- aar$pH_5week_st.0 - aar$pH_5week_st.4
aar$pH_5week.dif_5_0 <- aar$pH_5week_st.0 - aar$pH_5week_st.5

sqlQuery(sws, "DROP TABLE pH_spatial_gis ;")
sqlSave(sws, aar, "pH_spatial_gis")
sqlQuery(sws, "CREATE VIEW ")


aadif <- sqlQuery(sws, "SELECT * FROM aadif")

splom(aadif[,c(5:9,10,11,22,26,25,24,27,30)], panel = panel.lmci)

splom(aadif[,c(10,11,5:7,25,24,22,26)], panel = function(x,y, ...) {
	panel.xyplot(x,y, pch = 16, col = gray(.5), cex = 0.5, alpha = 0.5)
	panel.loess(x,y, col = "red", lwd = 2) } )




## MORE PH DATA ANALYSES

pH246 <- read.delim("/Users/anthony/sci/sws/rawdata/2008 pH/pH_246.txt")
pH246$block <- factor(pH246$block)
pH246 <- pH246[-15,] # eliminates an outlier where too much acid was likely added

N0.lme <- lme(exp(-pH_5w)~acidperg, data = pH246, subset = N == 0, random = ~1|block)
N2.lme <- lme(exp(-pH_5w)~acidperg, data = pH246, subset = N == 2, random = ~1|block)
N4.lme <- lme(exp(-pH_5w)~acidperg, data = pH246, subset = N == 4, random = ~1|block)
N6.lme <- lme(exp(-pH_5w)~acidperg, data = pH246, subset = N == 6, random = ~1|block)

xyplot(pH_5w ~ acidperg | factor(N), data = pH246, 
	panel = function(x,y) {
	pred <- range(pH246$acidperg)[2]/100*(1:100)
		panel.xyplot(x,y)
		llines(as.numeric(-log(predict(N0.lme, data.frame(acidperg = pred), level = 0))) ~ pred, col = "red")
		llines(as.numeric(-log(predict(N2.lme, data.frame(acidperg = pred), level = 0))) ~ pred)
		llines(as.numeric(-log(predict(N4.lme, data.frame(acidperg = pred), level = 0))) ~ pred)
		llines(as.numeric(-log(predict(N6.lme, data.frame(acidperg = pred), level = 0))) ~ pred)
		})
		
xyplot(exp(-pH_5w) ~ acidperg | factor(N), data = pH246, panel = function(x,y) {
	panel.xyplot(x,y)
	panel.lmline(x,y)})

N.lme <- lme(exp(-pH_5w)~acidperg:factor(N) - 1, data = pH246, random = ~1|block)

xyplot(pH_5w ~ acidperg | factor(N), data = pH246)
xyplot(pH_5w ~ acidperg | factor(N), g = block, data = pH246, col = "black", pch = c(1,2,8,16,17), cex = 1.5)

pred <- range(pH246$acidperg)[2]/100*(1:100)
p<- predict(N.lme, data.frame(N = rep(c(0,2,4,6), each = 100), acidperg = rep(pred,4)), level = 0)
p2 <- data.frame(N = rep(c(0,2,4,6), each = 100), acidperg = rep(pred,4), pH = -log(p))
xyplot(pH ~ acidperg | factor(N), data = p2, type = "l")

x0<-intervals(N0.lme)
x2<-intervals(N2.lme)
x4<-intervals(N4.lme)
x6<-intervals(N6.lme)
lbounds <-c(x0$fixed[2,1], x2$fixed[2,1], x4$fixed[2,1], x6$fixed[2,1])
slopes <-c(x0$fixed[2,2], x2$fixed[2,2], x4$fixed[2,2], x6$fixed[2,2])
ubounds <-c(x0$fixed[2,3], x2$fixed[2,3], x4$fixed[2,3], x6$fixed[2,3])
nvals <-c(0,2,4,6)

errs = slopes-lbounds
diff = c(slopes[2]-slopes[1],slopes[3]-slopes[1],slopes[4]-slopes[1])
differr = c(sqrt(errs[1]^2+errs[2]^2), sqrt(errs[1]^2+errs[3]^2), sqrt(errs[1]^2+errs[4]^2))
percentdiff = (1/slopes[1])*diff
percentdifferr = (1/slopes[1])*differr


## VEGETATION SURVEY DATA RESHAPE

v <- sqlQuery(sws, "SELECT bin, species, cover FROM vegetation_soilbins ORDER BY species")
	# should contain only bin, species, frequency
reshape(v, timevar = "species", idvar = "bin", direction = "wide") -> v2
write.table(v2, "/Users/anthony/sci/sws/rawdata/2009 Vegetation/veg_soilins_reshape.csv",
	sep = ",", na = "0")
	# go into the csv and rename the columns to just the species names
	# and move the column names one column over so they line up right
	# and delete the first column that R adds
	# Using create_veg_table.sql, go in and put all the columns from the csv into the new table
	# create the new table in terminal
	# Follow upload procedure for the csv table to load the data into the newly created table
	# LOAD DATA LOCAL INFILE '/home/users/anthony/bk/vr.csv'
	#	INTO TABLE veg_resins
	#	FIELDS TERMINATED BY ',' ENCLOSED BY '"'
	#	IGNORE 1 LINES ;
	# for the last part see ## ADDING DATA TO THE ARTIFEX MYSQL SERVER for more details (in otr)


v <- sqlQuery(sws, "SELECT resin, species, cover FROM vegetation_resins ORDER BY species")
reshape(v, timevar = "species", idvar = "resin", direction = "wide") -> v2
write.table(v2, "/Users/anthony/sci/sws/rawdata/2009 Vegetation/veg_resins_reshape.csv",
	sep = ",", na = "0")


## VEG NMDS MODELING
# the bulk of the creation of ordination variables is in another file, veg.ord.R

require(vegan)
require(labdsv)
source("/Users/anthony/sci/sws/code/R Code/gbm2.R")

vr <- sqlQuery(sws, "SELECT * FROM veg_resins_view")
vs <- sqlQuery(sws, "SELECT * FROM veg_soilbins_view")

# EDA
vrv <- vr[,c(5:129,134:138)] # no location data, plot numbers unident, unveg
abuocc(vrv) # GRAPHS EVERYTHING HERE
spc.pres.vr <- apply(vrv>0,2,sum) # presence absence (excluding unveg, unident)
sort(spc.pres.vr)
plot(sort(spc.pres.vr))
spc.mean.vr <- (apply(vrv,2,sum))/spc.pres.vr
plot(sort(spc.mean.vr),main="Cumulative Species Abundance Distribution",
	xlab="Cumulative Number of Species", ylab="Mean Abundance")
plt.pres.vr<-apply(vrv>0,1,sum) # number of species in each plot
plot(sort(plt.pres.vr))
plt.sum.vr<-apply(vrv,1,sum) #to calculate the total cover on each plot
plot(sort(plt.sum.vr))
plot(plt.pres.vr,plt.sum.vr) #to see the relationship between number of species/plot and total cover
spc.dat.vr<-list(spc.pres.vr,spc.mean.vr,plt.pres.vr,plt.sum.vr)
names(spc.dat.vr) <- c("spc.pres","spc.mean","plt.pres","plt.sum")

vr <- sqlQuery(sws, "SELECT * FROM veg_resins_view")
vrv <- vr[vr$cover == "open",]
vrv.open <- vrv[,c(5:129,134:138)]
spc.abun.open <- apply(vrv.open,2,sum)
vrv <- vr[vr$cover == "tree",]
vrv.tree <- vrv[,c(5:129,134:138)]
spc.abun.tree <- apply(vrv.tree,2,sum)
vrv <- vr[vr$cover == "willow" | vr$cover == "shrub",]
vrv.shrub <- vrv[,c(5:129,134:138)]
spc.abun.shrub <- apply(vrv.shrub,2,sum)
sort((spc.abun.open*0.72 + spc.abun.tree*0.24 + spc.abun.shrub*0.04)/
	(sum(spc.abun.open)*0.72+sum(spc.abun.tree)*0.24+sum(spc.abun.shrub)*0.04)*100,
	dec = T)[1:50] # estimate of total abundance of each species across study area

vs <- sqlQuery(sws, "SELECT * FROM veg_soilbins_view")
vsv <- vs[vs$cover == "open",]
vsv.open <- vsv[,c(5:123,128:132)]
spc.abun.open <- apply(vsv.open,2,sum)
vsv <- vs[vs$cover == "trees",]
vsv.tree <- vsv[,c(5:123,128:132)]
spc.abun.tree <- apply(vsv.tree,2,sum)
vsv <- vs[vs$cover == "willows",]
vsv.shrub <- vsv[,c(5:123,128:132)]
spc.abun.shrub <- apply(vsv.shrub,2,sum)
sort((spc.abun.open*0.72 + spc.abun.tree*0.24 + spc.abun.shrub*0.04)/
	(sum(spc.abun.open)*0.72+sum(spc.abun.tree)*0.24+sum(spc.abun.shrub)*0.04)*100,
	dec = T)[1:50] # estimate of total abundance of each species across study area


#ORDINATION
pca.1<-pca(vrv,cor=TRUE,dim=10)
summary(pca.1,dim=10)
varplot.pca(pca.1)
loadings.pca(pca)
plot(pca.1)

points(veg.ord[],vr$Carex_rupestris>0)
points(vrv.bray.nmds,vr$Caltha_leptosepala>0, col = "blue")
points(vrv.bray.nmds,vr$Trollius_laxus>0, col = "blue")
points(vrv.bray.nmds,vr$Picea_engelmannii>0, col = "green")
points(vrv.bray.nmds,vr$Abies_lasiocarpa>25, col = "green")
points(vrv.bray.nmds,vr$unvegetated>25, col = "purple")