SPAM_GSS_ADM1_ZONES.R

easypackages::packages(c("tidyverse", 'sf', 'raster','sp','dplyr','spatialEco',
                         'exactextractr','tibble', 'here', "readr",
                         "readxl", "foreign", "data.table", "rbenchmark", "benchmarkme",
                         "validate", "plyr"))

# Set working directory
setwd(here::here("Input_data"))

# Load data ----------------------------------------------------------------
# Yields ------------------------------------------------------------------
# dat_TA_yld <- read.csv("spam2010V2r0_global_Y_TA.csv", locale = locale(encoding = "LATIN2"))# Yields - all technologies together, ie complete crop - Data downloaded on 1. of September 2020
# dat_TI_yld <- read_csv("spam2010V2r0_global_Y_TI.csv", locale = locale(encoding = "LATIN2"))# Yields - Iirrigated portion of crop
# dat_TH_yld <- read_csv("spam2010V2r0_global_Y_TH.csv", locale = locale(encoding = "LATIN2"))# Yields - rainfed high inputs portion of crop
# dat_TL_yld <- read_csv("spam2010V2r0_global_Y_TL.csv", locale = locale(encoding = "LATIN2"))# Yields - rainfed low inputs portion of crop
# dat_TS_yld <- read_csv("spam2010V2r0_global_Y_TS.csv",locale = locale(encoding = "LATIN2"))# Yields - rainfed subsistence portion of crop
# dat_TR_yld <- read_csv("spam2010V2r0_global_Y_TR.csv", locale = locale(encoding = "LATIN2"))# Yields - rainfed portion of crop (= TA - TI, or TH + TL + TS)

# Production
dat_TA_pr <- read_csv("spam2010V2r0_global_P_TA.csv",locale = locale(encoding = "LATIN2")) # LATIN2 encoding is used for proper display of nordic and slavic names
# dat_TI_pr <- read_csv("spam2010V2r0_global_P_TI.csv",locale = locale(encoding = "LATIN2"))
# dat_TH_pr <- read_csv("spam2010V2r0_global_P_TH.csv",locale = locale(encoding = "LATIN2"))
# dat_TL_pr <- read_csv("spam2010V2r0_global_P_TL.csv",locale = locale(encoding = "LATIN2"))
# dat_TS_pr <- read_csv("spam2010V2r0_global_P_TS.csv",locale = locale(encoding = "LATIN2"))
# dat_TR_pr <- read_csv("spam2010V2r0_global_P_TR.csv",locale = locale(encoding = "LATIN2"))
# 
# # Physical area--------------------------------
# dat_TA_pa <- read_csv("spam2010V2r0_global_A_TA.csv",locale = locale(encoding = "LATIN2"))
# dat_TI_pa <- read_csv("spam2010V2r0_global_A_TI.csv",locale = locale(encoding = "LATIN2"))
# dat_TH_pa <- read_csv("spam2010V2r0_global_A_TH.csv",locale = locale(encoding = "LATIN2"))
# dat_TL_pa <- read_csv("spam2010V2r0_global_A_TL.csv",locale = locale(encoding = "LATIN2"))
# dat_TS_pa <- read_csv("spam2010V2r0_global_A_TS.csv",locale = locale(encoding = "LATIN2"))
# dat_TR_pa <- read_csv("spam2010V2r0_global_A_TR.csv",locale = locale(encoding = "LATIN2"))
# 
# # Harvesting area----------------------------------
# dat_TA_ha <- read_csv("spam2010V2r0_global_H_TA.csv",locale = locale(encoding = "LATIN2"))
# dat_TI_ha <- read_csv("spam2010V2r0_global_H_TI.csv",locale = locale(encoding = "LATIN2"))
# dat_TH_ha <- read_csv("spam2010V2r0_global_H_TH.csv",locale = locale(encoding = "LATIN2"))
# dat_TL_ha <- read_csv("spam2010V2r0_global_H_TL.csv",locale = locale(encoding = "LATIN2"))
# dat_TS_ha <- read_csv("spam2010V2r0_global_H_TS.csv",locale = locale(encoding = "LATIN2"))
# dat_TR_ha <- read_csv("spam2010V2r0_global_H_TR.csv",locale = locale(encoding = "LATIN2"))

test <- read_csv("Test_distinct.txt")


# Create a spatial SPAM file 
# geo_spam=dplyr::select(dat_TA_pr, x,y, cell5m, name_adm1, iso3) 
# %>% 
# distinct_at(vars(name_adm1))

# Make SPAM spatial 
# projcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
# proj_mercator = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
proj = "+proj=longlat +datum=WGS84 +no_defs +units=m" # taken from raster file below 
proj_meter = "+datum=WGS84 +units=m +no_defs"
geo_spam1 <- st_as_sf(x = geo_spam,                         
                      coords = c("x", "y"),
                      crs = proj)

# Load administrative unit shapefile which matches the SPAM adm1 zones
# Admin0 = st_read(here::here("Input_data","GAUL_GADM_SPAM_ULRIKE", "g2008_0_now","g2008_0.shp"))
# Admin1 = st_read(here::here("Input_data","GAUL_GADM_SPAM_ULRIKE", "g2008_1_now","g2008_1.shp")) #this is what we need
# Admin2 = st_read(here::here("Input_data","GAUL_GADM_SPAM_ULRIKE", "g2008_2_2020_11_02","g2008_2.shp")) #this is what we need

# Here I manually adjusted the NAME1 values to the once used in SPAM2010. see Documentation for more. 
Admin1Adj = st_read(here::here("Input_data","g2008_1_now", "g2008_1.shp"), options = "ENCODING=LATIN2")#WINDOWS-1252

# Now the values from GGS that are still not matching the SPAM2010 adm1 values are recoded
Admin1Adj2 = Admin1Adj %>% as_tibble() %>% 
  mutate(NAME1 = case_when(ADM0_CODE == 196 & #Philipines
                             NAME1 =="Autonomous region in Muslim Mindanao (ARMM)"~ "Autonomous region in Muslim Mi",