From 26835fc98cd1eb4ad9109443840947021d506cde Mon Sep 17 00:00:00 2001
From: mikewin-climsci <michael.g.windisch@gmail.com>
Date: Mon, 30 Nov 2020 16:42:00 +0100
Subject: [PATCH] Update to calcBphEffect to fill missing data

---
 R/calcBphEffect.R | 48 +++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 46 insertions(+), 2 deletions(-)

diff --git a/R/calcBphEffect.R b/R/calcBphEffect.R
index 244bcbb66..c4abb3f65 100644
--- a/R/calcBphEffect.R
+++ b/R/calcBphEffect.R
@@ -13,12 +13,56 @@
 calcBphEffect <-function(){
 
   x <- readSource("BphEffect", convert="onlycorrect")
+  k <- calcOutput("ClimateClass", aggregate = FALSE)
+  weight <- dimSums(calcOutput("LanduseInitialisation", aggregate=FALSE, cellular=TRUE, land="fao", input_magpie=TRUE, years="y1995", round=6), dim=3)
 
-  bph_mask_weight <- calcOutput("BphMask", aggregate="cluster", round=6)
+
+  #mapping to connect cell names with latitudes
+  map <- toolGetMapping(type="cell", name="CountryToCellMapping.csv")
+  map <- data.table(map)
+
+  #assuming 0 was NA before.
+  x[x==0] <- NA
+
+  #check: How many NA cells
+  length(x[is.na(x)])
+
+  #extract climate class names
+  cclass <- getNames(k)
+
+  #select dominating climate class
+  ccl <- apply(k,1,function(k) {
+    cclass[which(k==max(k))]
+  })
+  #add to mapping
+  map <- cbind(map,ccl)
+
+  #global mean, used in case of NA for all cells for a give lat
+  m_glo <- mean(x,na.rm=TRUE)
+  #m_glo <- mean(b,na.rm=FALSE)
+
+  #Loop over climate classes
+  for (sel in cclass) {
+    print(sel)
+    #get the magpie cells corresponding to cl
+    cells <- map[ccl %in% sel,cell]
+    #find cells with NA
+    cells_NA <- cells[is.na(x[cells,,])]
+    #If all cells are NA, use m_glo, otherwise calc mean based on the non NA cells.
+    if(identical(cells,cells_NA)) {
+      x[cells_NA,,] <- m_glo
+      print("Attention: Global mean used!")
+    } else {
+      m <- mean(x[cells,,],na.rm=TRUE)
+      x[cells_NA,,] <- m
+    }
+  }
+  #check: no more NA cells
+  length(x[is.na(x)])
 
   return(list(
     x=x,
-    weight=bph_mask_weight,
+    weight=weight,
     unit="degC",
     description="Biogeophysical temp change of afforestation in degC",
     isocountries=FALSE))