From b7551c6a0801f66b66a907dbd5fdba8200d900f1 Mon Sep 17 00:00:00 2001
From: Lisanne <lisanne.nauta@wur.nl>
Date: Tue, 16 Jun 2020 16:58:02 +0200
Subject: [PATCH] pathogenflow model can return emissions per sanitation type
 per landuse area per iso id from mapping tool input

---
 Model input/KLA/Input_file_new_KLA_div.csv |   6 ++
 Model scripts/Human_emissions.r            |  53 +---------
 Model scripts/Pathogenflows.R              | 111 +++++++++++++++++++--
 3 files changed, 111 insertions(+), 59 deletions(-)
 create mode 100644 Model input/KLA/Input_file_new_KLA_div.csv

diff --git a/Model input/KLA/Input_file_new_KLA_div.csv b/Model input/KLA/Input_file_new_KLA_div.csv
new file mode 100644
index 0000000..be3ec5b
--- /dev/null
+++ b/Model input/KLA/Input_file_new_KLA_div.csv	
@@ -0,0 +1,6 @@
+iso;subarea;hdi;population;fraction_urban_pop;fraction_pop_under5;sheddingRate;shedding_duration;incidence_urban_under5;incidence_urban_5plus;incidence_rural_under5;incidence_rural_5plus;flushSewer_urb;flushSeptic_urb;flushPit_urb;flushOpen_urb;flushUnknown_urb;pitSlab_urb;pitNoSlab_urb;compostingTwinSlab_urb;compostingTwinNoSlab_urb;compostingToilet_urb;bucketLatrine_urb;containerBased_urb;hangingToilet_urb;openDefecation_urb;other_urb;isShared_urb;sewerLeak_urb;emptied_urb;isWatertight_urb;hasLeach_urb;coverBury_urb;emptiedTreatment_urb;emptyFrequency_urb;pitAdditive_urb;flushElsewhere_urb;pitVIP_urb;pitTraditional_urb;otherLatrine_urb;otherImproved_urb;otherUnimproved_urb;dontKnow_urb;pitLined_urb;pitUnlined_urb;flushSewer_rur;flushSeptic_rur;flushPit_rur;flushOpen_rur;flushUnknown_rur;pitSlab_rur;pitNoSlab_rur;compostingTwinSlab_rur;compostingTwinNoSlab_rur;compostingToilet_rur;bucketLatrine_rur;containerBased_rur;hangingToilet_rur;openDefecation_rur;other_rur;isShared_rur;sewerLeak_rur;emptied_rur;isWatertight_rur;hasLeach_rur;coverBury_rur;emptiedTreatment_rur;emptyFrequency_rur;pitAdditive_rur;flushElsewhere_rur;pitVIP_rur;pitTraditional_rur;otherLatrine_rur;otherImproved_rur;otherUnimproved_rur;dontKnow_rur;pitLined_rur;pitUnlined_rur;FractionPrimarytreatment;FractionSecondarytreatment;FractionTertiarytreatment;FractionQuaternarytreatment;FractionPonds;FractionNontreatment
+1;Central;0.493;63206;1;0.177278887;1.00E+10;7;0.24;0.08;0.24;0.08;0.263426061;0.202395075;0.057237316;0.012221734;0.007753772;0.41097414;0.02447696;0;0;0.006835904;0;0;0;0.011578119;0.003100918;0.363873646;0;0.469944575;0;0.733457483;0;0.436780403;3;None;0;0;0;0;0;0;0;0.23077419;0.211512815;0.263426061;0.202395075;0.057237316;0.012221734;0.007753772;0.41097414;0.02447696;0;0;0.006835904;0;0;0;0.011578119;0.003100918;0.363873646;0;0.469944575;0;0.733457483;0;0.436780403;3;None;0;0;0;0;0;0;0;0.23077419;0.211512815;0.035;0.035;0;0;0.45;0.48
+2;Kawempe;0.493;333024;1;0.177278887;1.00E+10;7;0.24;0.08;0.24;0.08;0.036860622;0.142274476;0.057601288;0.001679583;0.00524847;0.736130931;0.014416174;0;0;0.00033069;0;0;0.0001061;0.00459866;0.000745624;0.439653607;0;0.218648771;0;0.705524168;0;0.211244894;3;None;0;0;0;0;0;0;0;0.321520783;0.429357014;0.036860622;0.142274476;0.057601288;0.001679583;0.00524847;0.736130931;0.014416174;0;0;0.00033069;0;0;0.0001061;0.00459866;0.000745624;0.439653607;0;0.218648771;0;0.705524168;0;0.211244894;3;None;0;0;0;0;0;0;0;0.321520783;0.429357014;0.035;0.035;0;0;0.45;0.48
+3;Makindye;0.493;385309;1;0.177278887;1.00E+10;7;0.24;0.08;0.24;0.08;0.028570368;0.259280693;0.082540034;0.003821692;0.017477095;0.574689697;0.024289707;0;0;0.001400786;0;0;3.45E-05;0.006082428;0.001834099;0.428658959;0;0.248507909;0;0.740031033;0;0.236410709;3;None;0;0;0;0;0;0;0;0.2588795;0.341499759;0.028570368;0.259280693;0.082540034;0.003821692;0.017477095;0.574689697;0.024289707;0;0;0.001400786;0;0;3.45E-05;0.006082428;0.001834099;0.428658959;0;0.248507909;0;0.740031033;0;0.236410709;3;None;0;0;0;0;0;0;0;0.2588795;0.341499759;0.035;0.035;0;0;0.45;0.48
+4;Nakawa;0.493;317023;1;0.177278887;1.00E+10;7;0.24;0.08;0.24;0.08;0.098606697;0.282278627;0.065945258;0.002724668;0.016024352;0.49795729;0.027461236;0;0;0.000559863;0;0;0.000678506;0.006932757;0.000831023;0.389882378;0;0.196882531;0;0.788037462;0;0.190629622;3;None;0;0;0;0;0;0;0;0.215936282;0.310042107;0.098606697;0.282278627;0.065945258;0.002724668;0.016024352;0.49795729;0.027461236;0;0;0.000559863;0;0;0.000678506;0.006932757;0.000831023;0.389882378;0;0.196882531;0;0.788037462;0;0.190629622;3;None;0;0;0;0;0;0;0;0.215936282;0.310042107;0.035;0.035;0;0;0.45;0.48
+5;Rubaga;0.493;383216;1;0.177278887;1.00E+10;7;0.24;0.08;0.24;0.08;0.004654201;0.164638786;0.089541423;0.001799026;0.010791608;0.688829209;0.031009589;0;0;0.000566465;0;0;2.01E-05;0.003735057;0.004392258;0.441673466;0;0.256076888;0;0.98927109;0;0.243857358;3;None;0;0;0;0;0;0;0;0.312687133;0.407713566;0.004654201;0.164638786;0.089541423;0.001799026;0.010791608;0.688829209;0.031009589;0;0;0.000566465;0;0;2.01E-05;0.003735057;0.004392258;0.441673466;0;0.256076888;0;0.98927109;0;0.243857358;3;None;0;0;0;0;0;0;0;0.312687133;0.407713566;0.035;0.035;0;0;0.45;0.48
diff --git a/Model scripts/Human_emissions.r b/Model scripts/Human_emissions.r
index 4ff602d..e0e4712 100644
--- a/Model scripts/Human_emissions.r	
+++ b/Model scripts/Human_emissions.r	
@@ -77,7 +77,6 @@ human.emissions.model.run <- function(scenario,pathogen,isoraster){
   # calc corrected population
   pop_corrected <- correct.population(popurban_grid,popurban_grid,isoraster,HumanData)
   
-  popurban_grid <- pop_corrected$urban
   poprural_grid <- pop_corrected$rural
   
   # CALCULATION OF EMISSIONS PER SANITION TYPE AND SUBAREA
@@ -245,6 +244,7 @@ calc.emissions <- function(scenario,pathogen,HumanData){
       
       int1<-integrate(f2,0,HumanData$storage_flush_urb[b[j]])
       HumanData$survival_pit_flush_urb[b[j]]<-int1$value/HumanData$storage_flush_urb[b[j]]
+      
       int1<-integrate(f2,0,HumanData$storage_flush_rur[b[j]])
       HumanData$survival_pit_flush_rur[b[j]]<-int1$value/HumanData$storage_flush_rur[b[j]]
       
@@ -507,60 +507,13 @@ calc.wttp.emissions.subarea <- function(emissions,WWTP_inputs,HumanData){
   return(emissions)
 }
 
-calc.water.emissions.grid <- function(emissions,isoraster){
-  pathogen_urban_water_pp<-data.frame(iso=emissions$iso,value=emissions$pathogen_urb_waterforgrid_pp)
-  pathogen_urban_water_pp_raster<-subs(isoraster, pathogen_urban_water_pp , subsWithNA=T)
-  
-  pathogen_rural_water_pp<-data.frame(iso=emissions$iso,value=emissions$pathogen_rur_waterforgrid_pp)
-  pathogen_rural_water_pp_raster<-subs(isoraster, pathogen_rural_water_pp , subsWithNA=T)
-  
-  pathogen_urban_water_grid<-pathogen_urban_water_pp_raster*popurban_grid
-  pathogen_rural_water_grid<-pathogen_rural_water_pp_raster*poprural_grid
-  
-  temp<-data.frame(bla=NA,value=0)
-  pathogen_urban_water_grid<-subs(pathogen_urban_water_grid,temp,subsWithNA=F)
-  pathogen_rural_water_grid<-subs(pathogen_rural_water_grid,temp,subsWithNA=F)
-  
-  temp<-data.frame(bla=NaN,value=0)
-  pathogen_urban_water_grid<-subs(pathogen_urban_water_grid,temp,subsWithNA=F)
-  pathogen_rural_water_grid<-subs(pathogen_rural_water_grid,temp,subsWithNA=F)
-  
-  pathogen_water_grid<-pathogen_urban_water_grid+pathogen_rural_water_grid
-  return(pathogen_water_grid)
-}
-
-calc.land.emissions.grid <- function(emissions,isoraster){
-  #for emissions to land
-  pathogen_urban_land_pp<-data.frame(iso=emissions$iso,value=emissions$pathogen_urb_landforgrid_pp)
-  pathogen_urban_land_pp_raster<-subs(isoraster, pathogen_urban_land_pp , subsWithNA=T)
-  
-  pathogen_rural_land_pp<-data.frame(iso=emissions$iso,value=emissions$pathogen_rur_landforgrid_pp)
-  pathogen_rural_land_pp_raster<-subs(isoraster, pathogen_rural_land_pp , subsWithNA=T)
-  
-  pathogen_urban_land_grid<-pathogen_urban_land_pp_raster*popurban_grid
-  pathogen_rural_land_grid<-pathogen_rural_land_pp_raster*poprural_grid
-  
-  temp<-data.frame(bla=NA,value=0)
-  pathogen_urban_land_grid<-subs(pathogen_urban_land_grid,temp,subsWithNA=F)
-  pathogen_rural_land_grid<-subs(pathogen_rural_land_grid,temp,subsWithNA=F)
-  
-  temp<-data.frame(bla=NaN,value=0)
-  pathogen_urban_land_grid<-subs(pathogen_urban_land_grid,temp,subsWithNA=F)
-  pathogen_rural_land_grid<-subs(pathogen_rural_land_grid,temp,subsWithNA=F)
-  
-  pathogen_land_grid<-pathogen_urban_land_grid+pathogen_rural_land_grid
-  return(pathogen_land_grid)
-  
-}
-
 calc.emissions.grid <- function(isoraster,iso,emissions_forgrid_pp,pop_grid){
   pathogen_pp<-data.frame(iso=iso,value=emissions_forgrid_pp)
   pathogen_pp_raster<-subs(isoraster, pathogen_pp , subsWithNA=T)
   pathogen_grid <- pathogen_pp_raster * pop_grid
-  # Why is this twice?
+
   temp<-data.frame(bla=NA,value=0)
-  pathogen_grid <- subs(pathogen_urban_water_grid,temp,subsWithNA=F)
-  pathogen_grid <- subs(pathogen_urban_water_grid,temp,subsWithNA=F)
+  pathogen_grid <- subs(pathogen_grid,temp,subsWithNA=F)
   return(pathogen_grid)
 }
 
diff --git a/Model scripts/Pathogenflows.R b/Model scripts/Pathogenflows.R
index c04b743..88d0676 100644
--- a/Model scripts/Pathogenflows.R	
+++ b/Model scripts/Pathogenflows.R	
@@ -2,20 +2,113 @@
 library(pathogenflows)
 
 pathogenflow.model.run <- function(input_data,pathogen){
-  emissions <- data.frame()
+  emissions <- data.frame(iso=input_data$iso,subarea=input_data$subarea)
   # 1) extract pathogen type from pathogen
   pathogenType <- "Virus"
-  # 2a) prepare data for urban and rural
-  
-  # 2b) call loadings function twice for urban and rural areas
-  onsite_urban = data.frame()
-  onsite_rural = data.frame()
-  rural_loadings <- pathogenflows::getLoadings(onsite_rural,pathogenType)
-  urban_loadings <- pathogenflows::getLoadings(onsite_urban,pathogenType)
-  # 2c) postprocess output
+  area_types <- c("urban","rural")
+  human_age_types <- c("child","adult")
+  all_loadings <- list()
+  # 2) call loadings function 4 times for urban-adult/child prevalence and rural-adult/child prevalence
+  for(area_type in area_types){
+    for(human_age_type in human_age_types ){
+      #2a) prepare data for urban and rural
+      onsite_data <- pathogenflow.model.get.input(input_data,area_type,human_age_type)
+      file_out <- file.path(working.path.in,sprintf("onsite_%s_%s.csv",area_type,human_age_type))
+      write.csv(onsite_data,file = file_out)
+      #call loadings
+      loadings <- pathogenflows::getLoadings(file_out,pathogenType)
+      # store in list
+      all_loadings[[area_type]][[human_age_type]] <- loadings
+    }
+  }
+
+  # 3) postprocess output
+  # add child  + adult loadings for sanitation and areas
+  # apply to area types urban/rural
+  for(area_type in names(all_loadings)){
+    area_loadings <- all_loadings[[area_type]]
+    adult <- area_loadings$adult$detailed
+    child <- area_loadings$child$detailed
+    # apply to all sub areas
+    for(i in 1:length(emissions$subarea)){
+      sub_area <- emissions$subarea[i]
+      ncols <- ncol(adult[[sub_area]])
+      sanitation_types <- adult[[sub_area]]$id
+      emissions_sanitation <- adult[[sub_area]][3:ncols] + child[[sub_area]][3:ncols]
+      
+      # apply to all sanitation types
+      for(j in 1:length(sanitation_types)){
+        sanitation_type <- sanitation_types[j]
+        emissions_col_name_postfix <- sprintf("%s_%s",sanitation_type,area_type)
+        emissions_col_to_surface <- sprintf("%s_%s","to_surface",emissions_col_name_postfix)
+        emissions_col_sewerage <- sprintf("%s_%s","sewerage",emissions_col_name_postfix)
+        to_surface <- emissions_sanitation$toSurface[j]
+        to_sewerage <- emissions_sanitation$sewerage[j]
+        # store only to_surface and sewerage
+        emissions[[emissions_col_to_surface]][i] <- to_surface
+        emissions[[emissions_col_sewerage]][i] <- to_sewerage
+      }
+    }
+  }
   
   # 3) create similair emissions data.frame as original human emissions
   
   # 4) return emissions
   return(emissions)
+}
+
+pathogenflow.model.get.input <- function(input,area_type,human_type){
+  #@Nynke: population is population per human type (adult/child) or total population?
+  col_names <- c("iso","scenario","region","population","sheddingRate","prevalence","flushSewer","flushSeptic","flushPit","flushOpen",
+                          "flushUnknown","pitSlab","pitNoSlab","compostingTwinSlab","compostingTwinNoSlab","compostingToilet","bucketLatrine",
+                          "containerBased","hangingToilet","openDefecation","other","isShared","sewerLeak","emptied","isWatertight","hasLeach","coverBury",	
+                          "emptiedTreatment","emptyFrequency","pitAdditive","flushElsewhere","pitVIP","pitTraditional","otherLatrine","otherImproved",
+                          "otherUnimproved","dontKnow","pitLined", "pitUnlined") 
+  onsite_data <- data.frame(matrix(ncol=length(col_names),nrow = dim(input)[1]))
+  colnames(onsite_data) <- col_names
+  onsite_data$iso <- input$iso
+  onsite_data$region <- input$subarea
+  onsite_data$sheddingRate <- input$sheddingRate
+
+  input_col_names <- colnames(input)
+  if(area_type=="urban"){
+    postfix <- "urb"
+    onsite_data$population <- input$population*input$fraction_urban_pop
+  }
+  else if(area_type=="rural"){
+    postfix <- "rur"
+    onsite_data$population <- input$population* (1- input$fraction_urban_pop)
+  }
+  # filter for urban and rural
+  filtered_cols <- input_col_names[str_ends(input_col_names,postfix)]
+  for(input_colname in filtered_cols){
+    # get corresponding colname for onsite_data
+    onsite_colname <- str_replace(input_colname,sprintf("_%s",postfix),"")
+    # check if it part of onsite_data
+    if(onsite_colname %in% col_names){
+      onsite_data[[onsite_colname]] <- input[[input_colname]]
+    }
+  }
+  incidence = 0
+  if(human_type == "child"){
+    onsite_data$population <- input$population * input$fraction_pop_under5
+    if(area_type == "urban"){
+      incidence <- input$incidence_urban_under5
+    }
+    else if(area_type == "rural"){
+      incidence <- input$incidence_rural_under5 
+    }
+  }
+  else if(human_type == "adult"){
+    onsite_data$population <- input$population * (1-input$fraction_pop_under5)
+    if(area_type == "urban"){
+      incidence <- input$incidence_urban_5plus
+    }
+    else if(area_type == "rural"){
+      incidence <- input$incidence_rural_5plus 
+    }
+  }
+  # @Nynke: Is this ok?
+  onsite_data$prevalence <- incidence*input$sheddingRate*input$shedding_duration
+  return(onsite_data)
 }
\ No newline at end of file
-- 
GitLab