Commit db456484 authored by Simon, Wolfram's avatar Simon, Wolfram
Browse files

In model data script: worked on newly discussed Import export approach where...

In model data script: worked on newly discussed Import export approach where we substrat the non food and stocks from the import to get more realistic results .In progress.
In grassland: the stocking rate approach is implemented where we have potential (highest) and realistic stocking rates expressed in LSU/ha (0-2)
parent 89c487e6
......@@ -111,15 +111,46 @@ dat_grass_crop =
# Subsetting per stocking rate --------------------------------------------
dat_grass_crop %>% group_by(adm0_cod, soil, aez, lut_cifos) %>% top_n(1, yield_DM_kgha) #selecting the stocking rate per zone with the maximum yield.
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)
# 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
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))
# Available land per lut --------------------------------------------------
available_land_area = dat_grass_crop %>%
mutate(area_ha = area_ha *100) %>%
filter(stock_rate=="00") %>%
filter(stock_rate == "00") %>%
dplyr::select(-stock_rate) %>%
left_join(., country_map %>% mutate(ADM0_CODE_GAUL = as.character(ADM0_CODE_GAUL)), by = c("adm0_cod" = "ADM0_CODE_GAUL"))
......@@ -181,10 +212,12 @@ 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)
grass_EU28_final =
grass_EU28_final =
grass_EU28 %>%
# group_by(soil,aez, stock_rate, lut_cifos) %>%
# dplyr::summarise(yield_FM_kgha = weighted.mean(x=yield_FM_kgha,
......@@ -209,7 +242,18 @@ 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) %>%
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",
"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") %>%
......
......@@ -2431,35 +2431,42 @@ write_csv(Baseline_food_supply_FoodGroups, "Input_data/Baseline_food_supply_FG.c
# Problem: Protein per day is too little! Maybe the mapping went wrong. I try again with a item list file that I download from FAO directly. https://www.fao.org/faostat/en/#data/FBS
# Latest update (21.2.22) FBS Baseline ----------------------------------------------
# The numbers here do not make fully sense
FBS_FAO_downloaded = read_csv("Input_data/FoodBalanceSheet.csv")
FBS_FAO_downloaded_EUgroup = read_csv("Input_data/FAOSTAT_data_2-9-2022.csv") #this sheet has a bit more options
# FBS_FAO_downloaded = read_csv("Input_data/FoodBalanceSheet.csv")
# FBS_FAO_downloaded_EUgroup = read_csv("Input_data/FAOSTAT_data_2-9-2022.csv") #this sheet has a bit more options
FBS_FAO_downloaded_bulk =
read_csv("Input_data/FoodBalanceSheets_E_All_Data/FoodBalanceSheets_E_All_Data_NOFLAG.csv") %>% #bulk downloaded from FAOSTAT on 21.2.22 with all options
pivot_longer(cols = Y2010:Y2019, names_to = "year",values_to = "value") %>% #year to long format
mutate(year = sub('Y', '', year)) %>% #removing the Y from year values
clean_names()
# FAO mapping with FOOD groups (downloaded the list files without --------
FBS_FAO_map =
FBS_FAO_downloaded_EUgroup %>% as_tibble() %>% clean_names()%>%
dplyr::filter(year %in% c("2010":"2019"), #for a year average
# area_code %in% country_eu, #filtering for EU+UKD countries
grepl('Protein|Import|Export',element))%>%
dplyr::select(-c(domain_code,domain,area_code, element_code,item_code,year_code,unit,flag,flag_description)) %>%
left_join(dat_FAO_map, by="item") %>%
distinct(item, .keep_all = T) %>%
select(item,FoodGroups)
# FBS_FAO_map =
# FBS_FAO_downloaded_bulk %>% as_tibble() %>% clean_names()%>%
# dplyr::filter(year %in% c("2010":"2019"), #for a year average
# # area_code %in% country_eu, #filtering for EU+UKD countries
# grepl('Protein|Import|Export',element))%>%
# dplyr::select(-c(area_code, element_code,item_code,unit)) %>%
# left_join(dat_FAO_map, by="item") %>%
# distinct(item, .keep_all = T) %>%
# dplyr::select(item,FoodGroups)
# fao food group mapping --------------------------------------------------
# write_csv(FBS_FAO_map,"Input_data/FBS_FAO_prot_food_supply_mapin.csv")
fao_foodgroup_map = read_csv("Input_data/FBS_FAO_prot_food_supply_mapout.csv") #this is the final mapping tp ALL fao and eat lancet food groups
# Protein supply ----------------------------------------------------------
FBS_FAO_prot_food_supply =
FBS_FAO_downloaded_EUgroup %>% as_tibble() %>% clean_names()%>%
# FBS_FAO_prot_food_supply =
FBS_FAO_downloaded_bulk %>% as_tibble() %>% clean_names()%>%
dplyr::filter(year %in% c("2010":"2019"), #for a year average
area_code %in% country_eu, #filtering for EU+UKD countries
grepl('Protein|(kg/capita/yr)',element)) %>%
dplyr::select(-c(domain_code,domain,area_code, element_code,item_code,year_code,unit,flag,flag_description)) %>%
dplyr::select(-c(area_code, element_code,item_code,unit)) #%>%
pivot_wider(names_from = element,
values_from = value) %>%
values_from = value) #%>%
clean_names() %>%
replace(is.na(.), 0) %>%
left_join(fao_foodgroup_map, by="item") %>%
......@@ -2508,7 +2515,20 @@ write_csv(FBS_FAO_food_supply, "Input_data/FBS_FAO_food_supply.csv")
# Import_Export -----------------------------------------------------------
# Meeting with Renee on the 21.02.22
# Renee: that was done for validation data
# - 2014-2018 average
# - EU 28 (group not list)
# - all product that were in the new food balance sheet in the FAO (2013 - )
# - Stocks: Calculated the net import and net export.
# Part of import is exported and part of import is also kept as stock
# import - export - stock - non_food_uses(feed, other, seed, ... ) > negative value net_export; vice versa
# For whole EU not for countries
FBS_FAO_import_export =
FBS_FAO_downloaded_EUgroup %>% as_tibble() %>% clean_names() %>%
dplyr::filter(year %in% c("2010":"2019"), #for a year average
......@@ -2524,7 +2544,7 @@ FBS_FAO_import_export =
Export_ton = mean(export_quantity)) %>%
dplyr::mutate_if(is.numeric, ~ . * 1000)
write_csv(FBS_FAO_import_export, "Input_data/FBS_FAO_import_export.csv")
write_csv(FBS_FAO_import_export, "Input_data/FBS_FAO_import_export.csv")
# Testing the amount of proteins we get: 76.4 g protein per day with the weighted mean calculation
FBS_FAO_prot_food_supply %>% summarise(prot = sum(protein_gppp))
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment