updated ghgvc_script.R to add crop defaults and new plots

ebimodeling · Oct 15, 2014 · 9e51e63 · 9e51e63
1 parent 79b5839
commit 9e51e63
Showing 1 changed file with 57 additions and 2 deletions.
diff --git a/src/ghgvc_script.R b/src/ghgvc_script.R
@@ -11,6 +11,37 @@ outdir <- args[2]
 config.xml <- file.path(rundir, "multisite_config.xml")
 config.list <- xmlToList(xmlParse(config.xml))
 
+#### Hack to put in some defaults for biophysical variables
+biophys.defaults <- structure(
+  list(ecosystem = structure(
+    c(7L, 4L, 5L, 3L, 6L, 1L, 2L),
+    .Label = c("BR_soy", "BR_sugarcane", "miscanthus", "soybean", 
+               "Spring Wheat", "switchgrass", "US_corn"), class = "factor"), 
+  Rnet = c(0.01, 0.06, 0.02, 0.02, 0.05, 0.06, 0.17), 
+  latent = c(0.25, 0.26, 0.13, 0.24, 0.25, 0.1, 0.33)), 
+  .Names = c("ecosystem", "sw_radiative_forcing", "latent"), 
+  class = "data.frame", row.names = c(NA, -7L))
+
+sites <- names(config.list)[-which(names(config.list)%in% "options")]
+for (site in sites){
+  for(ecosystem.idx in which(names(config.list[[site]]) %in% "pft")){
+    ecosystem <- config.list[[site]][[ecosystem.idx]]$name
+    eidx <- biophys.defaults$ecosystem == ecosystem
+    if (any(eidx)){
+#      if(!exists(config.list[[site]][[ecosystem.idx]]$sw_radiative_forcing) | 
+#           config.list[[site]][[ecosystem.idx]]$sw_radiative_forcing %in% c(NA, 0)){
+        config.list[[site]][[ecosystem.idx]]$sw_radiative_forcing <- biophys.defaults$sw_radiative_forcing[eidx]
+ #     }
+ #     if(!exists(config.list[[site]][[ecosystem.idx]]$latent) | 
+ #          config.list[[site]][[ecosystem.idx]]$latent %in% c(NA, 0)){
+        config.list[[site]][[ecosystem.idx]]$latent <- biophys.defaults$latent[eidx]
+ #     }
+
+    }
+  }
+} 
+
+
 x <- ghgvc2(config.list)
 
 writeLines(x, file.path(outdir, "output.json"))
@@ -35,6 +66,26 @@ for(site in names(outlist)){
         outdf <- rbind(outdf, tmpdf)
     }
 }
+
+
+
+## Hack to use crop defaults for demo
+crop.defaults <- structure(list(ecosystem = structure(c(7L, 4L, 5L, 3L, 6L, 1L, 
+                                       2L), .Label = c("BR_soy", "BR_sugarcane", "miscanthus", "soybean", 
+                                                       "Spring Wheat", "switchgrass", "US_corn"), class = "factor"), 
+               sw_radiative_forcing = c(3.84, 34.17, 10.47, 9.72, 24.88, 
+                                        30.12, 90.67), latent = c(136.27, 141.19, 70.59, 132.53, 
+                                                                  137.34, 55.43, 179.2)), .Names = c("name", "swRFV", 
+                                                                                                     "latent"), class = "data.frame", row.names = c(NA, -7L))
+
+for(name in outdf$name){
+  if(name %in% crop.defaults$name){
+    outdf$swRFV[outdf$name == name] <- crop.defaults[crop.defaults$name == name, "swRFV"]
+    outdf$latent[outdf$name == name] <- crop.defaults[crop.defaults$name == name, "latent"]*2
+  }
+
+}
+
 
 ## Cleaning up output for downloading
 
@@ -69,9 +120,13 @@ plotdata$CRV_BIOPHYS <- plotdata$Rnet + plotdata$LE
 plotdata$CRV_NET <- plotdata$CRV_BGC + plotdata$CRV_BIOPHYS
 plotdata[is.na(plotdata)] <- 0
 
+plotdata$CRV_BGC[plotdata$CRV_NET == 0] <- 0
+plotdata$CRV_BIOPHYS[plotdata$CRV_NET == 0] <- 0
+plotdata$CRV_NET[plotdata$CRV_NET == 0] <- NA
+
 nolabels <- theme(axis.title = element_blank(), axis.text.y = element_blank(), 
                   axis.ticks.y = element_blank(), legend.position = "top")
-baseplot <- ggplot(data = plotdata) + theme_minimal() + coord_flip() + nolabels 
+baseplot <- ggplot(data = plotdata) + theme_minimal() + coord_flip() + nolabels + theme(panel.grid.major = element_line("black", size = 0.14)) +  geom_hline(aes(yintercept =0)) 
 
 
 ## T_E is number of years considered
@@ -108,7 +163,7 @@ bgc.plot <- ghgvc.subplot(c("Storage", "Ongoing_Exchange"), data = longdata) +
   scale_fill_manual(values= brewer_pal(pal = "Greens")(6)[c(4,6)], labels = c("Storage", "Ongoing Exchange")) + labs(fill = "") +
   ggtitle("Biogeochemical") + theme(axis.text.y = element_text(size = 12, hjust = 1))
 
-biophys.plot <- ghgvc.subplot(c("Storage", "Ongoing_Exchange"), data = longdata) + 
+biophys.plot <- ghgvc.subplot(c("LE", "Rnet"), data = longdata) + 
   scale_fill_manual(values = brewer_pal(pal = "Blues")(6)[c(4,6)], labels = c(expression("LE", "R"["net"]))) + labs(fill = "") + 
   ggtitle("Biophysical")