#####################################################
### CODE FOR 
### "The greenhouse gas cost of agricultural intensification with 
### groundwater irrigation in a Midwest US row cropping system"
### McGill, Hamilton, Robertson, and Millar
### Global Change Biology
### Accepted 2018


### Research site information:
### Kellogg Biological Station Long Term Ecological Research 
### (KBS LTER) Resource Gradient Experiment (RGE)  
### and Main Cropping System Experiment (MCSE)
### more info @ https://lter.kbs.msu.edu/research/long-term-experiments/
### @KBSLTER


### 14 Sept 2018 ###
### Code by B. McGill
### bonniemcgill (at) gmail (dot) com
### @BonnSci
#####################################################


# contents
# STEP 1 Global Warming Impact (GWI) of gases dissolved in groundwater
# STEP 2 GWI of soil organic C (SOC) change due to irrigation
# STEP 3 GWI of fossil fuels used for irrigation
# STEP 4 GWI of changes in soil N2O emissions due to irrigation
# STEP 5 Calculate totals, graph for 2 irrigation scenarios:
			# dry year = 0.429 m irr
			# avg year = 0.178 m irr
		 # and add to no-till rainfed GWI (all farming practices), graph
# STEP 6 Calculate greenhouse gas intensity (GWI/yield)










### ### ### ###    setwd("your/directory/here")    ### ### ### ### 
setwd("~/Dropbox/Documents/Hamilton Lab/DATA")

# standard error function
se <- function(x){
	sd(x)/sqrt(length(x))
}









################### 1. Dissolved gases GWI


gas<-read.csv("Dissolved gas irrigation all compiled.csv", header=T, comment.char="#")
# head(gas)
# str(gas)

# note: #'d out 3 outliers whose CV was >0.2 and measurements different by an order of magnitude.  

gas$id <- factor(gas$id)
gas$date <- as.Date(gas$date, format="%m/%d/%y", origin="1980-01-01")
gas$rep <- factor(gas$rep)
gas$year <- factor(gas$year)
gas$site <- factor(gas$site, levels=c("DAIRY", "FSC", "MAGGIE", 
"GULL1", "GULL2", "GULL3", "GULL4", "GULL5", 
"SW1", "SW2", "SW3", "SW4", "SW5", "SW6", 
"SW7", "SW8", "SW9", "SW10", "SW11", "SW12",
"SW13", "SW14", "SW15", "SW16", "SW17"))

# starting from umol/L measured in original headspace 
# subtract expected concentration in atm equilibrated water
# expected value determined using Bunsen solubility coefficient, sample temp. Calculation not shown.
# difference is only gas that dissolved into the water in the subsurface
gas$co2.umsub <- gas$co2.umoll.samp-gas$co2.umoll.equil
gas$n2o.umsub <- gas$n2o.umoll.samp-gas$n2o.umoll.equil


# COMPARE TO SCHLESINGER 2000 (optional). ######
# W.H. Schlesinger / Agriculture, Ecosystems and Environment 82 (2000) 121–127. 
# https://doi.org/10.1016/S0167-8809(00)00221-8

# convert dissolved gas conc to atm: ppm= umol/L *RT, 1 atm=1000*1000 ppm
gas$co2.ppm <- gas$co2.umsub * 0.0821 * gas$t.k
gas$n2o.ppm <- gas$n2o.umsub * 0.0821 * gas$t.k
gas$n2oc.ppm <- gas$n2o.ppm*298

# aggregating ppm values THIS IS PPM IN WATER IN SUBSURFACE, NOT EMISSION
ppm.min.region<- aggregate(cbind(co2.ppm, n2o.ppm, n2oc.ppm)~region, data=gas, FUN="min")
colnames(ppm.min.region) <- c("region", "minco2", "minn2o", "minn2oc")
ppm.max.region<- aggregate(cbind(co2.ppm, n2o.ppm, n2oc.ppm)~region, data=gas, FUN="max")
colnames(ppm.max.region) <- c("region", "maxco2", "maxn2o", "maxn2oc")
ppm.se.region<- aggregate(cbind(co2.ppm, n2o.ppm)~region, data=gas, FUN="se")
colnames(ppm.se.region) <- c("region", "seco2", "sen2o")
ppm.mean.region<- aggregate(cbind(co2.ppm, n2o.ppm, n2oc.ppm)~region, data=gas, FUN="mean")
colnames(ppm.mean.region) <- c("region", "meanco2", "meann2o", "meann2oc")

ppm.means <- merge(ppm.mean.region, ppm.se.region)
ppm.means <- merge(ppm.means, ppm.min.region)
ppm.means <- merge(ppm.means, ppm.max.region)
ppm.means  # 


# subtract atmospheric conc. to get the amount of gas that evades the water when equilibrating with the atmosphere.
# need to convert atmospheric ppm to umol/l, use r and t
t <- mean(gas$t.k) # use average temp of measurements
co2.atm.umoll <- 406/(0.0821*t)
n2o.atm.umoll <- 0.328/(0.0821*t)
(10000/(0.0821*t))*(12+32)*(1/1000) # Schlesinger uses 10,000 ppmv CO2  ~ 18 mg/L CO2

# subtract the atm. conc. from measured conc. TO GET EMISSIONS
gas$co2.umoll.em <- gas$co2.umsub - co2.atm.umoll
gas$n2o.umoll.em <- gas$n2o.umsub - n2o.atm.umoll

# convert emissions to mg/l: umol/L * ug/umol * 1mg/1000ug 
gas$co2.mgl <- gas$co2.umoll.em*(12+32)*(1/1000)
gas$n2o.mgl <- gas$n2o.umoll.em*(28+16)*(1/1000)
# convert N2O mg/L to CO2 equivalents mg/L
gas$n2o.ceq <- gas$n2o.mgl*298


# aggregating mg/l values reported
mgl.min.region<- aggregate(cbind(co2.mgl, n2o.mgl, n2o.ceq)~region, data=gas, FUN="min")
colnames(mgl.min.region) <- c("region", "minco2", "minn2o", "minn2oc")
mgl.max.region<- aggregate(cbind(co2.mgl, n2o.mgl, n2o.ceq)~region, data=gas, FUN="max")
colnames(mgl.max.region) <- c("region", "maxco2", "maxn2o", "maxn2oc")
mgl.se.region<- aggregate(cbind(co2.mgl, n2o.mgl, n2o.ceq)~region, data=gas, FUN="se")
colnames(mgl.se.region) <- c("region", "seco2", "sen2o", "sen2oc")
mgl.mean.region<- aggregate(cbind(co2.mgl, n2o.mgl, n2o.ceq)~region, data=gas, FUN="mean")
colnames(mgl.mean.region) <- c("region", "meanco2", "meann2o", "meann2oc")

mgl.means <- merge(mgl.mean.region, mgl.se.region)
mgl.means <- merge(mgl.means, mgl.min.region)
mgl.means <- merge(mgl.means, mgl.max.region)
mgl.means$totceq <- mgl.means$meanco2 + mgl.means$meann2oc
mgl.means$totceq.se <- mgl.means$seco2 + mgl.means$sen2oc

mgl.means  # values reported come from here

### END SCHLESINGER COMPARISON


str(gas)
# only calculating GWI for KBS irrigation well at RGE
gas <- droplevels(data.frame(subset(gas, gas$site=="FSC")))
gas <- gas[,c(2,3,6, 21:23)]

# convert both to Global Warming Impact 
# g CO2e m-2 yr-1 = mg/L* (m irr water / (m2*y)) * 1000 L/m3 * 1g/1000mg * 1m/1000mm
# 2 of the 1000's cancel
# note that 429 mm of irrigation water = 429 L / m2


gas$co2.gwi.dry <- gas$co2.mgl*0.429			# in the paper we refer to 0.429 mm irrigation as the "high" scenario
gas$n2o.gwi.dry <- gas$n2o.ceq*0.429			# here we call it "dry", as in a dry summer means high irrigation


gas$co2.gwi.avg <-gas$co2.mgl*0.178
gas$n2o.gwi.avg <-gas$n2o.ceq*0.178	

gas$co2.gwi.1m <- gas$co2.mgl*1	
gas$n2o.gwi.1m <- gas$n2o.ceq*1
gas$tot.gwi.dry <- gas$co2.gwi.dry + gas$n2o.gwi.dry
gas$tot.gwi.avg <- gas$co2.gwi.avg + gas$n2o.gwi.avg
gas$tot.gwi.1m <- gas$co2.gwi.1m + gas$n2o.gwi.1m

mean.kbs <- aggregate(cbind(tot.gwi.avg, tot.gwi.dry, tot.gwi.1m) ~region, data=gas, FUN="mean")
colnames(mean.kbs) <- c("source", "avg", "dry", "per_m")
mean.kbs

# SE's for Table 2
se.kbs <- aggregate(cbind(tot.gwi.avg, tot.gwi.dry, tot.gwi.1m) ~region, data=gas, FUN="se")
se.kbs

tot.gwi <- mean.kbs
tot.gwi$source <- factor(tot.gwi$source)
levels(tot.gwi$source)[levels(tot.gwi$source)=="KBS"] <- "dissolved gases"
tot.gwi$source <- as.character(tot.gwi$source)


tot.gwi$avg.se <- se.kbs[1,2]
tot.gwi$dry.se <- se.kbs[1,3]
tot.gwi$per_m.se <- se.kbs[1,4]











################### 2. SOC GWI

# soils collected Nov 2016, 25 cm deep, split into 0-10cm, and 10-25 cm
# 3 cores per plot, 2 plots at F1, F5, F8 rainfed, and 2 plots at F1, F5, F8 irrigated
# 3 CHN replicates measurements per core segment
# samples re-run until CV among replicate measurements was <0.1


dat<-read.csv("LTER RG McGill CN dat.csv", header=T, comment.char="#") 
dat<- dat[,c(2:16, 18:20)]

# convert %C concentration to kg C/m2 using bulk density
# each core consists of the upper 0-10 cm layer and the lower 10-25 cm layer

dat$gCkgs <- dat$c_perc *10  # /100 to convert % to gC/gsoil then *1000 to convert to gC/kg soil, altogether *10


# c_perc/100 = g C / g soil
# bulk density = g soil /cm3 soil 
# 1 kg C/ 1000 g
# 10,000 cm2 / m2
# * 10 or 15 cm (depth)

# kg C per m2
for(i in 1:nrow(dat)){
	if(dat$depth[i]=="0-10"){
		dat$kgCm2[i] <- ((dat$c_perc[i]/100)*dat$bulk.density[i]*10*10000)/1000
		dat$kgNm2[i] <- ((dat$n_perc[i]/100)*dat$bulk.density[i]*10*10000)/1000
	}
	if(dat$depth[i]=="10-25"){
		dat$kgNm2[i] <- ((dat$n_perc[i]/100)*dat$bulk.density[i]*15*10000)/1000
		dat$kgNm2[i] <- ((dat$n_perc[i]/100)*dat$bulk.density[i]*15*10000)/1000
	}
}

# average across 3 CHN replicates
kgC <- aggregate(data=dat, cbind(kgCm2, gCkgs) ~plot + depth + core + irr + tmt + block + sand + silt + clay + yield11.100, FUN="mean")
kgC.se <- aggregate(data=dat, cbind(kgCm2, gCkgs)~plot + depth + core + irr + tmt + block + sand + silt + clay + yield11.100, FUN="se")
kgN <- aggregate(data=dat, kgNm2~plot+depth+core+irr+tmt+block+sand+silt+clay+yield16.100, FUN="mean")
core_ag <- merge(kgC, kgN)
core_ag$cn <- core_ag$kgCm2/core_ag$kgNm2

# for manuscript Table S2 SOC means and SE per tmt
kgc <- aggregate(data=dat, cbind(kgCm2, gCkgs)~depth+irr+tmt, FUN="mean")
kgc.se <- aggregate(data=dat, cbind(kgCm2, gCkgs)~depth+irr+tmt, FUN="se")
colnames(kgc.se)[c(4,5)] <- c("kgc.se", "gCkgs.se")
kgc <- merge(kgc, kgc.se)
kgc <- kgc[with(kgc, order(tmt, irr, depth)),]

# F5 0-25 % change due to irrigation (divided by 12 years for per year rate)
# conc (values from kgc above)
(18.7-16.7)/16.7
# pool
((1.67+0.86)-(1.33+0.88))/(1.33+0.88)
# Syswerda et al. % change in conc. due to no-till compared to conventionally tilled
(11.5-10.4)/ 10.4
# Syswerda et al. % change in conc. due to early succession compared to conventionally tilled
(14.3-10.4)/ 10.4
# Syswerda et al. % change in SOC pool due to no-till
(3.6-3.2)/3.2
# due to succession
(4.5-3.2)/3.2

# mean BD
mean(dat$bulk.density)
se(dat$bulk.density)

# loading pH data to get the lime treatment info (can also test with pH, but not)
ph<-read.csv("McGill_pH_20170908.csv", header=T, comment.char="#") 
# str(ph)

ph$block <- factor(ph$block)
ph$Nten <- ph$Nfertcorn/10
# otherwise get scaling error in models because N range is 0-246, while %C is like <1
ph16<- subset(ph, ph$year==2016)  # only interested in 2016 data here
# lime.yn only says "Y" in 2012 not every year, so 
ph12<- subset(ph, ph$year==2012)
ph12 <- ph12[,c(5,16)]
colnames(ph12) <- c("plot", "lime.yn12")
ph16 <- merge(ph16, ph12, all.x=TRUE)

ph16simp <- data.frame(ph16$plot, ph16$pH, ph16$Lime.full, ph16$Nfertcorn, ph16$Nten, ph16$lime.yn12)
colnames(ph16simp) <- c("plot", "pH", "lime", "Nfertcorn", "Nten", "lime.yn")


dat2<- merge(core_ag, ph16simp, all.x=TRUE)  # core_ag has only aggregated the replicate CHN measurements for kgCm2 per individual core (still have 3 cores per plot)

#### SOC regression after model selection, not shown here ####
dat2$lime.yn <- relevel(dat2$lime.yn, ref="N")
summary(mk <- lm(kgCm2~irr*lime.yn + irr*depth + Nfertcorn, dat=dat2))

# PREDICTIONS from model for calculating effect size and making plot CI:

# make new df for making predictions from model
n_corn <- c(0,134,168, 246)
irr2 <- c("N", "Y")
depth2 <- c("0-10", "10-25")
lime2 <- c("Y", "N")
block <- c("2", "3", "6", "7")
newdf <- data.frame(expand.grid( irr=irr2, depth=depth2,Nfertcorn=n_corn,lime.yn=lime2,block=block))

# calculate CI
conf.int <- as.data.frame(cbind(newdf, predict(mk, interval="conf", newdata=newdf)))

# aggregate blocks
xag <- aggregate(cbind(fit, lwr, upr)~irr+Nfertcorn+depth,conf.int, "mean")
xag2 <- aggregate(cbind(fit, lwr, upr)~irr+Nfertcorn+depth+lime.yn,conf.int, "mean")
subset(xag2, xag2$Nfertcorn=="134")

xag3 <- aggregate(cbind(fit, lwr, upr)~irr+Nfertcorn,conf.int, "mean")
subset(xag3, xag3$Nfertcorn=="134")

# checking predictions by hand:
# limed
	# IRR
		# 0-10
a<-		1.22+0.419+0.344-0.334-(0.00056*134)   # 1.57
		# 10-25
b<-		a-0.48-0.205   # 0.889
c<-	a+b   # 2.463
	# RAIN
		# 0-10
d<- 		1.22+0.344-(0.00056*134)  # 1.489
e <-  	d-0.48   # 1.009
f <- d+e		#2.498
lime.diff <- c-f   # -0.035

# not limed
	# IRR
		# 0-10
g<-		1.22+0.418-(0.00056*134)  #1.563
h <- 	g-0.48-0.205  # 0.878
i <- g+h  # 2.441
	# RAIN
		# 0-10
j<-		1.22-(0.00056*134)  #1.145
k <- 	j-0.48  # 0.665
l <- j+k  # 1.810
nolime.diff <- i-l  # 0.631
mean(c(lime.diff, nolime.diff))  # 0.298

# checking again
# USING AVERAGE OF LIMED AND UNLIMED (get same result if do lime, unlimed separately then avereage)
# predicted SOC (kg C / m2) (also calculated by hand in Schlesinger spreadsheet)
itop <- xag$fit[4] # IRR , F5 0-10 
rtop <- xag$fit[3] # RAIN , F5 0-10
ibottom <- xag$fit[12]  # IRR , F5 10-25
rbottom <- xag$fit[11]  # RAIN , F5 10-25
itot <- itop + ibottom
rtot <- rtop + rbottom
idiff <- itot-rtot  # 0.2995 kg C/m2 additional SOC due to 12 yr irrigation
idiff  # (old model was 0.277)

# difference in C accumulation between irrigated 0-10cm and 10-25 cm
itop-ibottom
# difference in C accumulation between rainfed 0-10cm and 10-25 cm
rtop-rbottom

# space for time substitution, where rainfed T5 is time 0, SOC measurements in 2016, 2016-2005 = 12 years of experiment
idiff.yr <- idiff/12  # 12 years
idiff.gyr <- idiff.yr*1000  # convert kg to g
irr.tot <- 1.9108  # total m irr applied 2005-2016 
idiff.gyr.mirr <- idiff.gyr/irr.tot  # g C /(m2 yr m irr.)
idiff.co2 <- idiff.gyr.mirr*(44/12)  # 47.9 g CO2/(m2 yr m irr.)  
soc.dry <- -1*idiff.co2*.429  # additional SOC gained in irrigated plots in a "dry" year in MI w/ 0.429 m irr
soc.avg <- -1*idiff.co2*.178  # additional SOC gained in irrigated plots in an "avg" year in MI w/ 0.178 m irr
soc.1m <- -1*idiff.co2*1  # additional SOC gained in irrigated plots in an "avg" year in MI per m of water (not realistic but needed for converting other things per m)

tot.gwi$source <- as.character(tot.gwi$source)
tot.gwi[2,1:4] <- c("soc", soc.avg, soc.dry, soc.1m)

# propogated SE = (SEa^2 + SEb^2 + SEc^2...)^(1/2)
# itot SE (SE of itop + ibottom)
(2*(0.058^2) + 2*(0.07^2) + (0.055^2) + (0.078^2) + (0.0003^2))^(1/2)  # 0.16
# rtot SE (SE of rtop + rbottom)
(2*(0.058^2) + (0.055^2) + (0.0003^2))^(1/2)  # 0.099
# So the SE of itot - rtot
soc.se <- ((0.16^2)+(0.099^2))^(1/2)  # 0.188

# treat SE like we treat idiff
# per 0.178 m irr
a<- (soc.se/12)*(1000/irr.tot)*(44/12)*0.178   # 5.35
# per 0.429 m irr
b<- (soc.se/12)*(1000/irr.tot)*(44/12)*0.429  # 12.9
# per m irr
c<- (soc.se/12)*(1000/irr.tot)*(44/12)*1
tot.gwi[2, 5:7] <- c(a,b,c)
















################### 3. Fossil fuel GWI



# available from EPA power profiler https://www.epa.gov/energy/power-profiler accessed 23 Mar 2018
# sources of fossil fuels, assuming coal is bituminous
source <- c("oil", "gas", "coal.bit")
# proportion of KBS electricity generated by fossil fuel source
perc.kbs <- c(0.012, 0.146, 0.596)
# lbs of CO2 emitted per million BTUs energy per source from EIA:  https://www.eia.gov/tools/faqs/faq.php?id=73&t=11  Accessed 10/23/17
co2.mbtu <- c(161.3, 117, 205.7)
ff <- as.data.frame(cbind(source, perc.kbs, co2.mbtu))
ff$perc.kbs <- as.numeric(as.character(ff$perc.kbs))
ff$co2.mbtu <- as.numeric(as.character(ff$co2.mbtu))

# irrigation water from West & McBride 2012 is 5.32 GJ/ha for 0.4 m water
# so per m of water that's
gj.ha <- 5.32/0.4
# 1 lb = 0.4536 kg
# 1 million BTU = 1.055 GJ so kg CO2 / GJ:
ff$kg.gj <- ff$co2.mbtu*0.4536/1.055

# dry year = 0.429 m irr
ff$kg.ha.dry <- ff$kg.gj*ff$perc.kbs*gj.ha*0.429
# avg year = 0.178 m irr
ff$kg.ha.avg <- ff$kg.gj*ff$perc.kbs*gj.ha*0.178
# per m 
ff$kg.ha.1m <- ff$kg.gj*ff$perc.kbs*gj.ha*1

# total of 3 ff's
tot.kg.ha.dry <- sum(ff$kg.ha.dry)
tot.kg.ha.avg <- sum(ff$kg.ha.avg)
tot.kg.ha.1m <- sum(ff$kg.ha.1m)

# convert kg/ha to g/m2 ... /10
ff.gwi.dry <- tot.kg.ha.dry/10
ff.gwi.avg <- tot.kg.ha.avg/10
ff.gwi.1m <- tot.kg.ha.1m/10

tot.gwi[3,] <- c("ff", ff.gwi.avg, ff.gwi.dry, ff.gwi.1m, NA, NA, NA)
# no measurement replication so no SEs

















################### 4. Soil N2O GWI

# note flux units are N2O-N g/ha/day

n2ox<-read.csv("mean n2o 2013 with precip, temp2.csv", header=T, comment.char="#") 
# daily mean of 4 n2o measurements per chamber >0.5 mg N2O-N/m2/hr--deemed cut off for outliers mostly from 8/14/13
# merged with daily total rainfall (mm) and daily mean air temp (degrees C)

n2ox<-n2ox[,2:15]
n2ox$date <- as.Date(n2ox$date, "%m/%d/%y")
n2ox$irr.in[is.na(n2ox$irr.in)] <- 0
n2ox$irr.mm <- n2ox$irr.in*25.4
n2ox$ip <- n2ox$precipitation + n2ox$irr.mm
n2ox$temp <- n2ox$air_temp_mean

# flux units are mg m-2 hr-1
# convert to g ha-1 d-1
n2ox$flux.ghad <- n2ox$flux*240


# create variable for if a sample date is on the day of N fertilizer application (13 May 2013) to ten days after, or not within that period (before/after = no)
for(i in 1:nrow(n2ox)){
	if(n2ox$date[i] >= (as.Date("2013-05-13")) & n2ox$date[i] < as.Date("2013-05-23")){
		n2ox$ndate[i] <- "yes"
	}else{
		n2ox$ndate[i] <- "no"
	}
}


# create variable for whether or not a sample is before or after crop harvest
for(i in 1:nrow(n2ox)){
	if(n2ox$date[i] >= as.Date("2013-07-19")){
		n2ox$harvest[i] <- "after"
	}else{
		n2ox$harvest[i] <- "before"
	}
}

n2ox$ndate <- factor(n2ox$ndate)
n2ox$harvest <- factor(n2ox$harvest)

# create variable for whether or not sample collected on a wet day (wet being >0 rainfall or, if an irrigated plot, >0 irrigation + rainfall)
# ip is the sum of precip + irr mm
n2ox$ynwet <- ifelse(n2ox$ip==0, 0, 1)
n2ox$ynwet <- factor(as.character(n2ox$ynwet))


# transform flux rate
n2ox$fluxrt3 <- n2ox$flux.ghad^(1/3)   

n2ox$ndate <- relevel(n2ox$ndate, ref="no")
n2ox$harvest <- relevel(n2ox$harvest, ref="before")
n2ox$irr <- relevel(n2ox$irr, ref="Rainfed")
n2ox$rowname <- as.numeric(row.names(n2ox))


### N2O regression after model selection, not shown here #####
summary(h <- lm(fluxrt3~ynwet+temp+irr+Nfert*harvest+ndate, data=n2ox))

# for summary table in supplementary info
fortable <- aggregate(flux~Nfert+irr+ynwet+harvest+ndate, data=n2ox, FUN="mean")
fortablese <- aggregate(flux~Nfert+irr+ynwet+harvest+ndate, data=n2ox, FUN="se")
colnames(fortablese)[6]<- "se"
fortable <- merge(fortable, fortablese)
fortable <- fortable[with(fortable, order(ndate, harvest, ynwet, irr, Nfert)),]

# make predictions
nf <- c(seq(-10,170,10), 112)  # include F6 for predictions
yn2<- c("0", "1")
temp2 <- seq(4, 30, 2)
irr2 <- c("Irrigated", "Rainfed")
harv <- c("before", "after")
ndat <- c("yes", "no")
newdf2 <- expand.grid(Nfert=nf, ynwet=yn2,  temp=temp2, irr=irr2, harvest=harv, ndate=ndat)
pred2 <- cbind(newdf2, predict(object=h, newdata=newdf2, type="response", interval="confidence"))

# back transform predictions
pred2$fit.bt <- pred2$fit^3
pred2$lwr.bt <- pred2$lwr^3
pred2$upr.bt <- pred2$upr^3

quantile(n2ox$flux.ghad, probs=c(.95))
# 95 % of the data are <24

# for irr vs. rainfed at F6 predictions:
ci.agb <- aggregate(cbind(fit.bt, lwr.bt, upr.bt)~irr+Nfert, data=pred2, mean)
F6 <- subset(ci.agb, ci.agb$Nfert==112)
#F6 <- F6[,1:3]

# diff between irr and rainfed predictions at F6
F6[1,3] - F6[2,3]

# SE (because CI = mean - 1.96*SE, mean + 1.96*SE, so SE= (upper-lower)/3.92
# see http://handbook-5-1.cochrane.org/chapter_7/7_7_7_2_obtaining_standard_errors_from_confidence_intervals_and.htm
i.se<- (F6[1,5]-F6[1,4])/3.92
r.se <-(F6[2,5]-F6[2,4])/3.92
(i.se^2 + r.se^2)^(1/2)  # 1.3

# by what proportion does irrigation change predicted N2O flux?
prop.nflux <- (F6[1,3] - F6[2,3])/F6[2,3]   # 55%

# in 2013 irrigated plots received 0.138 m irr water
# so proportion n flux change per m water:
prop.nflux.m <- prop.nflux / 0.138

# From Gelfand et al. 2016 global change biology doi: 10.1111/gcb.13426
# Table 3 mean measured kg N2O-N ha yr at the KBS LTER MCSE no-till for rainfed
corn <- 0.91
soy <- 0.71
wheat <- 0.92
nflux <- mean(c(corn,soy,wheat))
nflux/prop.nflux.m

scenario <- c("avg", "dry", "per m")
irr <- c(0.178, 0.429, 1)
soiln2o <- as.data.frame(cbind(scenario, irr))
soiln2o$irr <- as.numeric(as.character(soiln2o$irr))

# in an average or dry year, know estimated rainfed emissions (nflux above)
# how much does irrigation increase in mass n2o-N
soiln2o$kgn.ha.y <- nflux*soiln2o$irr*prop.nflux.m
# in mass of N2O
soiln2o$kgn2o.ha.y <- soiln2o$kgn.ha.y * (44/28)
# in co2 e
soiln2o$kgco2e.ha.y <- soiln2o$kgn2o.ha.y *298
# in g co2e per m2 per yr
soiln2o$gwi <- soiln2o$kgco2e.ha.y/10


tot.gwi[4,1:4] <- c("soil n2o", soiln2o[1,6], soiln2o[2,6], soiln2o[3,6])

# SE from Ilya's mean N2O from T2
((0.06^2)+(0.07^2)+(0.1^2))^(1/2)   # 0.13

# recall that CI = mean -2SE, mean + 2SE
# from model fit, lower, and upper CI ci.agb for 
# Irrigated N112 SE:
(0.056-0.048)/2   # 0.004
# Rainfed N112 SE:
(0.036-0.030)/2   # 0.003
# propogated error from calculating % change
((0.004^2)+(0.003^2))^(1/2)   # 0.005

# propagated error of combining Ilya's SE and the % change SE
e <- ((0.005^2)+(0.13^2))^(1/2)

# treat SE as we treat the mean 
# 0.178 m scenario
h<- (e/0.138)*0.178*(44/28)*298/10  # 7.9
# 0.429 m
i<- (e/0.138)*0.429*(44/28)*298/10  # 18.9
# 1 m
j<- (e/0.138)*1*(44/28)*298/10  # 18.9

tot.gwi[4,5:7] <- c(h,i,j)
















################### 5. Total GWI

tot.gwi$avg <- as.numeric(tot.gwi$avg)
tot.gwi$dry <- as.numeric(tot.gwi$dry)
tot.gwi$per_m <- as.numeric(tot.gwi$per_m)

tot.gwi[5,1] <- "tot"

tot.gwi[5,2:4] <- colSums(tot.gwi[2:4], na.rm=TRUE)
str(tot.gwi)

tot.gwi$avg.se <- as.numeric(tot.gwi$avg.se)
tot.gwi$dry.se <- as.numeric(tot.gwi$dry.se)
tot.gwi$per_m.se <- as.numeric(tot.gwi$per_m.se)

# tot SE
tot.gwi$avg.se2 <- tot.gwi$avg.se^2
tot.gwi$dry.se2 <- tot.gwi$dry.se^2
tot.gwi$per_m.se2 <- tot.gwi$per_m.se^2

tot.se <- as.data.frame(colSums(tot.gwi[8:10], na.rm=TRUE))
tot.se$colSums[1]
tot.se[2,1]
tot.gwi[5,5:7] <- c((tot.se[1,1])^(1/2), (tot.se[2,1])^(1/2), (tot.se[3,1])^(1/2))

tot.gwi.snapshot <- as.data.frame(tot.gwi[,1:7])

tot.gwi <- tot.gwi[,c(1:3, 5,6)]
tot.gwia <- tot.gwi


library(reshape2)

# melt for making plot:
tot.melt <- melt(tot.gwi, id="source")
colnames(tot.melt) <- c("source", "scenario", "gwi")
tot.melt.se <- tot.melt[11:20,]
colnames(tot.melt.se)[3] <- "gwi.se"
tot.melt <- tot.melt[1:10,]
tot.melt.se$scenario <- tot.melt$scenario
tot.melt <- merge(tot.melt, tot.melt.se)

tot.melt <- tot.melt[with(tot.melt, order(scenario)),]
tot.melt$order <- c(rep(c(1,2,4,3,5),2))
tot.melt <- tot.melt[with(tot.melt, order(scenario, order)),]
tot.melt$gwi <- as.numeric(tot.melt$gwi)

tot.gwi[is.na(tot.gwi)] <- 0

# to get SE's for plot below to line up appropriately on the stacked bars:
for(i in 1:nrow(tot.melt)){
	if(tot.melt$source[i] %in% "ff"){
		tot.melt$mean.stack[i] <- tot.melt$gwi[i] + tot.melt$gwi[i+1]
	}
	if(tot.melt$source[i] %in% "dissolved gases"){
		tot.melt$mean.stack[i] <- tot.melt$gwi[i] + tot.melt$gwi[i+1] + 
		tot.melt$gwi[i+2]
	}
	if(tot.melt$source[i] %in% "soc"){
		tot.melt$mean.stack[i] <- tot.melt$gwi[i] + tot.melt$gwi[i-1] 
	}
	if(tot.melt$source[i] %in% c("soil n2o", "soc", "tot")){
		tot.melt$mean.stack[i] <- tot.melt$gwi[i]
	}
}

tot.melt$grp <- c(rep(1,4), 2, rep(1,4), 2)
tot.melt$grp <- factor(tot.melt$grp)

###### make figure 5 stacked bar chart of irrigation gwi
# color palette:
# pal5 <- c(
# "#F46D43",  #1
# "#009E73",  #2
# "#74ADD1",  #4
# "burlywood4",#5
# "#E69F00" #3
# ) 

pal5 <- c(
"#56B4E9",  #1  light blue  GHG ev
"#999999",  #2  grey FF
#"white",  #4  IC
"#E69F00",#5  ORANGE SOC
"#CC79A7" #3  PINK N2O
) 


# # # to get subscripts in legend:
# # iddg <- expression(atop('Dissolved~GHG'), paste('Evasion'))
# # iddg <- paste(atop('Dissolved~GHG'), paste('Evasion'))
# iddg <- paste('Dissolved~GHG\nEvasion')
# idf <- paste('Fossil~Fuels') 
# #idic <- paste('Inorganic~C')
# idsoc <- paste("Soil~Organic~C")
# idsoil <- paste('Soil~N[2]*O')
# #cat(idsoil)






library(ggplot2)
library(stringr)  # for str_wrap

ggplot()+
# plot positive stacked bars separate from negative ones:
geom_bar(data=subset(tot.melt, tot.melt$order %in% c(1,2,3)), aes(x=scenario, y=gwi, fill=reorder(source, order)), stat='identity', position='stack', color='black', size=0.3, width=0.7, alpha=0.7) +
geom_bar(data=subset(tot.melt, tot.melt$order %in% c(4)), aes(x=scenario, y=gwi, fill=reorder(source, -order)), stat='identity', position='stack', color='black', size=0.3, width=0.7, alpha=0.7) +
geom_errorbar(data=subset(tot.melt, tot.melt$order %in% c(1,2,3,5)), aes(x=scenario, ymin=mean.stack-gwi.se, ymax=mean.stack+gwi.se,size=grp), width=0.2, position=position_dodge(width=0.5)) +
geom_errorbar(data=subset(tot.melt, tot.melt$order %in% c(4,5)), aes(x=scenario, ymin=mean.stack-gwi.se, ymax=mean.stack+gwi.se, size=grp), width=0.2,  position=position_dodge(width=0.5)) +
#geom_errorbar(data=subset(tot.melt, tot.melt$source=="tot"), aes(x=scenario, ymin=gwi-gwi.se, ymax=gwi+gwi.se),size=0.3, width=0.2, position=position_dodge(width=0.2)) + #position=position_nudge(x=0.2, y=0)) +
geom_point(data=subset(tot.melt, tot.melt$source=="tot"), aes(x=scenario, y=gwi, color=scenario),size=2, position=position_nudge(x=0.125)) +
scale_color_manual(name="", values=c("black", "black"), breaks="0.18 m", labels="Net GWI")+
# above is only so that we can label the point in a legend
#facet_grid(.~m, labeller=labeller(m=m.label))+
scale_size_manual(values=c(0.3, 0.6), guide=FALSE) +
geom_hline(yintercept=0, size=0.3, color='black')+
scale_x_discrete(labels=c("avg"="0.18", "dry"="0.43"), name="Irrigation (m)")+
scale_y_continuous(breaks=seq(-25, 125, 25))+
ylab(expression("Irrigation GWI (g CO"['2']*'e m'^'-2'~'yr'^'-1'*')'))+
#scale_fill_manual(name="GWI \nComponent", values=pal5,  breaks=c("dissolved gases", "ff", "soil n2o", "ic", "soc"), labels=c('Dissolved GHG\nEvasion', function(x) parse(text=c(idf, idsoil, idic, idsoc))))+
scale_fill_manual(name="GWI\nComponent", values=pal5,  breaks=c("dissolved gases", "ff", "soil n2o",  "soc"), labels=c(
		str_wrap(expression('Groundwater Degassing')), 
		expression('Fossil Fuels'),
		expression('Soil N'['2']*'O'),
		#expression('Inorganic C'),
		expression('SOC')))+
theme(legend.text.align=0,
	panel.grid.minor=element_blank(), 
	panel.grid.major=element_blank() ,
	#panel.grid.major.y=element_line(color="grey60", size=0.1),
	axis.text.x=element_text(size=10, color="black"),
	axis.text.y=element_text(size=10, color="black"),
	axis.title.x=element_text(size=12),
	axis.title.y=element_text(size=13, margin=margin(r=-2)),
	panel.background = element_rect(fill = NA) ,
	panel.border=element_rect(color="grey30", fill=NA, size=0.5),
	strip.text=element_text(size=10),
	legend.text=element_text(size=10,margin=margin(t=3, b=3)),
	legend.title=element_text(size=11),
	legend.key.width=unit(1, "line"),   
	legend.key.height=unit(1.2, 'line')  # size of whole box
)	+
guides(fill=guide_legend(override.aes=list(alpha=0.5)))

## ggsave("fig 5.png", path=("/your/folder/here"), dpi=300)
ggsave("IRR GWI totals 20180921.png", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)	
ggsave("IRR GWI totals 20180921.pdf", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)	




###### make Fig. 6 stacked bar chart for rainfed no-till, no-till with avg. irrigation and no-till with dry year irrigation

# rainfed numbers from Gelfand and Robertson (2015) KBS LTER synthesis book chapter
# available online: https://lter.kbs.msu.edu/open-access/citations/pdfs/3465/original/Gelfand-2015-Ecology-Agric-Landscapes.pdf
# see Gelfand table 12.3


soilc <- -122
nfert <- 33
lime <- 4
p<- 0.3
k<- 1.3
seed<- 7
pest<- 15.5
fuel<- 9
n2o<- 39
ch4<- -1


nt2 <- data.frame(cbind(soilc, nfert, lime, p, k, seed, pest, fuel, n2o, ch4))
nt2[2,] <- nt2[1,] # copying first row
nt2[3,] <- nt2[1,] # again
nt2$scenario <- c("0 m", "0.178 m", "0.429 m")
nt2 <- nt2[,c(11, 1:10)]


nt2$soilc.IRR <- c(0, tot.gwi[2,2], tot.gwi[2,3])
nt2$ff.IRR <- c(0, tot.gwi[3,2], tot.gwi[3,3])
nt2$n2o.IRR <- c(0, tot.gwi[4,2], tot.gwi[4,3])

# diss gases
nt2$diss.gas <- c(0, tot.gwi[1,2], tot.gwi[1,3])
# # ic
# nt2$ic <- c(0, tot.gwi[4,2], tot.gwi[4,3])


nt2$other <- nt2$lime + nt2$p + nt2$k + nt2$seed + nt2$pest	
nt2 <- nt2[, c(1:3, 9:16)]
nt2$tot <- c(sum(nt2[1,2:11]), sum(nt2[2,2:11]), sum(nt2[3,2:11]))	

nt2.melt <- melt(nt2, id="scenario")
colnames(nt2.melt) <- c("scenario", "practice", "gwi")

nt2.melt.tot <- subset(nt2.melt, nt2.melt$practice=="tot")
nt2.melt <- droplevels(as.data.frame(subset(nt2.melt, !(nt2.melt$practice=="tot"))))   #


nt2.melt$order <- as.numeric(c(
		rep(9,3), 
		rep(1,3), 
		rep(3,3), 
		rep(7,3), 
		rep(10,3), 
		rep(8,3), 
		rep(4,3), 
		rep(6,3), 
		rep(5,3), 
		rep(2,3)))
nt2.melt <- nt2.melt[with(nt2.melt, order(order)),]	
	

# dummy variable for highlighting new data in stacked bars
for(i in 1:nrow(nt2.melt)){
	if(nt2.melt$scenario[i] %in% "0 m"){
		nt2.melt$data[i] <- "IG"
		}
	else if(nt2.melt$practice[i] %in% c("diss.gas" ,"ff.IRR" , "n2o.IRR" , "soilc.IRR")){
	  nt2.melt$data[i] <- "BM"
	  }
	else{
		nt2.melt$data[i] <- "IG"
	} 
}
	
pal11 <- c(	
"#0072B2",	# ch4
"#56B4E9", 	# diss.gas
"grey45", 	# ff.IRR
"grey45", 	# fuel
#"white",  	# ic
"#CC79A7",	# n2o
"#CC79A7",	# n2o.IRR
"#009E73",	# nfert   
"#F0E442",	# other  
"#E69F00",	# soilc  
"#E69F00"	# soilc.IRR
	)


ggplot()+
# plot positive stacked bars separate from negative ones:
geom_bar(data=subset(nt2.melt, nt2.melt$order <8), aes(x=scenario, y=gwi, fill=reorder(practice, order), size=data, alpha=data), stat='identity', position='stack', color='black', width=0.7) +
geom_bar(data=subset(nt2.melt, nt2.melt$order >7), aes(x=scenario, y=gwi, fill=reorder(practice, -order), size=data, alpha=data), stat='identity', position='stack', color='black', width=0.7) +
geom_point(data=nt2.melt.tot, aes(x=scenario, y=gwi, color=scenario),shape=21, bg="white", size=1.2, stroke=0.7)+
#facet_grid(.~m, labeller=labeller(m=m.label))+
geom_hline(yintercept=0, size=0.2, color='black')+
scale_color_manual(name="", values=c("black", "black", "black"), breaks="0 m", labels="Net GWI")+
# above is only so that we can label the point in a legend
scale_size_manual(values=c(0.7, 0.2), guide=FALSE) + #name="Data sources", values=c(0.7, 0.2), labels=c("BM"="Present study", "IG"="Gelfand & Robertson\n2015"), guide=FALSE) + 
scale_alpha_manual(values=c(0.9, 0.5), guide=FALSE) + #name="Data sources", labels=c("BM"="Present study", "IG"="Gelfand & Robertson\n2015")) +
scale_x_discrete(labels=c("0 m"="0", "0.178 m"="0.18", "0.429 m"="0.43"), name="Irrigation (m)")+
scale_y_continuous(breaks=seq(-150, 250, 50))+
ylab(expression("Cropping system GWI (g CO"['2']*'e m'^'-2'~'yr'^'-1'*')'))+
#scale_fill_manual(name="GWI \nComponent", values=pal5,  breaks=c("dissolved gases", "ff", "soil n2o", "ic", "soc"), labels=c('Dissolved GHG\nEvasion', function(x) parse(text=c(idf, idsoil, idic, idsoc))))+
scale_fill_manual(name="GWI Component", values=pal11,  breaks=c("nfert", "other","fuel", "diss.gas", "n2o", "soilc", "ch4"), 
	labels=c(
		expression('N Fertilizer'),
		expression('Other Inputs'),
		expression('Fossil Fuels'),
		#expression('Dissolved GHG Evasion'), 
		str_wrap(expression('Groundwater Degassing')),  # 11), 
		expression('Soil N'['2']*'O'),
		expression('SOC'), 
		expression('Soil CH'['4'])))+
theme(legend.text.align=0,
	panel.grid.minor=element_blank(), 
	panel.grid.major=element_blank() ,
	#panel.grid.major.y=element_line(color="grey60", size=0.1),
	axis.text.x=element_text(size=10, color="black"),
	axis.text.y=element_text(size=10, color="black"),
	axis.title.x=element_text(size=12),
	axis.title.y=element_text(size=13, margin=margin(r=-1)),
	panel.background = element_rect(fill = NA) ,
	panel.border=element_rect(color="grey30", fill=NA, size=0.5),
	legend.text=element_text(size=9,margin=margin(t=3, b=3)),
	legend.title=element_text(size=11),
	legend.key.width=unit(1, "line"),   
	legend.key.height=unit(1.2, 'line'),  
	#legend.spacing.y=unit(-0.2, 'cm'),
	legend.key=element_rect(fill="white")
)	+
guides(
#size=guide_legend(override.aes=list(fill= 'grey60'), order=2),  # "#E69F00"
fill=guide_legend(override.aes=list( size=0.2), order=1),
#alpha=guide_legend(override.aes=list(alpha=c(0.9, 0.2)), order=2),
color=guide_legend(title=NULL, order=3))

ggsave("IRR GWI T2 rainfed vs IRR 20180921 other inputs lumped.png", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)		 ## ENTER
ggsave("IRR GWI T2 rainfed vs IRR 20180921 other inputs lumped.pdf", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)		

################### 7. GHG Intensity (GHGI)



yield<-read.csv("77-agronomic+yields.csv", header=T, comment.char="#") # resource gradient yields

yield$date <- as.Date(yield$date, format="%Y-%m-%d")

yield <- subset(yield, yield$year >2004)  # at current location in 2005

######### two entries for yield in 2005, harvested in different areas of plots  ###### 


## ** ## ** updated with 2016 data 


# NOTE: fert rate and yield are in kg/ha.
# CAUTION: 2007 wheat was fertilized at corn rates.
# CAUTION: 2014 irrigated plots seeded at 39,000 seeds/ac, rainfed at 30,000 - thus confounding irrigation with seeding rate.  (obscures irrigation effect)

yield$treatment <- factor(yield$treatment)
levels(yield$irrigated)[levels(yield$irrigated)=="f"] <- "Rainfed"
levels(yield$irrigated)[levels(yield$irrigated)=="t"] <- "Irrigated"
yield$yield_Mg <- (yield$yield/1000)

se <- function(x){
	sd(x)/sqrt(length(x))
}

# making Nfert column from fertilizer_...  for brevity
yield$Nfert <- yield$fertilizer_rate_kg_ha
yield$fertilizer_rate_kg_ha <- NULL

# making a "irr_rain" column from irrigated column for consistency with other code
yield$irr_rain <- yield$irrigated

# Nfert levels missing for 2015,2016
for (i in 1:nrow(yield)){
     
    if( yield$year[i]=='2015') {   # soy
        yield$Nfert[i] <- 0 
     
     }
    if( yield$year[i]=='2016' & yield$treatment[i]==1) {  #wheat
        yield$Nfert[i] <- 0
    }
   
   if( yield$year[i]=='2016' & yield$treatment[i]==2) {
        yield$Nfert[i] <- 23
    }    
      
    if( yield$year[i]=='2016' & yield$treatment[i]==3) {
        yield$Nfert[i] <- 45
    }  
	
	if( yield$year[i]=='2016' & yield$treatment[i]==4) {
        yield$Nfert[i] <- 67
    }     
   	
   	if( yield$year[i]=='2016' & yield$treatment[i]==5) {
        yield$Nfert[i] <- 90
    }  
   	
   	if( yield$year[i]=='2016' & yield$treatment[i]==6) {
        yield$Nfert[i] <- 112
    }  

	if( yield$year[i]=='2016' & yield$treatment[i]==7) {
        yield$Nfert[i] <- 135
    }  
   
   	if( yield$year[i]=='2016' & yield$treatment[i]==8) {
        yield$Nfert[i] <- 157
    }     
   
    if( yield$year[i]=='2016' & yield$treatment[i]==9) {
        yield$Nfert[i] <- 180
    }    
   
   if( yield$year[i]=='2017' & yield$treatment[i]==1) {  # corn
        yield$Nfert[i] <- 0
    }
   
   if( yield$year[i]=='2017' & yield$treatment[i]==2) {
        yield$Nfert[i] <- 34
    }    
      
    if( yield$year[i]=='2017' & yield$treatment[i]==3) {
        yield$Nfert[i] <- 67
    }  
	
	if( yield$year[i]=='2017' & yield$treatment[i]==4) {
        yield$Nfert[i] <- 101
    }     
   	
   	if( yield$year[i]=='2017' & yield$treatment[i]==5) {
        yield$Nfert[i] <- 134
    }  
   	
   	if( yield$year[i]=='2017' & yield$treatment[i]==6) {
        yield$Nfert[i] <- 168
    }  

	if( yield$year[i]=='2017' & yield$treatment[i]==7) {
        yield$Nfert[i] <- 202
    }  
   
   	if( yield$year[i]=='2017' & yield$treatment[i]==8) {
        yield$Nfert[i] <- 246
    }     
   
    if( yield$year[i]=='2017' & yield$treatment[i]==9) {
        yield$Nfert[i] <- 291
    }    
    
    
}



for(i in 1:nrow(yield)){
	if(yield$crop[i]=="Glycine max L. (*)"){
		yield$name[i] <- "soybean"
	}
	if(yield$crop[i]=="Triticum aestivum L. (*)"){
		yield$name[i] <- "wheat"
	}
	if(yield$crop[i]=="Zea mays L. (*)"){
		yield$name[i] <- "corn"
	}
}




for (i in 1:nrow(yield)){
		if(yield$year[i]=='2005'){
			yield$irr.m[i] <-0.217
		}
		 if(yield$year[i]=='2006'){
			yield$irr.m[i] <-0.070
		}
		if(yield$year[i]=='2007'){
			yield$irr.m[i] <-0.110
		}
		if(yield$year[i]=='2008'){
			yield$irr.m[i] <-0.145
		}
		 if(yield$year[i]=='2009'){
			 yield$irr.m[i] <-0.119
		}
		if(yield$year[i]=='2010'){
			yield$irr.m[i] <-0.019
		}
		if(yield$year[i]=='2011'){
			yield$irr.m[i] <-0.093
		}
		if(yield$year[i]=='2012'){
			yield$irr.m[i] <-0.429
		}
		if(yield$year[i]=='2013'){
			yield$irr.m[i] <-0.138
		}
		if(yield$year[i]=='2014'){
			yield$irr.m[i] <-0.227
		}
		if(yield$year[i]=='2015'){
			yield$irr.m[i]<-0.148
		}
		if(yield$year[i]=='2016'){
			yield$irr.m[i] <-0.195
		}
		if(yield$year[i]=='2017'){
			yield$irr.m[i] <-0.413
		}
		
	}	
	
	
for (i in 1:nrow(yield)){
		if(yield$irr_rain[i]=='Rainfed'){
			yield$irr.m[i] <-0
		}
	}





yield6 <- subset(yield, yield$treatment=="6")
# calculate gwi. Constants below are from the per m of water irrigation from tot.gwi.snapshot$per_m
# as of 8 Apr 2018 the below constants match tot.gwi.snapsht$per_m
yield6$gasev <- yield6$irr.m*39.13
#yield6$caco3 <- yield6$irr.m*1000*0.1*(44/40)*0.22
yield6$soc <- yield6$irr.m*(-47.895)
yield6$fuel <- yield6$irr.m*80.98
yield6$n2o <- yield6$irr.m*157.98
#yield6$ic <- yield6$irr.m*(-25.412)
yield6$t2.em <- rep(109.1, nrow(yield6))
yield6$t2.sk <- rep(-122, nrow(yield6))

yield6$gwi <- yield6$gasev + yield6$soc + yield6$fuel + yield6$n2o+ yield6$t2.em + yield6$t2.sk

yield6$gwi_kgha <- yield6$gwi*10

sub<- yield6[,c(8:22)]  

sub$net.intens <- sub$gwi_kgha/sub$yield_kg_ha


# # not necessary to calculate emission and sink intensity separately, get same answer. 



mean.yr <- aggregate(cbind(gwi_kgha,yield_kg_ha, net.intens)~name+year+irr_rain, sub, "mean")
se.yr <- aggregate(cbind(gwi_kgha,yield_kg_ha, net.intens)~name+year+irr_rain, sub, "se")
# note gwi_kgha.se should be 0, variation among four treatment reps comes from different yields, shows up in yield and net.intens.
colnames(se.yr)[c(4:6)] <-c("gwi_kgha.se", "yield_kgha.se", "net.intens.se")
yr.sum <- merge(mean.yr, se.yr)








# # 
# #### YIELD BY YEAR
# ggplot(data=yr.sum, aes(x=year, y=yield_Mg, color=name)) +
# #geom_errorbar(data=ymn2, aes(x=year, ymin=Mg-Mg.se, ymax=Mg+Mg.se, color=crop), width=0.3, size=0.6, alpha=0.8, show.legend=FALSE) +
# geom_errorbar(aes(ymin=yield_Mg-yield_Mg.se, ymax=yield_Mg+yield_Mg.se), width=1) +
# geom_line(aes(linetype=irr_rain), size=1, alpha=0.7)  +
# geom_point( size=1.2) +
# ylab(expression(bold("Yield (Mg ha"^'-1'*')')))+
# xlab("Year")+
# scale_x_continuous(breaks=seq(2005, 2017, 2), labels=seq(2005, 2017, 2), limits=c(2004, 2018)) +
# scale_y_continuous(breaks=seq(0,15, 2))+
# scale_linetype_discrete(name="Treatment") +
# scale_color_discrete(name="Crop")+
# theme(panel.grid.major.y=element_blank(),
# panel.grid.minor.y=element_blank(),
# panel.grid.minor.x=element_line(color="grey80", size=0.3),
# panel.grid.major.x=element_line(color="grey80", size=0.3),
# panel.background = element_rect(fill = "white"),
# panel.border=element_rect(color="grey50", size=1, fill=NA),
# axis.title.x=element_text(face="bold", size=11),
# axis.title.y=element_text(size=12),
# axis.text=element_text(size=10)
# )

# ggsave("Yield by year.png", path=("/Users/Bonnie/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)


# # a1 <- a + 
# # geom_line(data=rain, aes(x=year, y=value/40, group=varb), size=2, alpha=0.3, show.legend=FALSE) +
	# # scale_y_continuous(sec.axis=sec_axis(~.*40, name="Growing Season Water Inputs (mm)"))

# #### GHGI BY YEAR
# ggplot(data=yr.sum, aes(x=year, y=net.intens, color=name)) +
# geom_hline(yintercept=0, color="grey40", size=0.5)+
# geom_errorbar(aes(ymin=net.intens-net.intens.se, ymax=net.intens+net.intens.se)) +
# geom_line(aes(linetype=irr_rain), size=1, alpha=0.7)  +
# geom_point( size=1.2) +
# ylab(expression(bold("Yield (Mg ha"^'-1'*')')))+
# xlab("Year")+
# scale_x_continuous(breaks=seq(2005, 2017, 2), labels=seq(2005, 2017, 2), limits=c(2004, 2018)) +
# #scale_y_continuous(breaks=seq(0,15, 2))+
# scale_linetype_discrete(name="Treatment") +
# scale_color_discrete(name="Crop")+
# theme(panel.grid.major.y=element_blank(),
# panel.grid.minor.y=element_blank(),
# panel.grid.minor.x=element_line(color="grey80", size=0.3),
# panel.grid.major.x=element_line(color="grey80", size=0.3),
# panel.background = element_rect(fill = "white"),
# panel.border=element_rect(color="grey50", size=1, fill=NA),
# axis.title.x=element_text(face="bold", size=11),
# axis.title.y=element_text(size=12),
# axis.text=element_text(size=10)
# )

# ggsave("GHGI by year.png", path=("/Users/Bonnie/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)




# emsk.mean.yr <- aggregate(cbind(em.intens, sk.intens)~name+year+irr_rain, sub, "mean")
# emsk.se.yr <- aggregate(cbind(em.intens, sk.intens)~name+year+irr_rain, sub, "se")
# colnames(emsk.se.yr)[c(4,5)] <- c("em.se", "sk.se")
# emsk.yr.sum <- merge(emsk.mean.yr, emsk.se.yr)


# ### emissions vs. sink intensity by year
# ggplot(data=emsk.yr.sum, aes(x=year, color=name))+
# geom_errorbar(aes(ymin=em.intens-em.se, ymax=em.intens+em.se), width=0.3, size=0.6, alpha=0.8, show.legend=FALSE) +
# geom_errorbar(aes(ymin=sk.intens-sk.se, ymax=sk.intens+sk.se), width=0.3, size=0.6, alpha=0.8, show.legend=FALSE) +
# geom_hline(yintercept=0, color="grey40", size=0.5)+
# geom_line(aes(y=sk.intens, linetype=irr_rain), size=1, alpha=0.7)  +
# geom_point( aes(y=sk.intens), size=1.2) +
# geom_line(aes(y=em.intens, linetype=irr_rain), size=1, alpha=0.7)  +
# geom_point( aes(y=em.intens), size=1.2) +
# ylab(expression(bold("GHG Intensity (kg CO"['2']*'e Mg'^'-1'~'Yield)')))+
# xlab("Year")+
# scale_x_continuous(breaks=seq(2005, 2017, 2), labels=seq(2005, 2017, 2), limits=c(2004, 2018)) + 
# scale_y_continuous(breaks=seq(-60,60, 10))+
# scale_linetype_discrete(name="Treatment") +
# scale_color_discrete(name="Crop")+
# theme(panel.grid.major.y=element_blank(),
# panel.grid.minor.y=element_blank(),
# panel.grid.minor.x=element_line(color="grey80", size=0.3),
# panel.grid.major.x=element_line(color="grey80", size=0.3),
# panel.background = element_rect(fill = "white"),
# panel.border=element_rect(color="grey50", size=1, fill=NA),
# axis.title.x=element_text(face="bold", size=11),
# axis.title.y=element_text(size=12),
# axis.text=element_text(size=10)
# )

# ggsave("GHGI em, sink by year.png", path=("/Users/Bonnie/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)


# intensity summary table
sub$emissions <- (sub$gasev + sub$fuel + sub$n2o + sub$t2.em)*10  # kg/ha
sub$sinks <- (sub$soc + sub$t2.sk)*10  # kg/ha
sub$em.intens <- sub$emissions/ sub$yield_kg_ha    # kg CO2e / kg yield /yr
sub$sk.intens <- sub$sinks/ sub$yield_kg_ha
# use net.intens from above
intens.mean <- aggregate(cbind(em.intens, sk.intens, net.intens)~name+irr_rain, sub, "mean")
intens.se <- aggregate(cbind(em.intens, sk.intens, net.intens)~name+irr_rain, sub, "se")
colnames(intens.se)[c(3:5)] <- c("em.se", "sk.se", "net.se")
intens.sum <- merge(intens.mean, intens.se)

overall.mean <- aggregate(cbind(em.intens, sk.intens, net.intens)~irr_rain, sub, "mean")
overall.se <- aggregate(cbind(em.intens, sk.intens, net.intens)~irr_rain, sub, "se")
colnames(overall.se) <- c("irr_rain", "em.se", "sk.se", "net.se")
overall.mean <- merge(overall.mean, overall.se)

overall.mean$name <- c("tot", "tot")
overall.mean <- overall.mean[,c(8, 1:7)]

intens.sum <- rbind(intens.sum, overall.mean)

intens.melt <- melt(intens.sum, id.var=c("name", "irr_rain"))
intens <- intens.melt[1:24,]
i.se <- intens.melt[25:48,]
colnames(intens)[c(3,4)] <- c("type", "mean")
intens$se <-i.se$value

m.label <- c("corn"="Corn\n(5 yr)", "soybean"="Soybean\n(4 yr)", "wheat"="Wheat\n(4 yr)", "tot"="Rotational\nAverage\n(13 yr)")

intens$name <- factor(intens$name, levels=c("corn", "soybean", "wheat", "tot"), ordered=TRUE)

ggplot()+
# plot positive stacked bars separate from negative ones:
geom_bar(data=subset(intens, intens$type=="em.intens"), aes(x=name, y=mean, fill=irr_rain), stat='identity', position=position_dodge(width=0.9), color='black', width=0.7, size=0.3) +
geom_bar(data=subset(intens, intens$type=="sk.intens"), aes(x=name, y=mean, fill=irr_rain), stat='identity', position=position_dodge(width=0.9), color='black', width=0.7, size=0.3) +
geom_point(data=subset(intens, intens$type=="net.intens"), aes(x=name, y=mean, group=irr_rain), size=1, position=position_dodge(width=0.9)) + #, color="black")+
geom_errorbar(data=intens, aes(x=name, ymin=mean-se, ymax=mean+se, group=irr_rain), position=position_dodge(width=0.9), width=0.2, size=0.3) +
#facet_grid(.~name, labeller=labeller(name=m.label))+
geom_hline(yintercept=0, size=0.2, color='black')+
scale_x_discrete(name="Cropping system", labels=m.label)+
scale_y_continuous(breaks=seq(-0.4, 0.4, 0.10))+
ylab(expression("Cropping System GHGI (kg CO"['2']*'e kg'^'-1'~'yield'*')'))+
scale_fill_manual(values=c("white", "deepskyblue1"), name="")+

theme(legend.text.align=0,
	panel.grid.minor=element_blank(), 
	panel.grid.major=element_blank() ,
	#panel.grid.major.y=element_line(color="grey60", size=0.1),
	axis.text.x=element_text(size=10, color="black"), #angle=-30, vjust=1.0, hjust=0.2),
	axis.text.y=element_text(size=10, color="black"),
	axis.title.x=element_text(size=12),
	axis.title.y=element_text(size=13),
	panel.background = element_rect(fill = NA) ,
	panel.border=element_rect(color=NA, fill=NA, size=0.5),
	axis.line=element_line(size=0.2, color="grey30"),
	panel.spacing=unit(-0.7, "lines"),
	legend.text=element_text(size=10),
	legend.title=element_text(size=11),
	legend.key.width=unit(1, "line"),   
	legend.key.height=unit(1.2, 'line'),  
	legend.spacing.y=unit(-0.2, 'cm'),
	legend.key=element_rect(fill="white"),
	legend.justification=c(0,0), legend.position=c(0.75, 0.8)
	#legend.background=element_rect(fill="white", color="grey50", size=0.2)
	#strip.text=element_text(size=10, vjust=-0.2)
)	+
guides(fill=guide_legend(title.position="left", hjust=1))


## Fig. S3
ggsave("GHGI em, sink average.png", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)
ggsave("GHGI em, sink average.pdf", path=("~/Dropbox/Documents/Hamilton Lab/DATA/R/figures"), dpi=300)