finalized the grazing grassland approach where the prod_level is now selected...

finalized the grazing grassland approach where the prod_level is now selected based on the potential and the current grazing intensity. Current is based on EUROSTAT grazing intensity. While the maximum is based on the stocking rate with the highest yields  per zone.

Grassland/Grassland.R

@@ -112,12 +112,10 @@ dat_grass_crop =
 
 # Subsetting per stocking rate --------------------------------------------
 grassing_density_potential = 
-  dat_grass_crop %>%  
-  #selecting the stocking rate per zone with the maximum yield. 
-  group_by(adm0_cod, aez, soil, lut_cifos) %>% 
-  top_n(1, yield_DM_kgha)   
-
-
+  dat_grass_crop %>%
+  group_by(adm0_cod, aez, soil, lut_cifos) %>% #selecting the stocking rate per zone with the maximum yield.
+  top_n(1, yield_DM_kgha) %>%  #selects the  LSU (=grazing density) with the highest yields per group combination of adm0_cod, aez, soil, lut_cifos.
+  mutate(scenario = "potential") #tag it as potential to be able to distinguish it from the current later after rbind
 
 # Current stocking rate  --------------------------------------------------
 # For current grazing intensity we take country values from EUROSTAT: 
@@ -134,7 +132,8 @@ grassing_density_current = read_csv("Input_data/Garzing_density_2016_EUROSTAT.cs
                                        TRUE ~ grazing_livestock_density),
          adj_grazing_dens = adj_grazing_dens*10) %>% 
   
-  # If the stocking rate is uneven, we add 0.1 LSU to make it even and to match the LPJML
+  # If the stocking rate is uneven, we add 0.1 LSU to make it even and to match the LPJML. 
+  # This is done as we underestimate the current grazing density in all the countries with LSU > 2 to 2. 
   dplyr::mutate(adj_grazing_dens= case_when(adj_grazing_dens %% 2 != 0 ~ adj_grazing_dens+1,
                                             TRUE~adj_grazing_dens),
                 # Some formatting
@@ -144,7 +143,9 @@ grassing_density_current = read_csv("Input_data/Garzing_density_2016_EUROSTAT.cs
                                               adj_grazing_dens == "4" ~ "04",
                                               adj_grazing_dens == "6" ~ "06",
                                               adj_grazing_dens == "8" ~ "08",
-                                        TRUE ~ adj_grazing_dens))
+                                        TRUE ~ adj_grazing_dens),
+                scenario = "current") 
+
 
 # Available land per lut --------------------------------------------------
 
@@ -242,27 +243,34 @@ grass_EU28_final =
     relocate(crop_cifos, .after = lut_cifos) %>%
     relocate(country_cifos, .before = aez_code) %>% 
     relocate(lut_cifos, .before = prod_level) %>% 
-    relocate(yield_kgha, .before = area_ha) %>%
-    left_join(grassing_density_current, by =   c("prod_level"="adj_grazing_dens", "country_cifos"="eu_28")) %>% 
-    left_join(grassing_density_potential_cntr, by =   c("prod_level"="stock_rate", 
+    relocate(yield_kgha, .before = area_ha)
+  
+grass_current = grass_EU28_final %>% right_join(grassing_density_current %>% 
+                                                  select(-grazing_livestock_density), 
+                                                by =   c("prod_level"="adj_grazing_dens", "country_cifos"="eu_28"))%>% 
+  filter(!crop_cifos == "grass_other") %>% 
+  mutate(prod_level=paste0("cur_",prod_level)) %>% 
+  select(-scenario) %>% 
+  # THere is yields where there is no area. Then yields are set to 0
+  dplyr::mutate(yield_kgha = case_when(area_ha == 0 & yield_kgha > 0 ~ 0,
+                                       TRUE ~ yield_kgha))
+
+grass_potential = grass_EU28_final %>% ungroup() %>%  right_join(grassing_density_potential_cntr %>% ungroup() %>%
+                                                    select(-c(iso3_SPAM, area_ha,yield_FM_kgha,yield_DM_kgha,adm0_cod)), 
+                                                           by = c("prod_level"="stock_rate", 
                                                         "country_cifos"="name_cntr_CIFOS2021",
                                                         "aez_code"="aez", 
                                                         "soil_cifos"="soil",
-                                                        "lut_cifos"="lut_cifos")) %>% 
-  filter(!is.na(grazing_livestock_density) | !is.na(adm0_cod))
-  
-    
-  
-  
-    dplyr::filter(prod_level == "10") %>%
-    mutate(prod_level = "A") %>% 
-           #scenario = "current") %>% 
-  # relocate(scenario, .before = crop_cifos) %>% 
-  ungroup() %>% 
-  # dplyr::select(-adm0_cod) %>% 
+                                                        "lut_cifos")) %>% 
+  filter(!crop_cifos == "grass_other") %>% 
+  mutate(prod_level=paste0("pot_",prod_level)) %>% 
+  select(-scenario) %>% 
   # THere is yields where there is no area. Then yields are set to 0
   dplyr::mutate(yield_kgha = case_when(area_ha == 0 & yield_kgha > 0 ~ 0,
-                                TRUE ~ yield_kgha))
+                                       TRUE ~ yield_kgha))
+
+                                                        
+grass_final_EU28_stockrate_adj = bind_rows(grass_potential, grass_current)
     
 
 write_csv(grass_EU28_final, "Input_data/grass_EU28_final.csv")