From c508b5bf79808bd1e82fc7c50dfa089448c56423 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Mon, 2 Nov 2015 20:30:31 -0600
Subject: [PATCH] fixed run.biocro code to properly sample over years, with
 replacement, and return results in order of sample

---
 models/biocro/R/run.biocro.R                  | 25 ++++++++++++-------
 .../regionalbiocro_met_uncertainty.Rscript    |  7 +++---
 2 files changed, 19 insertions(+), 13 deletions(-)

diff --git a/models/biocro/R/run.biocro.R b/models/biocro/R/run.biocro.R
index f0879079aa4..b6f8f2faf54 100644
--- a/models/biocro/R/run.biocro.R
+++ b/models/biocro/R/run.biocro.R
@@ -23,14 +23,20 @@ run.biocro <- function(lat, lon, met.nc = met.nc,
   ## Meteorology
   
   if(met.uncertainty == TRUE){
-    met <- load.cfmet(met.nc, lat = lat, lon = lon, start.date = "1979-01-01", 
-                      end.date = "2010-12-31")
-    year.sample <- met[,list(year = sample(unique(year), size = 15))] 
-    met <- met[year %in% year.sample$year]
+    start.date <- "1979-01-01"
+    end.date <- "2010-12-31"
+    years <- sample(year(start.date):year(end.date), 
+                                 size = 15, replace = TRUE)
+    
+
   } else {
     met <- load.cfmet(met.nc, lat = lat, lon = lon, start.date = start.date, end.date = end.date)
+    years <- year(start.date):year(end.date)
   }
-  
+  met <- load.cfmet(met.nc, lat = lat, lon = lon, 
+                    start.date = start.date, end.date = end.date)
+  met <- met[year %in% years]
+
   met.hr <- cfmet.downscale.time(cfmet = met, output.dt = 1)
   biocro.met <- cf2biocro(met.hr)
   
@@ -46,7 +52,7 @@ run.biocro <- function(lat, lon, met.nc = met.nc,
   }
   soil.parms <- lapply(config$pft$soilControl, as.numeric)
   
-  years <- year(start.date):year(end.date)
+  
   for(yeari in years){
     yearchar <- as.character(yeari)
     WetDat <- biocro.met[biocro.met$year == yeari, ]
@@ -94,7 +100,7 @@ run.biocro <- function(lat, lon, met.nc = met.nc,
                               photoControl=config$pft$photoParms)
       
     } else if (genus == "Miscanthus"){
-      if(yeari == min(years)){
+      if(yeari == years[1]){
         iRhizome <- config$pft$iPlantControl$iRhizome
       } else {
         iRhizome <- last(tmp.result$Rhizome)
@@ -119,9 +125,9 @@ run.biocro <- function(lat, lon, met.nc = met.nc,
                                            Stem, Leaf, Root, Rhizome, Grain, LAI,
                                            SoilEvaporation, CanopyTrans, 
                                            key = c("year", "doy", "hour")))
-    if(yeari == min(years)){
+    if(yeari == years[1]){
       hourly.results <- result.yeari.hourly
-    } else if (yeari > min(years)){
+    } else if (!yeari == years[1]){
       hourly.results <- rbind(hourly.results, result.yeari.hourly)
     }
   }
@@ -130,6 +136,7 @@ run.biocro <- function(lat, lon, met.nc = met.nc,
   setkeyv(hourly.results, c("year", "doy", "hour"))
 
   hourly.results <- merge(biocro.met.dt, hourly.results) ## right join
+  hourly.results <- hourly.results[year %in% years]
   daily.results <- hourly.results[,list(Stem = max(Stem), Leaf = max(Leaf),
                                         Root = max(Root), Rhizome = max(Rhizome),
                                         SoilEvaporation = sum(SoilEvaporation),
diff --git a/models/biocro/inst/regionalbiocro_met_uncertainty.Rscript b/models/biocro/inst/regionalbiocro_met_uncertainty.Rscript
index db74adc4b57..4f49c8f096b 100755
--- a/models/biocro/inst/regionalbiocro_met_uncertainty.Rscript
+++ b/models/biocro/inst/regionalbiocro_met_uncertainty.Rscript
@@ -14,7 +14,6 @@ require(PEcAn.data.atmosphere)
 require(PEcAn.utils)
 require(lubridate)
 
-options(error = browser)
 
 if(interactive()) {
 
@@ -29,8 +28,8 @@ if(interactive()) {
 
 config <- read.biocro.config(file.path(rundir, "config.xml"))
 
-met.nc  <- nc_open(config$run$met.file)
-soil.nc <- nc_open(config$run$soil.file)
+met.nc  <- nc_open("~/.pecan/dbfiles/met/narr/threehourly/illinois.nc")#config$run$met.file)
+soil.nc <- nc_open("~/.pecan/dbfiles/soil/hwsd.nc")#config$run$soil.file)
 # atmco2.nc <- nc_open(file.path(inputdir, "co2/CO2_Global_HD_v1.nc"))
 
 lat <- ncvar_get(met.nc, "latitude")
@@ -45,7 +44,7 @@ for(i in 1:500){
   set.seed(i)
   lat <- latlon$lat[point]
   lon <- latlon$lon[point]
-  out <- run.biocro(lat, lon, met.nc = met.nc, soil.nc = soil.nc, config = config)
+  out <- run.biocro(lat, lon, met.nc = met.nc, soil.nc = soil.nc, config = config, met.uncertainty = TRUE)
   annual.result <- rbind(annual.result, out$annually)
   save(annual.result, file = file.path(outdir, "annual.result.RData"))