Fertilization rates are done

Fertilization_rates/Fertilizatio_rates.R

@@ -106,70 +106,83 @@ crops_NP_map = crops_NP %>% left_join(., SPAM2010_CropList, by="short_spam2010")
 # --> DONE in excel
 
 # Calculate total NP content per crop  - NP content DM * production
-# Calculating the crop requirements (country, crop, all intensity levels combined (A)) > assign the FAO consumption to the crops based on requirements
+# Calculating the crop requirements (country, crop, all intensity levels combined (A))
 NP_content = crops_NP_map %>% 
   right_join(., SPAM_EU28_dat, by=c("short_spam2010"="crops_spam")) %>%
-  mutate(Nutr_req_harv_N = Ntotal/100*(DM_perc)*Production_mt,
-         Nutr_req_harv_P = Ptotal/100*(DM_perc)*Production_mt) %>% 
-  # rename(Country = Country.x) %>% 
-  dplyr::select(short_spam2010, code_fao, iso3, country_cifos, peat, HarvestedArea_ha, Production_mt, PhysicalArea_ha, Yield_kgha ,  Nutr_req_harv_N , Nutr_req_harv_P)
+  mutate(Nutr_req_harv_N = (Ntotal/100)*(DM_perc/100)*Production_mt, #NP removal from production
+         Nutr_req_harv_P = (Ptotal/100)*(DM_perc/100)*Production_mt) %>% 
+  mutate(N_removal_kgha = Nutr_req_harv_N*1000/HarvestedArea_ha, #NP removal (kg/ha) from yields (transforming tonnes to kg)
+         P_removal_kgha = Nutr_req_harv_P*1000/HarvestedArea_ha) %>% 
+  mutate_at(vars(N_removal_kgha,P_removal_kgha), ~replace(., is.nan(.), 0))%>% 
+  dplyr::select(short_spam2010, code_fao, iso3, country_cifos, peat, HarvestedArea_ha, Production_mt, PhysicalArea_ha, Yield_kgha ,  Nutr_req_harv_N , Nutr_req_harv_P, N_removal_kgha, P_removal_kgha)
+
+# ggplot(NP_content, aes(y=N_removal_kgha))+
+#   geom_histogram()
          
 # Ratio of crops per defined zones (here per crop and country) ------------
 NP_content_ratio = NP_content %>%  
-  pivot_longer(cols= "Nutr_req_harv_N":"Nutr_req_harv_P", names_to = "nutr_req") %>%
+  pivot_longer(cols= "Nutr_req_harv_N":"Nutr_req_harv_P", names_to = "nutr_req_production") %>%
+  # pivot_longer(cols= "N_removal_kgha":"P_removal_kgha", names_to = "nutr_req_yield", names_repair = "unique")
   mutate_at(vars(value), ~replace(., is.nan(.), 0))%>% 
-  group_by(iso3, nutr_req) %>% 
-  mutate(perc = value/sum(value)) %>%
-  mutate(label = case_when(stringi::stri_detect_fixed(nutr_req, "_N") ~ "N",
-                           stringi::stri_detect_fixed(nutr_req, "_P")  ~ "P"))
+  group_by(iso3, nutr_req_production) %>% 
+  mutate(NP_content_ratio = value/sum(value)) %>%
+  # mutate(NP_content_ratio_yld = value...13/sum(value...13)) %>%
+  mutate(label = case_when(stringi::stri_detect_fixed(nutr_req_production, "_N") ~ "N",
+                           stringi::stri_detect_fixed(nutr_req_production, "_P")  ~ "P"))
 
 # eval - all correct (28.04.2021)
 # val=NP_content_ratio %>% group_by(iso3, nutr_req) %>% 
 #   summarize(val = sum(perc))
 
-# FAO NP consumption data: Creating three year average and cleaning variables
-avg_yr_all_EU28 = fertilizer_FAO_EU28 %>%
-  group_by(Area, Element, Domain, Item) %>%
-  mutate(Value = mean(Value)) %>%  #creating three year average
-  filter(Year==2010) %>%
-  dplyr::select(-c(`Domain Code`, `Year Code`)) %>% 
-    mutate(Value = case_when(Item == "Nutrient phosphate P2O5 (total)" ~ Value/2.29, #converting P2O5 to P > conversion facrotr 2.29 https://www.researchgate.net/post/How-to-calculate-phosphorous-from-P2O5-which-has-been-recommended-at-the-standard-dose-of-45-kg-ha, 
-                             Item == "Nutrient nitrogen N (total)" ~ Value),
-           Item = recode(Item,'Nutrient phosphate P2O5 (total)' =  "Nutrient phosphate P (total)"))
-  
+# # FAO NP consumption data: Creating three year average and cleaning variables
+# avg_yr_all_EU28 = fertilizer_FAO_EU28 %>%
+#   group_by(Area, Element, Domain, Item) %>%
+#   mutate(Value = mean(Value)) %>%  #creating three year average
+#   filter(Year==2010) %>%
+#   dplyr::select(-c(`Domain Code`, `Year Code`)) %>% 
+#     mutate(Value = case_when(Item == "Nutrient phosphate P2O5 (total)" ~ Value/2.29, #converting P2O5 to P > conversion facrotr 2.29 https://www.researchgate.net/post/How-to-calculate-phosphorous-from-P2O5-which-has-been-recommended-at-the-standard-dose-of-45-kg-ha, https://plants.usda.gov/npk/fertilizer_equivalents_def.html 
+#                              Item == "Nutrient nitrogen N (total)" ~ Value),
+#            Item = recode(Item,'Nutrient phosphate P2O5 (total)' =  "Nutrient phosphate P (total)"))
+#   
 
 nutrients_FAO_2010_EU28 = fertilizer_FAO_EU28 %>%
   group_by(Area, Element, Domain, Item) %>%
   mutate(Value = mean(Value)) %>%  #creating three year average
   filter(Year==2010) %>%
   # dplyr::select(-c(`Domain Code`, `Year Code`)) %>% 
-  mutate(Value = case_when(Item == "Nutrient phosphate P2O5 (total)" ~ Value/2.29, #converting P2O5 to P > conversion facrotr 2.29 https://www.researchgate.net/post/How-to-calculate-phosphorous-from-P2O5-which-has-been-recommended-at-the-standard-dose-of-45-kg-ha, 
+  mutate(Value = case_when(Item == "Nutrient phosphate P2O5 (total)" ~ Value/2.29, #converting P2O5 to P > conversion facrotr 2.29 https://www.researchgate.net/post/How-to-calculate-phosphorous-from-P2O5-which-has-been-recommended-at-the-standard-dose-of-45-kg-ha, https://plants.usda.gov/npk/fertilizer_equivalents_def.html
                            Item == "Nutrient nitrogen N (total)" ~ Value),
          Item = recode(Item,'Nutrient phosphate P2O5 (total)' =  "Nutrient phosphate P (total)")) %>% 
   dplyr::select(Area, iso3_code_FAO, `Area Code`, Domain, Item, `Item Code`, Year, Value)
 
-# Identify the gap between FAO and Requirements
+# Adding NP labels, multiplying FAO consumption with crop share ratios and Identify the gap between FAO and Requirements
 gapFAO_SPAM = nutrients_FAO_2010_EU28 %>% 
   mutate(label = case_when(stringi::stri_detect_fixed(Item, "nitrogen") ~ "N",
                            stringi::stri_detect_fixed(Item, "phosphate")  ~ "P")) %>% 
-  # pivot_wider(names_from = Item, values_from = Value)
-  
-  # dplyr::select(-c(Domain, `Area Code`, Item))# %>%
   right_join(.,NP_content_ratio, by=c("iso3_code_FAO"="iso3", "label"="label")) %>%
-  arrange(Area)
-  # filter(Element == "Production") %>%
-  # distinct_at(vars(Area, short_spam2010), .keep_all = T) %>%
-    mutate(N_gap = Value-Nutr_req_harv_N,
-           P_gap = Value-Nutr_req_harv_P)
-
-
-# What is with negative gaps?
-
+  mutate(FAO_nutr = NP_content_ratio *Value) %>% 
+  mutate(nutr_gap = case_when(value > 0 ~ FAO_nutr-value, # We look at the gap btw the ratio of fao and the ratio of crop-contetnt related fertilizer use
+                              TRUE ~ 0)) %>% 
+  mutate(fert_rate_FAO_SPAM = case_when(nutr_gap > 0 & label == "N" ~ nutr_gap/HarvestedArea_ha+N_removal_kgha,
+                                        nutr_gap > 0 & label == "P" ~ nutr_gap/HarvestedArea_ha+P_removal_kgha,
+                                        nutr_gap <= 0 & label == "P" ~ nutr_gap/HarvestedArea_ha+P_removal_kgha,
+                                        nutr_gap <= 0  & label == "N" ~ nutr_gap/HarvestedArea_ha+N_removal_kgha)) %>% 
+  mutate_at(vars(fert_rate_FAO_SPAM), ~replace(., is.nan(.), 0)) z
 
+mean(gapFAO_SPAM$fert_rate_FAO_SPAM)
 
+# gapFAO_SPAM %>%
+#   ggplot(aes(y=fert_rate))+
+#   geom_histogram(bins = 5)
 
 # Calculate the ratio of each crop per Zone, country based on requirements (Allocation factor based) > Scaling
 # Ratio from NP values of the crops in EU27
+NP_ratio = NP_content_ratio %>% dplyr::select(iso3, peat, label, short_spam2010, NP_content_ratio)
+
+fert_rates_alloc = gapFAO_SPAM %>% 
+  mutate(FAO_fert_alloc = NP_content_ratio*Value)
+  
+  
 
 
 # allocate the amount of NP fertilizer (GAP between requirement and FAO) to the crops based on the allocation factor based on point 2