From ad24682ae5bca8ee132a583e7e74942c5e9f69c6 Mon Sep 17 00:00:00 2001 From: "WUR\\simon083" <wolfram.simon@wur.nl> Date: Thu, 29 Jul 2021 22:37:22 +0200 Subject: [PATCH] New zonation based on ADM0 - AEZ - DOMSOIL_IPCC(mineral and Organic) - Validation with EUROSTAT and FAO. Data transformation. Different types of agricultural land extractions --- Zones/Zones_AEZ.R | 255 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 255 insertions(+) create mode 100644 Zones/Zones_AEZ.R diff --git a/Zones/Zones_AEZ.R b/Zones/Zones_AEZ.R new file mode 100644 index 0000000..4efedac --- /dev/null +++ b/Zones/Zones_AEZ.R @@ -0,0 +1,255 @@ +easypackages::packages(c("tidyverse", 'sf', 'raster','sp','dplyr','spatialEco', + 'exactextractr','tibble', 'here', "readr", "margrittr","tydr", + "readxl", "foreign", "data.table", "spData","spDataLarge", "purrr", "rgdal", "lwgeom" +)) + + +# Loading data ----------------------------------------------------------- +dat = st_read("Input_data/ADM_AEZ_IPCCDOMS.shp", agr= c(AEZ="identity", gridcode_o="identity", DOM_IPCCso="identity", ADM1_CODE="identity", ADM0_CODE="identity", Adm0_Adm1="identity")) + +# Loading GAUL shape files +shape_gaul=read_sf(here::here("Input_data", "g2015_2014_1.shp")) %>% st_set_geometry(NULL) %>% as_tibble() %>% + dplyr::select(ADM0_NAME, ADM0_CODE) + + +# Creating the final zones ------------------------------------------------------ +dat_zon = dat %>% dplyr::select(gridcode_o, DOM_IPCCso, ADM1_CODE, ADM0_CODE, Adm0_Adm1, geometry) %>% + rename(DOMS = DOM_IPCCso, + AEZ = gridcode_o) + +# Checking different levels of the IPCC soil variables -------------------- + + +dat_sol = dat %>% dplyr::select(DOM_IPCCso,SEQ1_IPCCs, DOM_Number, DOM_IPCC_1,ADM1_CODE, ADM0_CODE, Adm0_Adm1) %>% st_set_geometry(NULL) %>% as_tibble() + +# ADm1 level +adm1 = dat_sol %>% dplyr::select(Adm0_Adm1) %>% distinct() +# 3288 + +# ADM0 level +dist = dat_sol %>% dplyr::select(ADM0_CODE) %>% distinct() +# 256 + +dist = dat_sol %>% dplyr::select(DOM_IPCCso) %>% distinct() + +# DOM_IPCCso +# 1 HAC - Low activity clay soils +# 2 WR - Water +# 3 POD - Spodic soils +# 4 ORG - Organic Soils (Histosols) +# 5 LAC - Low activity clay soils +# 6 SAN - Sandy soils +# 7 ND - No data +# 8 WET - Wetlands +# 9 GG - Glaciers / Land ice +# 10 RK - Rock outcrops +# 11 VOL - Volcanic soils +# 12 ST - Salt flats +# 13 <NA> - NA + +dist2 = dat_sol %>% dplyr::select(SEQ1_IPCCs) %>% distinct() +# > dist2 +# SEQ1_IPCCs +# 1 HAC +# 2 WR +# 3 POD +# 4 ORG +# 5 LAC +# 6 san +# 7 ND +# 8 WET +# 9 GG +# 10 RK +# 11 SAN +# 12 VOL +# 13 HAc +# 14 ST +# 15 <NA> + +dist3 = dat_sol %>% dplyr::select(DOM_Number) %>% distinct() +# DOM_Number +# 1 1 +# 2 3 +# 3 2 +# 4 7 +# 5 4 +# 6 5 +# 7 6 +# 8 8 +# 9 9 +# 10 0 + +dist4 = dat_sol %>% dplyr::select(DOM_IPCC_1) %>% distinct() +# DOM_IPCC_1 +# 1 HAC +# 2 WR +# 3 POD +# 4 ORG +# 5 LAC +# 6 san +# 7 ND +# 8 WET +# 9 GG +# 10 RK +# 11 SAN +# 12 VOL +# 13 HAc +# 14 ST +# 15 <NA> + + +# Creating Moisture factors based on AEZ-33 ----------------------------- +# Translation from AEZ54 to AEZ33 done in excel: C:\Users\simon083\OneDrive - Wageningen University & Research\PhD_WJS\Academic\RQ1\Text_docs\Method_desc\Documentation\Documentation_paper1.xlsx +dat_moist = dat_zon %>% + mutate(Moist = case_when(AEZ == "2"~ "wet", + AEZ == "3"~ "wet", + AEZ == "5"~ "wet", + AEZ =="6"~ "wet", + AEZ =="8"~ "wet", + AEZ =="9"~ "wet", + AEZ =="11"~ "wet", + AEZ =="12"~ "wet", + AEZ =="14"~ "wet", + AEZ =="15"~ "wet", + AEZ =="17"~ "wet", + AEZ =="18"~ "wet", + AEZ =="20"~ "wet", + AEZ =="21"~ "wet", + AEZ =="23"~ "wet", + AEZ =="24"~ "wet", + AEZ =="30" ~ "wet"), + Moist = case_when(AEZ =="1" ~ "dry", + AEZ =="4" ~ "dry", + AEZ =="7" ~ "dry", + AEZ =="10" ~ "dry", + AEZ =="13" ~ "dry", + AEZ =="16" ~ "dry", + AEZ =="19" ~ "dry", + AEZ =="22" ~ "dry", + AEZ =="29" ~ "dry",TRUE ~ Moist), + Moist = case_when(AEZ =="25" ~ "other", + AEZ =="26" ~ "other", + AEZ =="28" ~ "other", + AEZ =="31" ~ "other", + AEZ =="32" ~ "other",TRUE ~ Moist), + Moist = case_when(AEZ =="33" ~ "water",TRUE ~ Moist), + Moist = case_when(AEZ =="27" ~ "irrigated",TRUE ~ Moist)) + + +# Creating Soil factors based on Dominant Soil Types from IPCC2009. +dat_soil = dat_moist %>% + mutate(Soil = case_when(DOMS =="HAC"~ "mineral", + DOMS =="LAC"~ "mineral", + DOMS =="SAN"~ "mineral", + DOMS =="POD"~ "mineral", + DOMS =="VOL"~ "mineral", + DOMS =="WET"~ "mineral", + DOMS =="WR"~ "other", + DOMS =="ND"~ "other", + DOMS =="GG"~ "other", + DOMS =="RK"~ "other"), + Soil = case_when(DOMS =="ORG"~ "organic",TRUE ~ Soil)) %>% + dplyr::select(AEZ,DOMS, ADM1_CODE, ADM0_CODE, Adm0_Adm1, Moist, Soil,geometry) + + +# 1 HAC - Low activity clay soils +# 2 WR - Water +# 3 POD - Spodic soils +# 4 ORG - Organic Soils (Histosols) +# 5 LAC - Low activity clay soils +# 6 SAN - Sandy soils +# 7 ND - No data +# 8 WET - Wetlands +# 9 GG - Glaciers / Land ice +# 10 RK - Rock outcrops +# 11 VOL - Volcanic soils +# 12 ST - Salt flats +# 13 <NA> - NA + +# Write and read shp files to set the identity for all vars. Otherwise, I cannot use the summarize function below. +write_sf(dat_soil, here("Input_data", "dat_soil_moist.shp")) +dat_soil_moist = st_read("Input_data/dat_soil_moist.shp", + agr= c(AEZ="identity",DOMS="identity", ADM1_CODE="identity", + Adm0_Adm1="identity", ADM0_CODE="identity", + Moist="identity", Soil="identity")) +# st_agr(dat_soil_moist) # Check if identity for all vars (--> allows aggregation below using summarize()) + + +# Create different zonation data sets ---------------------------------------- +# ADM0 - SOIL - AEZ +# Sys.time() +Zone_ADM0_DOMS_AEZ = dat_soil_moist %>% group_by(ADM0_CODE, DOMS, AEZ) %>% + summarise() # 6,317 more rows +write_sf(Zone_ADM0_DOMS_AEZ, here("Input_data", "Zone_ADM0_DOMS_AEZ.shp")) + +# ADM1 - SOIL - AEZ +Zone_ADM1_DOMS_AEZ = dat_soil_moist %>% group_by(Adm0_Adm1, DOMS, AEZ) %>% + summarise() +write_sf(Zone_ADM1_DOMS_AEZ, here("Input_data", "Zone_ADM1_DOMS_AEZ.shp")) + +# ADM0 - SOILsimp - AEZ +Zone_ADM0_SOIL_AEZ = dat_soil_moist %>% group_by(ADM0_CODE, Soil, AEZ) %>% + summarise() +write_sf(Zone_ADM0_SOIL_AEZ, here("Input_data", "Zone_ADM0_SOIL_AEZ.shp")) + +# ADM1 - SOILsimp - AEZ +Zone_ADM1_SOIL_AEZ = dat_soil_moist %>% group_by(Adm0_Adm1, Soil, AEZ) %>% + summarise() +write_sf(Zone_ADM1_SOIL_AEZ, here("Input_data", "Zone_ADM1_SOIL_AEZ.shp")) + +# ADM0 - SOILsimp - MOIST +Zone_ADM0_SOIL_Moist = dat_soil_moist %>% group_by(ADM0_CODE, Soil, Moist) %>% + summarise() +write_sf(Zone_ADM0_SOIL_Moist, here("Input_data", "Zone_ADM0_SOIL_Moist.shp")) + +# ADM1 - SOILsimp - MOIST +Zone_ADM1_SOIL_Moist = dat_soil_moist %>% group_by(Adm0_Adm1, Soil, Moist) %>% + summarise() +write_sf(Zone_ADM1_SOIL_Moist, here("Input_data", "Zone_ADM1_SOIL_Moist.shp")) +Sys.time() + +# Extraction of leandros agri land files (10km) based on the ADM0-AEZ-SOI -------- +# Loading the data +zon_cropland = st_read("Input_data/Agricultural_land/cropland.shp")#, agr= c(AEZ="identity", gridcode_o="identity", DOM_IPCCso="identity", ADM1_CODE="identity", ADM0_CODE="identity", Adm0_Adm1="identity")) +zon_pasture = st_read("Input_data/Agricultural_land/pasture.shp")#, agr= c(AEZ="identity", gridcode_o="identity", DOM_IPCCso="identity", ADM1_CODE="identity", ADM0_CODE="identity", Adm0_Adm1="identity")) +zon_pasturecrop = st_read("Input_data/Agricultural_land/pasturecrop.shp")#, agr= c(AEZ="identity", gridcode_o="identity", DOM_IPCCso="identity", ADM1_CODE="identity", ADM0_CODE="identity", Adm0_Adm1="identity")) +zon_rangeland = st_read("Input_data/Agricultural_land/rangeland.shp")#, agr= c(AEZ="identity", gridcode_o="identity", DOM_IPCCso="identity", ADM1_CODE="identity", ADM0_CODE="identity", Adm0_Adm1="identity")) + +# Combining agricultural land layers ------------------------------------- +Comb_Zon = cbind(zon_cropland, zon_pasture, zon_pasturecrop,zon_rangeland) %>% + dplyr::select(ADM0_CODE,Soil,AEZ,X_croplands, X_pasturesu, X_pasturecr, X_rangeland) %>% + rename(crop_areaha = X_croplands, + pasture_areaha= X_pasturesu, + pastureCrop_areaha = X_pasturecr, + rangeland_areaha = X_rangeland) %>% + mutate(across(c("crop_areaha", "pasture_areaha", "pastureCrop_areaha", "rangeland_areaha"), ~ (.*100))) + + +# Checking the area per country +# adm0 = shape_gaul %>% dplyr::select(ADM0_CODE, ADM0_NAME) %>% distinct_at(vars(ADM0_CODE, ADM0_NAME)) +# +# Chk_country = Comb_Zon %>% st_set_geometry(NULL) %>% as_tibble() %>% +# group_by(ADM0_CODE) %>% +# summarise(across(c("crop_areaha", "pasture_areaha", "pastureCrop_areaha", "rangeland_areaha"), ~ sum(., na.rm = TRUE))) %>% +# left_join(.,adm0,by="ADM0_CODE") +# +# Germany_cropland = Chk_country %>% filter(ADM0_NAME == "Germany") +# Netherlands_cropland = Chk_country %>% filter(ADM0_NAME == "Netherlands") +# Switzerland_cropland = Chk_country %>% filter(ADM0_NAME == "Switzerland") +# France_cropland = Chk_country %>% filter(ADM0_NAME == "France") +# Brazil_cropland = Chk_country %>% filter(ADM0_NAME == "Brazil") +# China_cropland = Chk_country %>% filter(ADM0_NAME == "China") +# India_cropland = Chk_country %>% filter(ADM0_NAME == "India") +# Mongolia_cropland = Chk_country %>% filter(ADM0_NAME == "Mongolia") +# Australia_cropland = Chk_country %>% filter(ADM0_NAME == "Australia") +# NewZealand_cropland = Chk_country %>% filter(ADM0_NAME == "New Zealand") +# USA_cropland = Chk_country %>% filter(ADM0_NAME == "United States of America") +# Russia_cropland = Chk_country %>% filter(ADM0_NAME == "Russia") + + + + + + + + -- GitLab