New zonation based on ADM0 - AEZ - DOMSOIL_IPCC(mineral and Organic) -...

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

+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")
+
+
+
+
+
+
+
+