regional_pecan_workflow.Rmd

---
title: "Miscanthus regional"
author: "David LeBauer"
date: "3/18/2015"
output: html_document
---


```{r, echo=FALSE,warning=FALSE}

library(PEcAn.all)

```
### Load PEcAn settings file.

Open and read in settings file for PEcAn run.


```{r, echo=FALSE, eval = FALSE, warning=FALSE}

settings.file <- c("/home/dlebauer/dev/biocro_regional/vignettes/illinois_mxg_settings.xml")

settings <- read.settings(settings.file)
settings$database$bety$dbname <- "ebi_production"
settings$pfts <- get.trait.data(pfts = settings$pfts, 
                                modeltype = "BIOCRO", 
                                dbfiles = settings$run$dbfiles ,
                                database = settings$database$bety, 
                                forceupdate = "AUTO") 
run.meta.analysis(pfts = settings$pfts, random = TRUE, threshold = 1.2,
                  iterations = settings$meta.analysis$iter, 
                  dbfiles = settings$run$dbfiles, 
                  database = settings$database$bety)

run.write.configs(settings, settings$database$bety$write)

settings$run$host$qsub <- "qsub -l walltime=03:00:00,nodes=1:ppn=10 -N @NAME@ -o @STDOUT@ -e @STDERR@"
start.model.runs(settings, settings$database$bety$write) # Start ecosystem model runs

#settings$sensitivity.analysis$variable <- settings$ensemble$variable <- "TVeg"
get.results(settings = settings)

run.sensitivity.analysis()      # Run sensitivity analysis and variance decomposition on model output

run.ensemble.analysis(plot.timeseries=TRUE)  	      # Run ensemble analysis on model output. 
```

### Analysis for CyberGIS Demo

```{r}
library(data.table)
library(ggplot2)


files <- data.table(run = 9982:10481)
#files <- fread("../run/runs.txt")[, list(filename = file.path(V1, "biocro_output.RData"))]

allout <- list()
for(i in 1:length(files$run)){
  run <- files$run[i]
  filename = file.path("~/pecan_remote/205/out", run, "annual.result.RData")
  load(filename)
  allout[[i]] <- cbind(run = run, annual.result[, list(lat, lon, year, yield = Stem + Leaf)])
  rm(annual.result)
}


all <- rbindlist(allout)
save(all, file = "output/illinois_ensemble.RData")
write.csv(all, file = "output/illinois_ensemble.csv", row.names = FALSE)

range <- all[,list(mean = mean(yield)), by = 'lat,lon,year'][,list(range = max(mean) - min(mean)), by = 'lat,lon']

ranks <- all[, list(mean = mean(yield)), by = run][order(mean), ]
#ranks[5*(1:100) - 2]
sampled.runs <- ranks[5*(1:100)-2]

merge(sampled.runs[, list(run)], all, by = 'run', all.x = TRUE, all.y = FALSE)

```

### Climate

```{r climate-summary}
library(data.table)
library(ggplot2)

# copy on roger.ncsa.illinois.edu:/gpfs_scratch/biocro_output.RData 
load("~/biocro_output.RData")
theme_set(theme_bw())

z <- biocro_result[lat == 40.875 & lon == -87.875 ]
zz <- z[,list(lowyields = ifelse(year %in% c(2005,2007),TRUE,FALSE), doy, p = cumsum(precip), gdd = cumsum((tmax + tmin)/2)), by = 'year']

pdf("il_climate.pdf")

ggplot(data = zz, 
       aes(doy, gdd, group = year, color = lowyields)) + 
  geom_line() + 
  ggtitle("Cumulative GDD")

ggplot(data = zz, 
       aes(doy, p, group = year, color = lowyields)) + 
  geom_line() + 
  ggtitle("Cumulative Precipitation")



zz <- biocro_result[ , list(
  p = max(cumsum(precip)), 
  gdd = max(cumsum((tmax + tmin)/2))), 
  by = 'year,lat,lon']

zz[,`:=`( 
  # that frown `:=`( 
  #    is crazy syntax for 'add new variable'
  #    this is the data.table package. 
  #    Its great for this and many other reasons, 
  #    but dplyr is easier to learn and read
  lowyield = ifelse(year %in% c(2005,2007),TRUE,FALSE))]

ggplot(data = zzz, aes(x = year, y = p)) + 
  geom_line(aes(color = factor(lat))) + 
  ggtitle("total annual precip, by latitude")

ggplot(data = zzz, aes(x = year, y = gdd)) + 
  geom_line(aes(color = factor(lat))) + 
  ggtitle("gdd, by latitude")

dev.off()



```



## Yield Plots

```{r}


library(ggplot2)
library(maps)
illinois <- borders("county", "illinois", colour = "white")


theme_set(theme_minimal())
ggplot(data=stats) + geom_tile(aes(lon,lat, fill = mean)) +   scale_fill_gradientn(colours = rainbow(7)) + illinois

load("annual.result.RData")
y <- annual.result[, list(year = rank(yearindex)), by = 'lat,lon']

mean = mean(Stem + Leaf)* 0.67, se = sd(Stem + Leaf)* 0.67), by = 'lat,lon,']
annual.result[, `:=`(mean = mean(Stem + Leaf)* 0.67, se = sd(Stem + Leaf)* 0.67), by = 'lat,lon,']
ggplot(data=range) + geom_tile(aes(lon,lat, fill = range)) +   scale_fill_gradientn(colours = rainbow(7)) + illinois

mean_timeseries <- all[,list(median = median(yield)), by = 'lat,lon,year']
ggplot(data = mean_timeseries)+ geom_tile(aes(lon,lat, fill = median)) +   scale_fill_gradientn(colours = rainbow(4)) + illinois + facet_wrap(~year)
```