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

Restored the lost grassland script. All is ready for scaling now.

parent 259bee71
......@@ -11,7 +11,6 @@ country_map = read.csv(here::here("Input_data", "Mappings",
"dat_map_SPAM_GAUL_CIFOS_FAO.csv")) %>%
as_tibble()%>% dplyr::select(ADM0_CODE_GAUL, iso3_SPAM, name_cntr_CIFOS2021)
# Loading the dataset ----------------------------------------------------
# Th`is dataset was created by Leandro based on all the different LPJML raster datasets
dat_NPP = read_csv("Input_data/NPP_grassland/output/zone_adm0_soilagri_aez_analysis.csv" )
......@@ -19,7 +18,7 @@ dat_NPP = read_csv("Input_data/NPP_grassland/output/zone_adm0_soilagri_aez_analy
# Quality - Pasture and rangelands (Table from Pablo Modernel's grassland paper)
# arable grassland
pasture_arable_quality_EU = read_excel("Input_data/Pasture_rangeland_quality.xlsx", # Ryegrass + Fescue + Dactylis + White clover
sheet ="pasture_quality") %>%
sheet ="pasture_quality") %>%
janitor::clean_names() %>%
dplyr::filter(continent =="Europe and Russia") %>%
dplyr::select(-climate) %>% dplyr::select(-c(1)) %>%
......@@ -29,10 +28,10 @@ pasture_arable_quality_EU = read_excel("Input_data/Pasture_rangeland_quality.xls
# marginal grassland
pasture_marginal_quality_EU = pasture_arable_quality_EU %>%
dplyr::mutate(lut_cifos = dplyr::recode(lut_cifos, "pasture_arable" = "pasture_marginal"))
# Rangeland
rangeland_quality_EU =read_excel("Input_data/Pasture_rangeland_quality.xlsx",
sheet ="rangeland_quality") %>%
sheet ="rangeland_quality") %>%
janitor::clean_names() %>%
dplyr::filter(continent1=="Europe") %>%
dplyr::select(-ash_percent) %>%
......@@ -45,10 +44,10 @@ rangeland_quality_EU =read_excel("Input_data/Pasture_rangeland_quality.xlsx",
grass_quality = bind_rows(pasture_arable_quality_EU,
pasture_marginal_quality_EU,
rangeland_quality_EU) #%>%
# tidyr::pivot_longer(cols = dry_matter_percent:metabolizable_energy_mj_kg_dm,
# names_to = "content")
# tidyr::pivot_longer(cols = dry_matter_percent:metabolizable_energy_mj_kg_dm,
# names_to = "content")
write_csv(grass_quality, "Input_data/grass_quality.csv")
# Creating scale factors -------------------------------------------------
# With Pablo I discussed that there is several factors that lead to differences of uptake
# Stocking rate is positively related to the amount of grass that is taken up.
......@@ -92,68 +91,36 @@ dat_grass_crop =
# Calculating the above and below ground biomass
# Note: on the 12/11/2021 I talked to Tom Schut (PPS). He told me that LPJML refers to ABG and BG Biomass for NPP. 1 m ROOTING ZONE IS ASSUMED.
dplyr::mutate(yield_DM_kgha = (value/0.45*10), #Conversion factor from https://edepot.wur.nl/156130 -> Grassland simulation with the LPJmL model Version 3.4.018
# , #
# Grass_ABG_DM_kgha = Tot_grass_DM_kgha*(1-(1/(1+mean(c(1.54,2.54))))), # Shoot:root ratio grassland -> https://www.sciencedirect.com/science/article/pii/S0378429017301132?via%3Dihub
# Grass_BG_DM_kgha = Tot_grass_DM_kgha-Grass_ABG_DM_kgha,
# #converting DM to FM by dividing the DM/frac_DM of grass
# , #
# Grass_ABG_DM_kgha = Tot_grass_DM_kgha*(1-(1/(1+mean(c(1.54,2.54))))), # Shoot:root ratio grassland -> https://www.sciencedirect.com/science/article/pii/S0378429017301132?via%3Dihub
# Grass_BG_DM_kgha = Tot_grass_DM_kgha-Grass_ABG_DM_kgha,
# #converting DM to FM by dividing the DM/frac_DM of grass
yield_FM_kgha = yield_DM_kgha/(1-mean(c(0.70,0.86)))) %>% # Assumption: LWC (0.70–0.86 g g−1) -- Source: DOI: https://doi.org/10.1111/j.1469-8137.1994.tb04036.x; Leaf anatomy, specific mass and water content in congeneric annual and perennial grass species
dplyr::select(-first,-three ,-second,-five, -six, -seven, -eight, -value, -...1) %>%
dplyr::rename(stock_rate = four)
# # adding the ratios per lut and stocking rates
# # dplyr::mutate(lut_fct = case_when(lut_cifos == "pasture_arable" ~ 0.75,
# # lut_cifos == "pasture_marginal" ~ 0.5,
# # lut_cifos == "rangeland" ~ 0.25)) %>%
# # Adding the scales for stocking intensity
# dplyr::left_join(.,stock_tbl, by="stock_rate") %>%
# # Calculating the adjusted Yields
# dplyr::mutate(Grass_ABG_DM_kgha_adj = Grass_ABG_DM_kgha* stock_rate_fct, #*lut_fct
# FM_ABG_grass_kgha_adj = FM_ABG_grass_kgha * stock_rate_fct)
# 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)
# # adding the ratios per lut and stocking rates
# # dplyr::mutate(lut_fct = case_when(lut_cifos == "pasture_arable" ~ 0.75,
# # lut_cifos == "pasture_marginal" ~ 0.5,
# # lut_cifos == "rangeland" ~ 0.25)) %>%
# # Adding the scales for stocking intensity
# dplyr::left_join(.,stock_tbl, by="stock_rate") %>%
# # Calculating the adjusted Yields
# dplyr::mutate(Grass_ABG_DM_kgha_adj = Grass_ABG_DM_kgha* stock_rate_fct, #*lut_fct
# FM_ABG_grass_kgha_adj = FM_ABG_grass_kgha * stock_rate_fct)
# 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") %>%
dplyr::mutate(area_ha = area_ha *100) %>%
dplyr::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"))
# 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")
......@@ -164,11 +131,11 @@ lut_area_av =
available_land_area %>% filter(., iso3_SPAM %in% EU_iso3) %>%
dplyr::select(-c(yield_DM_kgha, yield_FM_kgha, iso3_SPAM, adm0_cod )) %>%
relocate(name_cntr_CIFOS2021, .before = soil) %>%
relocate(soil, .before = lut_cifos) %>%
rename(country_cifos = name_cntr_CIFOS2021,
soil_cifos = soil,
aez_code =aez)
relocate(soil, .before = lut_cifos) %>%
rename(country_cifos = name_cntr_CIFOS2021,
soil_cifos = soil,
aez_code =aez)
write_csv(lut_area_av, "Input_data/land_av_lut.csv")
......@@ -182,8 +149,8 @@ lut_area_av_absolute =
soil_cifos = soil,
aez_code =aez) %>%
filter(lut_cifos != "pasture_arable")
write_csv(lut_area_av_absolute, "Input_data/land_av_lut_absolute.csv")
write_csv(lut_area_av_absolute, "Input_data/land_av_lut_absolute.csv")
# Pasture_arable only
......@@ -201,7 +168,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
......@@ -212,8 +180,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)
......@@ -221,54 +188,320 @@ grass_EU28_final =
grass_EU28 %>%
# group_by(soil,aez, stock_rate, lut_cifos) %>%
# dplyr::summarise(yield_FM_kgha = weighted.mean(x=yield_FM_kgha,
# w = area_ha, #weighted mean
# yield_FM_kgha[!is.na(yield_FM_kgha) & yield_FM_kgha!= 0], na.rm = TRUE),
# area_ha=sum(area_ha)) %>%
# w = area_ha, #weighted mean
# yield_FM_kgha[!is.na(yield_FM_kgha) & yield_FM_kgha!= 0], na.rm = TRUE),
# area_ha=sum(area_ha)) %>%
dplyr::mutate(yield_FM_kgha=replace_na(yield_FM_kgha ,0)) %>%
# dplyr::select(-yield_FM_kgha) %>%
dplyr::rename(prod_level = stock_rate) %>%
dplyr::rename(prod_level = stock_rate) %>%
dplyr::mutate(cifos_crop = case_when(lut_cifos == "pasture_arable" ~ "grass_arable",
lut_cifos == "pasture_marginal" ~ "grass_marginal",
lut_cifos == "rangeland" ~ "grass_rangeland",
TRUE ~ "grass_other"),
aez = as.character(aez)) %>%
dplyr::select(-c(adm0_cod , yield_DM_kgha, iso3_SPAM)) %>%
dplyr::rename(soil_cifos = soil,
country_cifos = name_cntr_CIFOS2021,
crop_cifos = cifos_crop,
aez_code = aez,
yield_kgha = yield_FM_kgha) %>%
relocate(soil_cifos, .after = aez_code) %>%
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",
"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::select(-c(adm0_cod , yield_DM_kgha, iso3_SPAM)) %>%
dplyr::rename(soil_cifos = soil,
country_cifos = name_cntr_CIFOS2021,
crop_cifos = cifos_crop,
aez_code = aez,
yield_kgha = yield_FM_kgha) %>%
relocate(soil_cifos, .after = aez_code) %>%
relocate(crop_cifos, .after = lut_cifos) %>%
relocate(country_cifos, .before = aez_code) %>%
relocate(lut_cifos, .before = prod_level) %>%
relocate(yield_kgha, .before = area_ha)
# Adding stocking rate current and potential ------------------------------
# Subsetting per stocking rate --------------------------------------------
grassing_density_potential_init =
dat_grass_crop %>%
dplyr::mutate(area_ha = area_ha *100) %>%
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) %>%
dplyr::filter(prod_level == "10") %>%
mutate(prod_level = "A") %>%
#scenario = "current") %>%
# relocate(scenario, .before = crop_cifos) %>%
ungroup() %>%
# dplyr::select(-adm0_cod) %>%
# 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 %>%
dplyr::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)) %>%
dplyr::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_potential = grass_EU28_final %>% ungroup() %>% right_join(grassing_density_potential_cntr %>% ungroup() %>%
dplyr::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")) %>%
filter(!crop_cifos == "grass_other") %>%
mutate(prod_level=paste0("pot_",prod_level)) %>%
dplyr::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_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", "pasture_arable"))
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
gg_leandro_eurostat_arable=Fun_ggplot_area_deviation(arable_area_cifos, lut_vec="arable", author_vec = "Leandro", vec_source = "EUROSTAT")
gg_leandro_eurostat_grass=Fun_ggplot_area_deviation(grass_area_cifos,lut_vec="grass", author_vec = "Leandro" , vec_source = "EUROSTAT")
# Comparing the own extraction from R script 'Zones_AEZ.R' -----------------
Comb_Zon = read_csv(here::here("Input_data", "AEZ_ADM0_SOIL_agri_global.csv")) #global zonation with land type spatial area
# SUbsetting zonation variables with land type area for EU28
Comb_Zon_EU28 = Comb_Zon %>% filter(., iso3_SPAM %in% EU_iso3)
Fun_area_join_WJSextract = function(dat,vec_lut){
Comb_Zon_EU28 %>% dplyr::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) %>%
dplyr::group_by(country_cifos) %>%
dplyr::summarise(area_ha_cifos = sum(area_ha,na.rm=T)) %>%
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
gg_WJS_eurostat_arable=Fun_ggplot_area_deviation(grass_area_cifos_WJex, lut_vec="arable", author_vec = "WJS",vec_source = "EUROSTAT")
gg_WJS_eurostat_grass=Fun_ggplot_area_deviation(arable_area_cifos_WJex,lut_vec="grass", author_vec = "WJS", vec_source = "EUROSTAT")
# comparing with FAO area ---------------------------------------------------------------------
EU_cntr_Cifos = c("Austria","Belgium","Bulgaria","Croatia","Cyprus", "Czechia","Denmark",
"Estonia","Finland","France","Germany","Greece","Hungary","Ireland",
"Italy","Latvia","Lithuania","Luxembourg","Malta","Netherlands","Poland",
"Portugal","Romania","Slovakia","Slovenia","Spain","Sweden",'United Kingdom')
landuse_fao_bulk = read_csv("INput_data/Inputs_LandUse_E_All_Data_NOFLAG.csv")
# For checking
check_dat_LUT = landuse_fao_bulk %>% select(Area, Item,Y2008,Y2009,Y2010,Y2011,Y2012) %>%
filter(Item %in% c("Arable land", "Land under permanent crops",
"Land under perm. meadows and pastures")) %>%
clean_names() %>% type.convert() %>%
mutate(area = recode(area, "United Kingdom of Great Britain and Northern Ireland"= "United Kingdom")) %>%
filter(., area %in% EU_cntr_Cifos) %>%
pivot_longer(cols = "y2008":"y2012",
names_to = "year",
values_to = "area_ha") %>%
dplyr::mutate(area_ha = area_ha*1000) %>%
rename("country" = "area") %>%
group_by(country, item) %>%
summarise(area_ha = mean(area_ha,na.rm=T))
fao_arable_grass_area =
landuse_fao_bulk %>%
clean_names() %>% type.convert() %>%
mutate(area = recode(area, "United Kingdom of Great Britain and Northern Ireland"= "United Kingdom")) %>%
filter(., area %in% EU_cntr_Cifos) %>%
pivot_longer(cols = "y1961":"y2019",
names_to = "year",
values_to = "area_ha") %>%
dplyr::mutate(area_ha = area_ha*1000) %>%
# HYDE paper: We started with country totals for the FAO categories of
# “arable land and permanent crops” and “permanent mead-
# ows and pastures”, further referred to here as “cropland”
# and “grazing land” respectively - https://doi.org/10.5194/essd-9-927-2017
dplyr::filter(year %in% c("y2009","y2010","y2011"),
item %in% c("Arable land", "Land under permanent crops",
"Land under perm. meadows and pastures")) %>%
rename("country" = "area") %>%
group_by(country, item) %>%
summarise(area_ha = mean(area_ha,na.rm=T))
# Summing over items
# arable
arable_area_fao = fao_arable_grass_area %>% dplyr::filter(item %in% c("Arable land","Land under permanent crops")) %>%
group_by(country) %>%
summarise(area_ha = sum(area_ha,na.rm=T))
# grass
grass_area_fao = fao_arable_grass_area %>%
dplyr::filter(item %in% c("Land under perm. meadows and pastures")) %>%
group_by(country) %>%
summarise(area_ha = sum(area_ha,na.rm=T))
# Applying the function with FAO arable area data to the extracted datasets of Leandro and Wolfram
arable_area_cifos_WJex_fao = Fun_area_join_WJSextract(dat = arable_area_fao, vec_lut = c("arable"))
arable_area_cifos_fao = Fun_area_join(dat = grass_area_fao, vec_lut = c("arable"))
grass_area_cifos_WJex_fao = Fun_area_join_WJSextract(dat = grass_area_fao, vec_lut = c("pasture_marginal", "rangeland","pasture_arable"))
grass_area_cifos_fao = Fun_area_join(dat = grass_area_fao, vec_lut = c("pasture_marginal", "rangeland","pasture_arable"))
# Plotting the two extraction types from Leandro and Wolfram
gg_WJS_fao_arable = Fun_ggplot_area_deviation(arable_area_cifos_WJex_fao, lut_vec="arable", author_vec = "WJS",vec_source = "FAOSTAT - 3 year average around 2010")
gg_leandro_fao_arable = Fun_ggplot_area_deviation(arable_area_cifos_fao, lut_vec="arable", author_vec = "Leandro",vec_source = "FAOSTAT - 3 year average around 2010")
gg_WJS_fao_grass = Fun_ggplot_area_deviation(grass_area_cifos_WJex_fao, lut_vec="grass", author_vec = "WJS",vec_source = "FAOSTAT - 3 year average around 2010")
gg_leandro_fao_grass = Fun_ggplot_area_deviation(grass_area_cifos_fao, lut_vec="grass", author_vec = "Leandro",vec_source = "FAOSTAT - 3 year average around 2010")
# Conclusion: Leandros extraction is way to low (average around -60% of the FAOSTAT reference)
# Arable land
ggarrange(gg_leandro_eurostat_arable, gg_leandro_fao_arable,
gg_WJS_eurostat_arable, gg_WJS_fao_arable)
# Grass area
ggarrange(gg_leandro_eurostat_grass, gg_WJS_eurostat_grass, nrow = 2)
# Grass area
ggarrange(gg_WJS_fao_grass, gg_leandro_fao_grass, nrow = 2)
# 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
......@@ -278,8 +511,8 @@ weighted_yld_agg =
dplyr::ungroup() %>%
dplyr::group_by(aez_code, soil_cifos, crop_cifos) %>%
dplyr::summarise(yield_kgha = weighted.mean(x=yield_FM_kgha, #aggregating the stocking rate
w = area_ha, #weighted mean
yield_FM_kgha [!is.na(yield_FM_kgha ) & yield_FM_kgha != 0])) %>% #adding uncertanties SE) #%>% #exluding 0's.
w = area_ha, #weighted mean
yield_FM_kgha [!is.na(yield_FM_kgha ) & yield_FM_kgha != 0])) %>% #adding uncertanties SE) #%>% #exluding 0's.
dplyr::mutate(yield_kgha=replace_na(yield_kgha ,0))
grassEU_stock_avg =
......@@ -310,16 +543,16 @@ grassEU_stock_avg_mean=
group_by(lut_cifos) %>%
summarise(yield_kgha_mean = mean(yield_kgha),
yield_kgha_wmean = weighted.mean(x=yield_kgha ,
w = area_ha, #weighted mean
yield_kgha [!is.na(yield_kgha ) & yield_kgha != 0]))
w = area_ha, #weighted mean
yield_kgha [!is.na(yield_kgha ) & yield_kgha != 0]))
# Cretating an html table
kbl(grassEU_stock_avg_mean) %>% kable_classic(full_width = F, html_font = "Cambria")
# column_spec(6, color = "white",
# background = spec_color(grassEU_stock_avg_mean$yield_kgha_wmean[1:18], end = 0.7))
# popover = paste("am:", grassEU_stock_avg_mean$yield_kgha_wmean[1:18]))
# column_spec(6, color = "white",
# background = spec_color(grassEU_stock_avg_mean$yield_kgha_wmean[1:18], end = 0.7))
# popover = paste("am:", grassEU_stock_avg_mean$yield_kgha_wmean[1:18]))
# kable_styling("hover", full_width = F)
# kable_styling("hover", full_width = F)
# Clustering stock rates -------------------------------------------------
# K-Means Cluster Analysis
......@@ -334,7 +567,7 @@ kbl(grassEU_stock_avg_mean) %>% kable_classic(full_width = F, html_font = "Camb
# Visualization ----------------------------------------------------------
# Comparing the Stocking rates per country - weigthed mean - adjusted--------------------------------
# DM_yields_lines_dot_per_cnt_meanadj =
grass_EU28 %>%
grass_EU28 %>%
group_by(iso3_SPAM,stock_rate) %>%
# dplyr::summarise(# Grass_ABG_DM_kgha_mean = median(Grass_ABG_DM_kgha_adj, na.rm = TRUE)) %>%
dplyr::summarise(
......@@ -343,14 +576,14 @@ kbl(grassEU_stock_avg_mean) %>% kable_classic(full_width = F, html_font = "Camb
w = area_ha, #weighted mean
yield_DM_kgha[!is.na(yield_DM_kgha) &
yield_DM_kgha != 0])) %>%
ggplot(aes(x=stock_rate, y=yield_DM_kgha , group=iso3_SPAM, col=iso3_SPAM)) +
ggplot(aes(x=stock_rate, y=yield_DM_kgha , group=iso3_SPAM, col=iso3_SPAM)) +
geom_line()+
geom_point()+
facet_wrap(~iso3_SPAM)
# Fresh weight
# FM_yields_lines_dot_per_cnt_meanadj =