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

require(RODBC)
require(lattice)
sws <- odbcConnect("swsartifex") # Tunnel to artifex port 3306

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

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

pdf("~/sci/sws/svg/prelimsoils2006.pdf", width = 8.5, height = 11)
source("http://anthony.darrouzet-nardi.net/works/adnlattice.R")

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))

dev.off()

## WATER SAMPLES IN EPHEMERAL FOREST STREAMS

pwater <- sqlQuery(sws, "SELECT * FROM prelimwater")

pdf("~/sci/sws/svg/prelimwater.pdf", width = 8.5, height = 11)
source("http://anthony.darrouzet-nardi.net/works/adnlattice.R")

plot(northing ~ easting, data = pwater, type = "n", main = "location", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$location), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "TN mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$TN), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "PN mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$PN), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "TDN mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$TDN), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "DON mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$DON), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "InN mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$InN), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "NO3 mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$NO3), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "NH4 mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$NH4), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "SO4 mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$SO4), cex = 0.7)

plot(northing ~ easting, data = pwater, type = "n", main = "CL mg/L", asp = 1)
text(pwater$easting, pwater$northing, as.character(pwater$CL), cex = 0.7)

dev.off()


## SOILS

t <- read.delim("/Users/anthony/sci/sws/data/gps.txt")
t$datetime <- as.POSIXct(strptime(t$datetime, "%m/%d/%y %H:%M"))

t$col[t$cover == "open"] <- "green"
t$col[t$cover == "trees"] <- "gold2"
t$col[t$cover == "willows"] <- "red"
t$col[t$cover == "stream"] <- "blue"
t$col[t$cover == "road"] <- "black"

pdf("~/sci/sws/svg/allpoints.pdf", width = 15.153, height = 11.673)
plot(t$northing ~ t$easting, type = "n", asp = 1, axes = F, ann = F)
text(t$northing ~ t$easting, as.character(t$number), 
	col = t$col, cex = 0.6)
dev.off()

pdf("~/sci/sws/svg/openpoints.pdf", width = 15.153, height = 11.673)
plot(t$northing ~ t$easting, type = "n", asp = 1, axes = F, ann = F)
text(t$northing[t$cover == "open"] ~ t$easting[t$cover == "open"],
	as.character(t$number[t$cover == "open"]), col = "green", cex = 0.6)
dev.off()

pdf("~/sci/sws/svg/treepoints.pdf", width = 15.153, height = 11.673)
plot(t$northing ~ t$easting, type = "n", asp = 1, axes = F, ann = F)
text(t$northing[t$cover == "trees"] ~ t$easting[t$cover == "trees"],
	as.character(t$number[t$cover == "trees"]), col = "gold2", cex = 0.6)
dev.off()

pdf("~/sci/sws/svg/willowpoints.pdf", width = 15.153, height = 11.673)
plot(t$northing ~ t$easting, type = "n", asp = 1, axes = F, ann = F)
text(t$northing[t$cover == "willows"] ~ t$easting[t$cover == "willows"],
	as.character(t$number[t$cover == "willows"]), col = "red", cex = 0.6)
dev.off()


## RESINS

resins <- read.delim("/Users/anthony/sci/sws/data/resins.txt")
t$GPS_time <- as.POSIXct(strptime(t$time, "%m/%d/%y %H:%M"))

pdf("~/sci/sws/svg/resins.pdf", width = 15.153, height = 11.673)
xyplot(resins$northing ~ resins$easting, asp = 1, col = "green",
	pch = 16, groups = "location")
dev.off()


## MODELING

require(RODBC)
require(lattice)
sws <- odbcConnect("swsartifex") # Tunnel to artifex port 3306

sad <- sqlQuery(sws, paste("SELECT datetime, temp, soil_temp, Rh,",
	"soil_moisture, wind_vector, wind_dir, wind_max, solar_radiation FROM saddlemet2006"))
sod <- sqlQuery(sws, paste("SELECT datetime, temp, soiltemp5cm,",
	"soiltemp10cm, soiltemp25cm, soiltemp50cm, ETIprecip, wind_vector,",
	"wind_dir, Rh, radiation_short_atm, wind_max FROM soddiemet2006"))
soil <- sqlQuery(sws, paste("SELECT samplepoint_id, soilpoints2006.bin,",
	"datetime, elevation,",
	"northing, easting, temp, cover, pH, moisture FROM soilpoints2006", 
	"INNER JOIN soilbins2006 ON soilpoints2006.bin = soilbins2006.bin"))

# adding weather data to soils dataframe
sod$rain36h <- 0
for(x in 217:26134) {
sod$rain36h[x] <- sum(sod$ETIprecip[(x-136*6):x]) }

soil$sadneartime <- as.POSIXct("2000-01-01")
for (x in 1:1939) sad[[1]][as.numeric(soil$datetime[x] - sad$datetime) 
	< 0][1] -> soil$sadneartime[x]
for (n in 2:9) {
soil[10+n] <- 1
names(soil)[10+n] <- paste("sad.",names(sad)[n], sep = "")
for (x in 1:1939) sad[[n]][as.numeric(soil$sadneartime)[x] ==
	as.numeric(sad$datetime)] -> soil[[10+n]][x]   }

soil$sodneartime <- as.POSIXct("2000-01-01")
for (x in 1:1939) sod[[1]][as.numeric(soil$datetime[x] - sod$datetime) 
	< 0][1] -> soil$sodneartime[x]
for (n in 2:12) {
soil[19+n] <- 1
names(soil)[19+n] <- paste("sod.",names(sod)[n], sep = "")
for (x in 1:1939) sod[[n]][as.numeric(soil$sodneartime)[x] ==
	as.numeric(sod$datetime)] -> soil[[19+n]][x]   }

# air temp and all soil temps
xyplot(sod$temp + sod$soiltemp5cm + sod$soiltemp10cm + sod$soiltemp25cm +
	sod$soiltemp50cm ~ sod$datetime, type = "l", col = c("black", "brown",
	"red", "orange", "yellow"), au = T, aspect = "xy", lty = 1)

splom(soil[c(4, 7, 9, 10, 12, 13, 15, 16, 17, 19, 21, 22, 23, 24, 25, 27, 32)])


lm(temp ~ datetime+elevation+cover+pH+moisture+sad.temp+sad.soil_temp+
	sad.Rh+sad.soil_moisture+sad.wind_vector+sad.solar_radiation+sod.temp+
	sod.soiltemp5cm+sod.soiltemp10cm+sod.soiltemp25cm+sod.soiltemp50cm+
	sod.wind_vector+sod.Rh+sod.radiation_short_atm+sod.rain36h, data = soil)