Finalized the grassland approach with current and potential stocking rates.

Also started to compare the area extracted (by Leandro and Me) with the EUROSTAT and FAOSTAT values.
Found a major bug: country_map seems to be not matching the country codes from the FAO sheet. @WJS: check again what is going wrong here!

Grassland/Grassland.R

@@ -110,42 +110,6 @@ dat_grass_crop =
   #               FM_ABG_grass_kgha_adj = FM_ABG_grass_kgha * stock_rate_fct)
 
 
-# Subsetting per stocking rate --------------------------------------------
-grassing_density_potential = 
-  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: 
-# https://ec.europa.eu/eurostat/statistics-explained/images/1/1a/Livestock_patterns_statistics_2016.xlsx
-# The grazing livestock density index measures the stock of grazing animals 
-# (cattle, sheep, goats and equidae) expressed in livestock units (LSU) per hectare of fodder area. 
-
-grassing_density_current = read_csv("Input_data/Garzing_density_2016_EUROSTAT.csv") %>%
-  clean_names() %>% 
-  
-    # when grazing intensity is larger 2, it is set to 2. This is done
-    # because we have no more than grazing density of 2 in our grassland data from LPJML
-  mutate(adj_grazing_dens  = case_when(grazing_livestock_density > 2 ~ 2,
-                                       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. 
-  # 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
-                adj_grazing_dens = as.character(adj_grazing_dens),
-                adj_grazing_dens  = case_when(adj_grazing_dens == "0" ~ "00",
-                                              adj_grazing_dens == "2" ~ "02",
-                                              adj_grazing_dens == "4" ~ "04",
-                                              adj_grazing_dens == "6" ~ "06",
-                                              adj_grazing_dens == "8" ~ "08",
-                                        TRUE ~ adj_grazing_dens),
-                scenario = "current") 
-
 
 # Available land per lut --------------------------------------------------
 
@@ -155,6 +119,10 @@ available_land_area = dat_grass_crop %>%
   dplyr::select(-stock_rate) %>% 
   left_join(., country_map %>% mutate(ADM0_CODE_GAUL = as.character(ADM0_CODE_GAUL)), by = c("adm0_cod" = "ADM0_CODE_GAUL"))
 
+
+
+
+
 # EU28 - Subsetting -------------------------------------------------------
 EU_iso3 = c("AUT","BEL", "BGR", "HRV", "CZE", "DNK", "EST","FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA","LVA", 
             "LTU", "LUX", "MLT", "NLD", "POL", "PRT","ROU", "SVK","SVN", "ESP", "SWE", "GBR", "CYP")
@@ -202,7 +170,8 @@ write_csv(lut_area_av_pasture_arable, "Input_data/lut_area_av_pasture_arable.csv
 
 # Adding iso3  ------------------------------------------------------------
 # Subset EU28 - ADding yields, renaming, mapping names
-dat_grass = dat_grass_crop %>% dplyr::filter(lut_cifos != "cropland") %>% 
+dat_grass = dat_grass_crop %>% 
+  dplyr::filter(lut_cifos != "cropland") %>% 
   mutate(area_ha = area_ha *100) #conversion from km2 to ha
 
 # Joining countries
@@ -213,8 +182,7 @@ dat_grass_cntr = dat_grass %>%
 EU_iso3 = c("AUT","BEL", "BGR", "HRV", "CZE", "DNK", "EST","FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA","LVA", 
             "LTU", "LUX", "MLT", "NLD", "POL", "PRT","ROU", "SVK","SVN", "ESP", "SWE", "GBR", "CYP")
 
-grassing_density_potential_cntr = grassing_density_potential %>% left_join(., country_map %>% mutate(ADM0_CODE_GAUL = as.character(ADM0_CODE_GAUL)), by = c("adm0_cod" = "ADM0_CODE_GAUL")) %>% 
-  filter(., iso3_SPAM  %in% EU_iso3)
+
 
 grass_EU28 = dat_grass_cntr %>% filter(., iso3_SPAM  %in% EU_iso3)
 
@@ -245,9 +213,69 @@ grass_EU28_final =
     relocate(lut_cifos, .before = prod_level) %>% 
     relocate(yield_kgha, .before = area_ha)
   
-grass_current = grass_EU28_final %>% right_join(grassing_density_current %>% 
+
+
+
+# Adding stocking rate current and potential ------------------------------
+
+# Subsetting per stocking rate --------------------------------------------
+grassing_density_potential_init = 
+  dat_grass_crop %>%
+  dplyr::group_by(adm0_cod, aez, soil, lut_cifos) %>% #selecting the stocking rate per zone with the maximum yield.
+  dplyr::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
+
+dat_pot_zero = grassing_density_potential_init %>% filter(any(yield_DM_kgha == 0 & 
+                                      area_ha == 0)) %>% 
+  dplyr::top_n(-1, stock_rate)
+
+dat_pot_nonzero = grassing_density_potential_init %>% 
+  anti_join(dat_pot_zero, by=c("adm0_cod", "soil", "aez","stock_rate", "lut_cifos")) %>% 
+  dplyr::top_n(1, stock_rate)
+
+grassing_density_potential = bind_rows(dat_pot_zero,dat_pot_nonzero) %>% 
+  dplyr::top_n(-1, stock_rate)
+
+# Classification of coutnries
+grassing_density_potential_cntr = grassing_density_potential %>% left_join(., country_map %>% mutate(ADM0_CODE_GAUL = as.character(ADM0_CODE_GAUL)), by = c("adm0_cod" = "ADM0_CODE_GAUL")) %>% 
+  filter(., iso3_SPAM  %in% EU_iso3)
+  
+
+# Current stocking rate  --------------------------------------------------
+# For current grazing intensity we take country values from EUROSTAT: 
+# https://ec.europa.eu/eurostat/statistics-explained/images/1/1a/Livestock_patterns_statistics_2016.xlsx
+# The grazing livestock density index measures the stock of grazing animals 
+# (cattle, sheep, goats and equidae) expressed in livestock units (LSU) per hectare of fodder area. 
+
+grassing_density_current = read_csv("Input_data/Garzing_density_2016_EUROSTAT.csv") %>%
+  clean_names() %>% 
+  
+  # when grazing intensity is larger 2, it is set to 2. This is done
+  # because we have no more than grazing density of 2 in our grassland data from LPJML
+  mutate(adj_grazing_dens  = case_when(grazing_livestock_density > 2 ~ 2,
+                                       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. 
+  # 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
+                adj_grazing_dens = as.character(adj_grazing_dens),
+                adj_grazing_dens  = case_when(adj_grazing_dens == "0" ~ "00",
+                                              adj_grazing_dens == "2" ~ "02",
+                                              adj_grazing_dens == "4" ~ "04",
+                                              adj_grazing_dens == "6" ~ "06",
+                                              adj_grazing_dens == "8" ~ "08",
+                                              TRUE ~ adj_grazing_dens),
+                scenario = "current")
+
+
+
+
+grass_current = grass_EU28_final%>% ungroup() %>% right_join(grassing_density_current %>% 
                                                   select(-grazing_livestock_density), 
-                                                by =   c("prod_level"="adj_grazing_dens", "country_cifos"="eu_28"))%>% 
+                                                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) %>% 
@@ -267,16 +295,148 @@ grass_potential = grass_EU28_final %>% ungroup() %>%  right_join(grassing_densit
   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")
+# Checking the data
+# Country
+grass_potential %>% distinct(country_cifos)
+grass_current %>% distinct(country_cifos)
+
+
+write_csv(grass_final_EU28_stockrate_adj, "Input_data/grass_EU28_final.csv")
 grass_EU28_final = read_csv("Input_data/grass_EU28_final.csv")
 
 
+
+
+# Comparing area with EUROSTAT --------------------------------------------
+# Function to extract data from EUROSTAT sheets (source: https://ec.europa.eu/eurostat/databrowser/view/tag00025/default/table?lang=en)
+# https://ec.europa.eu/eurostat/web/main/data/database
+Fun_Eurostat_area = function(Eurostat_data){
+  Eurostat_data %>% dplyr::slice(-c(1:5)) %>% rename("country" = "TIME") %>% 
+    dplyr::select(-contains("...")) %>% clean_names() %>% 
+      dplyr::filter(!country %in% c(":", "e", "Available flags:","p","Special value"),
+                    !is.na(country)) %>% 
+    pivot_longer(cols = "x2009":"x2020", 
+                 names_to = "year",
+                 values_to = "area_ha") %>% 
+  dplyr::mutate(area_ha=as.numeric(area_ha),
+         area_ha=area_ha*1000) %>% 
+  dplyr::group_by(country) %>% 
+  dplyr::summarise(area_ha = mean(area_ha, na.rm = TRUE)) %>% 
+    mutate(country = recode(country, "Germany (until 1990 former territory of the FRG)"="Germany"))
+}
+
+# Loading different croptypes
+eurostat_area_grass = read_excel("Input_data/EUROSTAT_area_agri_utilisiation.xlsx", sheet = "Sheet 3", skip = 8)
+eurostat_area_arable = read_excel("Input_data/EUROSTAT_area_agri_utilisiation.xlsx", sheet = "Sheet 2", skip = 8)
+eurostat_area_permanentcrp = read_excel("Input_data/EUROSTAT_area_agri_utilisiation.xlsx", sheet = "Sheet 4", skip = 8)
+
+# Running function on all dfs
+grass_area = Fun_Eurostat_area(eurostat_area_grass)
+arable_area = Fun_Eurostat_area(eurostat_area_arable)
+permanent_area = Fun_Eurostat_area(eurostat_area_permanentcrp)
+
+# join it to our data and calculate the relative difference in percentage
+Fun_area_join = function(dat, vec_lut){
+  lut_area_av %>% 
+  dplyr::filter(lut_cifos %in% vec_lut) %>%
+  # group_by(country_cifos,aez_code, soil_cifos, lut_cifos) %>%
+  # summarise(area_ha_cifos = max(area_ha)) %>% #summarizing over stocking rates (all same area)
+  # ungroup() %>%
+  group_by(country_cifos) %>% 
+  summarise(area_ha_cifos = sum(area_ha)) %>%
+  left_join(dat, by=c("country_cifos"="country")) %>% 
+  mutate(diff_perc = (100/area_ha*area_ha_cifos)-100)
+}
+
+grass_area_cifos = Fun_area_join(dat = grass_area, vec_lut = c("pasture_marginal", "rangeland"))
+arable_area_cifos =  Fun_area_join(dat = arable_area, 
+                                   vec_lut = c("arable"))
+
+# Visualising the deviation results
+Fun_ggplot_area_deviation = function(dat,lut_vec,author_vec,vec_source){
+  dat$dev_type <- ifelse(dat$diff_perc < 0, "below", "above")  # above / below avg flag
+  ggplot(dat, aes(x=reorder(country_cifos, -diff_perc), y=diff_perc, label=diff_perc)) + 
+    geom_bar(stat='identity', aes(fill=dev_type), width=.5)  +
+    scale_fill_manual(name="Deviation", 
+                      labels = c("Above Eurostat", "Below Eurostat"), 
+                      values = c("above"="#00ba38", "below"="#f8766d")) + 
+    labs(title= paste0("Area ",lut_vec, " comparison - " , author_vec, " extraction"),
+         subtitle = paste0("Reference source: ", vec_source)) + 
+    coord_flip() + 
+    theme_bw()
+}
+
+# Plotting - leandro extraction
+Fun_ggplot_area_deviation(arable_area_cifos, lut_vec="arable", author_vec = "WJS", vec_source = "EUROSTAT")
+Fun_ggplot_area_deviation(grass_area_cifos,lut_vec="grass", author_vec = "WJS" , vec_source = "EUROSTAT")
+
+# Comparing the own extraction from R script 'Zones_AEZ.R' -----------------
+Fun_area_join_WJSextract = function(dat,vec_lut){
+Comb_Zon_EU28 %>% select(-c(ADM0_CODE,iso3_SPAM)) %>% 
+  rename(country_cifos =  name_cntr_CIFOS2021,
+         soil_cifos = Soil_agri,
+         aez_code = AEZ, 
+         arable = crop_areaha, 
+         pasture_arable = pasture_areaha, 
+         pasture_marginal = pastureCrop_areaha, 
+         rangeland = rangeland_areaha) %>% 
+  pivot_longer(cols = arable:rangeland,
+               names_to = "lut_cifos", 
+               values_to = "area_ha") %>% 
+  relocate(soil_cifos, .before = lut_cifos) %>% 
+    dplyr::filter(lut_cifos %in% vec_lut) %>%
+    group_by(country_cifos) %>% 
+    summarise(area_ha_cifos = sum(area_ha)) %>%
+    left_join(dat, by=c("country_cifos"="country")) %>% 
+    mutate(diff_perc = (100/area_ha*area_ha_cifos)-100)
+  }
+
+# Final comparison tables for grass and arable land
+grass_area_cifos_WJex = Fun_area_join_WJSextract(dat = grass_area, vec_lut = c("pasture_marginal", "rangeland", "pasture_arable"))
+arable_area_cifos_WJex =  Fun_area_join_WJSextract(dat = permanent_area %>% bind_rows(arable_area) %>% 
+                                                group_by(country) %>% summarise(area_ha = sum(area_ha)), 
+                                              vec_lut = c("arable"))
+
+# Plotting - WJS extraction
+Fun_ggplot_area_deviation(grass_area_cifos_WJex, lut_vec="arable", author_vec = "WJS",vec_source = "EUROSTAT")
+Fun_ggplot_area_deviation(arable_area_cifos_WJex,lut_vec="grass", author_vec = "WJS", vec_source = "EUROSTAT")
+
+
+# comparing with FAO area ---------------------------------------------------------------------
+fao_arable_area =
+  read_csv("INput_data/area_harvested_FAO_arablecrops.csv") %>% 
+  clean_names() %>% type.convert()  %>% 
+  left_join(country_map, by = c("area_code_fao" = "ADM0_CODE_GAUL")) #%>% 
+  filter(., iso3_SPAM  %in% EU_iso3) %>% 
+  rename("country" = "name_cntr_CIFOS2021",
+         "area_ha"="value") %>%
+  group_by(year_code, country) %>%
+  summarise(area_ha = mean(area_ha,na.rm=T)) %>% 
+  group_by(country) %>%
+  summarise(area_ha = sum(area_ha,na.rm=T)) 
+
+grass_area_cifos_WJex_fao = Fun_area_join_WJSextract(dat = fao_arable_area, vec_lut = c("arable"))
+
+
+Comb_Zon_EU28 %>%
+  left_join(fao_arable_area, by=c("country_cifos"="country")) %>% 
+  mutate(diff_perc = 100-(100/area_ha*area_ha_cifos))
+
+# EU28 - Subsetting -------------------------------------------------------
+EU_iso3 = c("AUT","BEL", "BGR", "HRV", "CZE", "DNK", "EST","FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA","LVA", 
+            "LTU", "LUX", "MLT", "NLD", "POL", "PRT","ROU", "SVK","SVN", "ESP", "SWE", "GBR", "CYP")
+
+
+
+grass_EU28 = dat_grass_cntr %>% filter(., iso3_SPAM  %in% EU_iso3)
+
+
+
 # Aggreagted by stocking rate ---------------------------------------------
 #idea: to aggregate by stocking rate weigthed mean by stocking rate 08,10,12
 

R_to_SQL_database/Cifos_Model_Data_EU.R

@@ -2938,3 +2938,5 @@ CropNutr$crop_cifos[!CropNutr$crop_cifos %in% procout_raw_dist$Product] #charact
 
 
 
+
+

Zones/Zones_AEZ.R

@@ -266,7 +266,7 @@ Comb_Zon_EU28 = read_csv(here::here("Input_data", "AEZ_ADM0_SOILagri_EU28.csv"))
 # Checking the area per country
 adm0_check = shape_gaul %>% dplyr::select(ADM0_CODE, ADM0_NAME) %>% distinct_at(vars(ADM0_CODE, ADM0_NAME))
 # 
-Chk_country = Comb_Zon %>% 
+Chk_country = Comb_Zon_EU28 %>% 
   # st_set_geometry(NULL) %>% 
   as_tibble() %>%
   group_by(ADM0_CODE) %>%