AWESOME.Rmd

---
title: "GlobalErosionDB"
author: "Jinshi"
date: "2/2/2020"
output: html_document
---

# Load package
```{r load package, message=FALSE, include=FALSE, echo=FALSE}
# https://stackoverflow.com/questions/4090169/elegant-way-to-check-for-missing-packages-and-install-them/19873732

package_list <- c("cowplot","data.table","dplyr","ggplot2", "lubridate", "leaflet"
                  ,"kableExtra","knitr","ggmap","maps","mapdata","tidyr","sp","ggpubr"
                  ,"readxl")
package_new <- package_list[!(package_list %in% installed.packages()[,"Package"])]
if(length(package_new)) install.packages(package_new)

library(cowplot)
library(data.table)
library(dplyr)
library(ggplot2)
theme_set(theme_bw())
library(lubridate)
library(kableExtra)
library(knitr)
library("ggpubr")
library(tidyr)
library("leaflet")
library(sp)
library(readxl)
# install.packages("ggridges")
library(ggridges)
library(patchwork)
require(foreign)
require(MASS)
# devtools::install_github('kongdd/Ipaper')
library(Ipaper)
source("functions.R")
```

# Setup
```{r preliminaries, message=FALSE, include=FALSE, echo=FALSE, cache=TRUE}
# Set chunks defaults, these options will be applied to all subsequent chunks
knitr::opts_chunk$set(message = TRUE, include = TRUE, echo = FALSE,
                      fig.height = 4, fig.width = 8)

```

# load data
```{r load data}
wos_summary <- read_excel('data/Number_Studies_byYear.xlsx', sheet = 2, skip = 1)
colnames(wos_summary)
SEDB_del <- drake::readd(SEDB_del)
counties <- drake::readd(counties)
GlobalMATMAP <- drake::readd(GlobalMATMAP)
IGBP_MODIS <- drake::readd(IGBP_MODIS)

## join and get climate classification, vegetation information
SEDB_del %>% mutate(Latitude2 = round(Latitude*2)/2,
                 Longitude2 = round(Longitude*2)/2,
                 Lat_dif = ifelse(Latitude2 - Latitude >=0, 0.25, -0.25),
                 Lon_dif = ifelse(Longitude2 - Longitude >=0, 0.25, -0.25),
                 Latitude2 = Latitude2 - Lat_dif,
                 Longitude2 = Longitude2 - Lon_dif) %>% 
  dplyr::select(-Lat_dif, -Lon_dif) -> 
  SEDB_del

# Get Ecosystem class, MAT and MAP for srdb data
left_join(SEDB_del, IGBP_MODIS, by=c("Latitude2"="Latitude", "Longitude2"="Longitude")) ->
  SEDB_del

# Get MAT and MAP
left_join(SEDB_del, GlobalMATMAP, by=c("Latitude2"="Latitude", "Longitude2"="Longitude")) ->
  SEDB_del

SEDB_del %>% 
  dplyr::select(Study_midyear, Latitude, Longitude, MAT, MAP, Study_temp, Study_precip, MAT_Del, MAP_Del, Tannual_del, Pannual_del) 
# Change Latitude and Longitude to 0.5 resolution for IGBP

```

```{r SEDB processing}
# vegetation group
SEDB_del %>% dplyr::select(IGBP.x) %>% unique()

SEDB_del %>%
  mutate(IGBP = case_when(IGBP.x %in% c("OSH", "CSH") ~ "CSH",
                          IGBP.x %in% c("EBF", "EDF", "ENF", "MF", "BDF", "FOR", "DBF", "BEF") ~ "FOR",
                          IGBP.x %in% c("CRO") ~ "CRO",
                          IGBP.x %in% c("GRA", "SAV", "PST") ~ "GRA",
                          IGBP.x %in% c("BSV") ~ "BSV",
                          IGBP.x %in% c("PLT") ~ "PLT",
                          IGBP.x %in% c("ORC", "OCD") ~ "OCD",
                          IGBP.x %in% c("MIX") ~ "MIX",
                          TRUE ~ "OTH")) -> SEDB_del

SEDB_del %>% dplyr::select(IGBP.x, IGBP) %>% filter(IGBP == "OTH" & IGBP.x !="")

# soil group
SEDB_del %>% 
  dplyr::select(Soil_texture, ER_annual) %>% 
  filter(!is.na(ER_annual)) %>% 
  group_by(Soil_texture) %>% 
  summarise(ER_mean = mean(ER_annual),
            obs = n(),
            se = sd(ER_annual) / sqrt(obs))

SEDB_del %>%
  mutate(Soil_Group = case_when(Soil_texture %in% c("clay", "Clay") ~ "C",
                          Soil_texture %in% c("Sandy clay") ~ "SC",
                          Soil_texture %in% c("clay loam", "Clay loam", "Clay Loam") ~ "CL",
                          Soil_texture %in% c("Loamy clay", "Silty clay") ~ "SiC",
                          Soil_texture %in% c("Silty clay loam", "Silty Clay Loam") ~ "SiCL",
                          Soil_texture %in% c("Sandy clay loam", "Sandy Clay Loam") ~ "SCL",
                          Soil_texture %in% c("loam", "Loam") ~ "L",
                          Soil_texture %in% c("Silt loam", "Silt Loam") ~ "SiL",
                          Soil_texture %in% c("Sand") ~ "S",
                          Soil_texture %in% c("sandy loam", "Sandy loam", "Sandy Loam") ~ "SL",
                          Soil_texture %in% c("Loamy sand", "Loamy Sand") ~ "LS",
                          Soil_texture %in% c("Silt") ~ "Si",
                          TRUE ~ "OTH")) -> SEDB_del

SEDB_del %>% dplyr::select(Soil_texture, Soil_Group) %>% filter(Soil_Group == "OTH" & Soil_texture !="")


# Measure_method group
SEDB_del %>% 
  dplyr::select(Meas_method, ER_annual) %>% 
  filter(!is.na(ER_annual)) %>% 
  group_by(Meas_method) %>% 
  summarise(ER_mean = mean(ER_annual),
            obs = n(),
            se = sd(ER_annual) / sqrt(obs))

SEDB_del %>%
  mutate(Meas_method_group = case_when(Meas_method %in% c("Runoff plot") ~ "Plot",
    Meas_method %in% c("Erosion pins", "Topography", "Limnigraph", "Microprofiler", "Profilometer") ~ "Pins",
    Meas_method %in% c("Caesium 137", "Be7", "Pb210", "Radioisotopic", "Tracer") ~ "Isotopic",
    Meas_method %in% c("Catchment") ~ "Catchment",
    Meas_method %in% c("Bathymetry", "Gauging station") ~ "Bathymetry",
    Meas_method %in% c("Field measurement") ~ "Field measurement",
    Meas_method %in% c("Field rainfall simulation") ~ "Field rainfall simulation",
    Meas_method %in% c("AGNPS model", "ASEAT", "ASWAT", "ATEM/SEDEM model", "Model", "MUSLE", "RMMF", "RUSLE", "RUSLE2",
                       "SWAT", "USLE", "WATEM/SEDEM model", "WEPP") ~ "Model",
    TRUE ~ "OTH"  )) -> SEDB_del

SEDB_del %>% dplyr::select(Meas_method, Meas_method_group) %>% filter(Meas_method_group == "OTH" & Meas_method !="")

SEDB_del %>% dplyr::select(Study_number, Meas_method_group) %>% filter(Meas_method_group == "Model") %>% count(Study_number)
SEDB_del %>% dplyr::select(Study_number, Meas_method_group) %>% filter(Meas_method_group == "Isotopic") %>% count(Study_number)

```


```{r output SEDB}
writ_file(SEDB_del, "SoilErosion_out.csv")
```

# plot number of study by year
```{r number of study by year, fig.width=8, fig.height=9}
bind_rows(
  wos_summary %>% 
    mutate(Source = "(a) Number of papers published between 1960 and 2021"),
  
  SEDB_del %>% 
    filter(!is.na(Paper_year) & Study_number != 259) %>% 
    dplyr::select(Paper_year) %>% 
    count(Paper_year) %>% 
    mutate(Year = floor(Paper_year), Num_study = n, Source = "(b) Number of measurements by publication year in SoilErosionDB"),
  
  SEDB_del %>% 
    filter(!is.na(Study_midyear) & Study_number != 259) %>% 
    dplyr::select(Study_midyear) %>% 
    count(Study_midyear) %>% 
    mutate(Year = floor(Study_midyear), Num_study = n, Source = "(c) Number of measurements by measure year in SoilErosionDB")
  ) %>% 
  ggplot(aes(x=Year, y=Num_study)) +
  geom_bar(stat = "sum") +
  facet_wrap(Source ~ ., scales = "free", ncol = 1) +
  labs(x=expression(Year), y=expression(Count~(n))) +
  scale_x_continuous(breaks = seq(1900,2020,10)) +
  theme(legend.position = "none")
ggsave("outputs/Figure 1 wos summary.jpg", width = 8, height = 7)
```

```{r prepare data for site map}
SEDB_del %>% 
  dplyr::select(Paper_year, Study_midyear) %>% 
  na.omit() %>% 
  ggplot(aes(Paper_year, Study_midyear)) +
  geom_point()

SEDB_del %>% 
  dplyr::select(Latitude, Longitude) %>% 
  na.omit() %>% 
  unique()

# site_plot using leaflet function
bind_rows(
  SEDB_del %>% 
    # mutate(ER_annual = coalesce(ER_annual, ER_g_m2_h)) %>% 
    dplyr::select(Longitude, Latitude, ER_annual) %>% 
    filter(!is.na(Longitude) & !is.na(ER_annual)) %>%
    group_by(Longitude, Latitude) %>% 
    summarise(count = n()) %>% 
    mutate(Indicator = "Erosion"),
  
  SEDB_del %>% 
    dplyr::select(Longitude, Latitude, Sediment_annual) %>% 
    filter(!is.na(Longitude) & !is.na(Sediment_annual)) %>%
    group_by(Longitude, Latitude) %>% 
    summarise(count = n()) %>% 
    mutate(Indicator = "Sediment"),
  
  SEDB_del %>% dplyr::select(Longitude, Latitude, Runoff_annual) %>% 
    filter(!is.na(Longitude) & !is.na(Runoff_annual)) %>%
    group_by(Longitude, Latitude) %>% 
    summarise(count = n()) %>% 
    mutate(Indicator = "Runoff"),
  
  SEDB_del %>% dplyr::select(Longitude, Latitude, Leaching) %>% 
    filter(!is.na(Longitude) & !is.na(Leaching)) %>%
    group_by(Longitude, Latitude) %>% 
    summarise(count = n()) %>% 
    mutate(Indicator = "Leaching") ) ->  SEDB_site

# SEDB_site %>% dplyr::select(Longitude, Latitude) %>% 
#   na.omit() %>% 
#   unique() %>% leaflet() %>% 
#   addMarkers() %>% 
#   addTiles(options = providerTileOptions(minZoom = 1, maxZoom = 3)) 
```

# site map
```{r plot site map, fig.height=4, fig.width=8}
# global map
# sort(unique(counties$region))
ggplot(data = counties) + 
  geom_polygon(aes(x = long, y = lat, group = group),
               color = "white", fill = 'gray', alpha = 0.25) + 
  guides(fill=FALSE) +
  geom_point(data = SEDB_site,
             aes(x=Longitude, y=Latitude,
                 size = count, 
                 col = Indicator,
                 shape = Indicator),
             alpha = 0.5,
             stroke = 1 ) +
  scale_shape_manual(values = c(1, 16, 3, 83)) +
  scale_color_manual(values = c("black", "blue", "red", "skyblue")) +
  scale_x_continuous(name="Longitude", breaks=seq(-180,180, 60),labels = seq(-180,180, 60))+
  scale_y_continuous(limits = c(-60, 80),name="Latitude", breaks=seq(-60,80,30),labels = seq(-60,80,30)) +
  scale_size_continuous(name = "Obs (n)") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# country_lat_long_check(SEDB_del %>% filter(Country %in% c("China", "USA")))

ggsave("outputs/Figure 2. Sites distribution.png", width = 8, height = 4, dpi = 300, units = "in" )
```

# Scatter plot
```{r plot erosion vs runoff}
SEDB_del %>% dplyr::select(ER_annual, Runoff_annual, IGBP.y) %>% na.omit() %>% 
  ggplot(aes(x = Runoff_annual, y = ER_annual, col = IGBP.y)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = expression(Annual~runoff~"(mm)"),
      y = expression(Annual~erosion~"(t ha"^{-1}~"yr"^{-1}~")"))

SEDB_del %>% dplyr::select(Sediment_annual, Runoff_annual, IGBP.y) %>% na.omit() %>% 
  ggplot(aes(x = Runoff_annual, y = Sediment_annual, col = IGBP.y)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = expression(Runoff~"rate (mm hr"^{-1}~")"),
      y = expression(Annual~erosion~"(g m"^{-2}~"hr"^{-1}~")"))

SEDB_del %>% dplyr::select(ER_annual, Rainfall_amount, IGBP.y) %>% na.omit() %>% 
  ggplot(aes(x = Rainfall_amount, y = ER_annual, col = IGBP.y)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = expression(Annual~rainfall~"(mm)"),
      y = expression(Annual~erosion~"(t ha"^{-1}~"yr"^{-1}~")"))

```

# density plot
```{r plot density, fig.height=4, fig.width=8}
plot_grid(
  SEDB_del %>% dplyr::select(Runoff_annual) %>% 
    na.omit() %>% 
    ggplot(aes(Runoff_annual)) +
    geom_histogram(fill = "gray", col = "black") +
    labs(x = expression(Annual~Runoff~"(mm yr"^{-1}~")"),
         y = expression("Density")),
  
  SEDB_del %>% dplyr::select(ER_annual) %>% 
    na.omit() %>% 
    filter(ER_annual <= 100) %>% 
    ggplot(aes(ER_annual)) +
    geom_histogram(fill = "gray", col = "black") +
    labs(x = expression(Annual~erosion~"(t ha"^{-1}~"yr"^{-1}~")"),
         y = expression("Density")),
  ncol = 2,
  hjust = -4, vjust = 3,
  labels = c("(a)", "(b)")
)

ggplot(iris, aes(x = Sepal.Length, y = Species, fill = 0.5 - abs(0.5 - stat(ecdf)))) +
  stat_density_ridges(geom = "density_ridges_gradient", calc_ecdf = TRUE) +
  scale_fill_viridis_c(name = "Tail probability", direction = -1)

SEDB_del %>% dplyr::select(Runoff_annual, Ecosystem2) %>% 
  dplyr::filter(Runoff_annual > 0) %>% 
  na.omit() %>% 
  ggplot(aes(Runoff_annual, Ecosystem2, fill = 0.5 - abs(0.5 - stat(ecdf)))) +
  stat_density_ridges(geom = "density_ridges_gradient", calc_ecdf = TRUE) +
  scale_fill_viridis_c(name = "Tail probability", direction = -1) +
  labs(x = expression(Annual~Runoff~"(mm yr"^{-1}~")"),
       y = expression("Density")) +
  # scale_y_discrete(breaks = c("Wetland", "T&B FOR", "Shrubland", "Savanna", "Grassland", "Agriculture"),
  #                  labels = c("Wetland", "Forest", "Shrubland", "Savanna", "Grassland", "Agriculture")) +
  theme(legend.position = "none") ->
  p1

SEDB_del %>% dplyr::select(ER_annual, Ecosystem2) %>% 
  na.omit() %>% 
  filter(ER_annual <= 50 & Ecosystem2 != "Urban" & ER_annual > 0) %>%
  ggplot(aes(ER_annual, Ecosystem2, fill = 0.5 - abs(0.5 - stat(ecdf)))) +
  stat_density_ridges(geom = "density_ridges_gradient", calc_ecdf = TRUE) +
  scale_fill_viridis_c(name = "Probability", direction = -1) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "right") +
  labs(x = expression(Annual~erosion~"(t ha"^{-1}~"yr"^{-1}~")"),
       y = expression("Density")) ->
  p2

p1 + p2

# use boxplot
SEDB_del %>% 
  filter(!is.na(ER_annual)) %>% 
  ggplot(aes(x=IGBP.x, y=ER_annual)) +
  # geom_hline(yintercept = 0, linetype = 2, col = "red") + 
  # geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
  geom_boxplot2(lwd = 0.5, width = 0.5, alpha=1, fatten = NULL,
               fill = "orange", col = "orange") +
  stat_summary(fun.y=median, geom="point", shape=16, size=4.5, col = "gray") +
  stat_summary(fun.y=median, geom="point", shape=16, size=1.5, col = "black") +
  labs(y = expression(Annual~soil~erosion~(Mg/ha/yr))) 
  # facet_grid(cols = vars(data)) 

```



```{r Figure 3, plot erosion rate by different group}
# by vegetation type
SEDB_del %>% 
  dplyr::select(IGBP, ER_annual) %>% 
  filter(!is.na(ER_annual)) %>% 
  group_by(IGBP) %>% 
  summarise(ER_mean = mean(ER_annual),
            obs = n(),
            se = sd(ER_annual) / sqrt(obs)) %>% 
  filter(IGBP != "OTH") -> SEDB_IGBP

SEDB_IGBP %>% 
  ggplot(aes(IGBP, ER_mean)) +
  geom_bar(stat = "identity", col = "black", fill = "white", width = 0.75) +
  geom_errorbar(aes(ymin=ER_mean - se, ymax = ER_mean + se), col = "black", width = 0.05) +
  coord_cartesian(ylim = c(0, 100)) +
  annotate("text", x = c(1:8), y = 80,
           label = paste0("n=",SEDB_IGBP$obs), angle = 90, hjust = 0) +
  labs(x=expression(Vegetation~type),
       y=expression(Soil~erosion~(t~ha^{-1}~yr^{-1}))) -> plot_igbp
  
# by soil texture
SEDB_del %>% 
  dplyr::select(Soil_Group, ER_annual) %>% 
  filter(!is.na(ER_annual)) %>% 
  group_by(Soil_Group) %>% 
  summarise(ER_mean = mean(ER_annual),
            obs = n(),
            se = sd(ER_annual) / sqrt(obs)) %>% 
  filter(Soil_Group != "SC") %>% 
  filter(Soil_Group != "OTH") -> SEDB_soil

SEDB_soil %>% 
  ggplot(aes(Soil_Group, ER_mean)) +
  geom_bar(stat = "identity", col = "black", fill = "white", width = 0.75) +
  geom_errorbar(aes(ymin=ER_mean - se, ymax = ER_mean + se), col = "black", width = 0.05) +
  coord_cartesian(ylim = c(0, 100)) +
  annotate("text", x = c(1:11), y = 80,
           label = paste0("n=",SEDB_soil$obs), angle = 90, hjust = 0) +
  labs(x=expression(Soil~type),
       y=expression(Soil~erosion~(t~ha^{-1}~yr^{-1}))) -> plot_soil

# measure method
SEDB_del %>% 
  dplyr::select(Meas_method_group, ER_annual) %>% 
  filter(!is.na(ER_annual)) %>% 
  group_by(Meas_method_group) %>% 
  summarise(ER_mean = mean(ER_annual),
            obs = n(),
            se = sd(ER_annual) / sqrt(obs)) %>% 
  filter(Meas_method_group != "Field rainfall simulation") %>% 
  filter(Meas_method_group != "OTH") -> SEDB_meas


SEDB_meas %>% 
  ggplot(aes(Meas_method_group, ER_mean)) +
  geom_bar(stat = "identity", col = "black", fill = "white", width = 0.75) +
  geom_errorbar(aes(ymin=ER_mean - se, ymax = ER_mean + se), col = "black", width = 0.05) +
  coord_cartesian(ylim = c(0, 310)) +
  annotate("text", x = c(1:5), y = 250,
           label = paste0("n=", SEDB_meas$obs), angle = 90, hjust = 0) +
  labs(x=expression(Measure~method),
       y=expression(Soil~erosion~(t~ha^{-1}~yr^{-1}))) -> plot_meas



SEDB_scale %>% 
  ggplot(aes(Field_scale, ER_mean)) +
  geom_bar(stat = "identity", col = "black", fill = "white", width = 0.75) +
  geom_errorbar(aes(ymin=ER_mean - se, ymax = ER_mean + se), col = "black", width = 0.05) +
  coord_cartesian(ylim = c(0, 80)) +
  annotate("text", x = c(1:3), y = 62,
           label = paste0("n=", SEDB_scale$obs), angle = 90, hjust = 0) +
  labs(x=expression(Field~scale),
       y=expression(Soil~erosion~(t~ha^{-1}~yr^{-1}))) -> plot_scale

plot_igbp / plot_soil / plot_meas -> plot_catg
plot_catg + plot_annotation(tag_levels = 'a', tag_prefix = '(',  tag_suffix = ')')

ggsave("outputs/Figure 3. Sistribution of erosion and runoff.png", width = 6, height = 10, dpi = 300, units = "in" )
```

## Make tables for the ESSD manuscript
```{r Table 1}
# Table 1
SEDB_meas
knitr::kable(SEDB_meas)
```


```{r Table 4}
# Table 4
SEDB_del %>% 
  dplyr::select(Quality_flag) %>% 
  count(Quality_flag)
```

```{r Table 5}
# Table 5
SEDB_del %>% 
  dplyr::select(Manipulation) %>% 
  count(Manipulation)
```




```{r Table S1}
# Table X
SEDB_del %>% 
  dplyr::select(Leaching_type) %>% 
  count(Leaching_type)
```


```{r Figure 6}
# Link with external data and plot
# SEDB_del$Study_precip[SEDB_del$Study_number == 82] <- 1632
SEDB_del %>% 
  dplyr::select(Study_midyear, Latitude, Longitude, MAT, MAP, Study_temp, Study_precip, MAT_Del, MAP_Del, Tannual_del, Pannual_del) 

bind_rows(
  SEDB_del %>% 
    dplyr::select(MAT, MAT_Del) %>% 
    rename(Study = MAT, Del = MAT_Del) %>% 
    na.omit() %>% 
    mutate(Type = "(a) MAT"),
  
  SEDB_del %>% 
    dplyr::select(Study_temp, Tannual_del) %>% 
    rename(Study = Study_temp, Del = Tannual_del) %>% 
    na.omit() %>% 
    mutate(Type = "(b) Annual temperature"),
  
  SEDB_del %>% 
    dplyr::select(MAP, MAP_Del) %>% 
    filter(MAP < 6000) %>% 
    rename(Study = MAP, Del = MAP_Del) %>% 
    na.omit() %>% 
    mutate(Type = "(c) MAP (mm)"),
  
  SEDB_del %>% 
    dplyr::select(Study_precip, Pannual_del) %>% 
    rename(Study = Study_precip, Del = Pannual_del) %>% 
    na.omit() %>% 
    mutate(Type = "(d) Annual precipitation (mm)")) %>% 
  ggplot(aes(Del, Study)) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_abline(col = "red", linetype = "dotdash", size = 1.15) +
  facet_wrap(~ Type, ncol = 2, scales = "free") +
  labs(x = expression(Values~from~the~global~climate~dataset),
       y = expression(Values~from~papers))

# ggsave("outputs/Figure 6. Tm Pm relationship.png", width = 8, height = 6, dpi = 300, units = "in" )

SEDB_del %>% 
  dplyr::select(Study_temp, Tannual_del) %>% 
  rename(Study = Study_temp, Del = Tannual_del) %>% 
  na.omit() %>% 
  ggplot(aes(Del, Study)) +
  geom_point() +
  geom_smooth(method = "lm")
    
```

## Soil erosion vs variables relationship

```{r Table 6}
# linear regression results between soil erosion and numeric factors
lm_er_var <- function(sdata, var_column) {
  
  sdata %>%
    dplyr::select(ER_annual, {{var_column}}) %>% 
    filter(ER_annual > 0) %>% na.omit() -> sub_data
  
  var_column = colnames (sub_data)[2]
  
  colnames (sub_data) = c("ER_annual", "var")
  
  sub_data %>% 
    group_by(var) %>% 
    summarise(Obs = n(), ER_annual = mean(ER_annual)) %>% 
    mutate(ER_annual = log(ER_annual)) -> sdata_agg
  
  # linear regression
  slm <- lm(sdata_agg$ER_annual ~ sdata_agg$var, weights = sdata_agg$Obs) 
  summary(slm) -> sum_lm
  
  outputs = tibble(var_num = var_column,
                   intercept_a <- summary(slm)$coefficients[1,1] %>% round(6),
                   slope_b <- summary(slm)$coefficients[2,1] %>% round(6),
                   p_slope_b <- summary(slm)$coefficients[2,4]%>% round(6),
                   R2 = sum_lm$r.squared,
                   n = nrow(sdata_agg) )
  colnames(outputs) <- c("Var", "intercept", "slope", "p_slope", "R2", "n")
  return(outputs)
}

bind_rows(
  lm_er_var(SEDB_del, Pannual_del),
  lm_er_var(SEDB_del, Tannual_del),
  lm_er_var(SEDB_del, MPET),
  lm_er_var(SEDB_del, Soil_sand),
  lm_er_var(SEDB_del, Soil_silt),
  lm_er_var(SEDB_del, Soil_clay),
  lm_er_var(SEDB_del, Soil_BD),
  lm_er_var(SEDB_del, Soil_pH),
  lm_er_var(SEDB_del, Soil_SOC),
  lm_er_var(SEDB_del, Soil_N),
  lm_er_var(SEDB_del, Soil_CN),
  lm_er_var(SEDB_del, Ksat),
  lm_er_var(SEDB_del, L),
  lm_er_var(SEDB_del, Slope),
  lm_er_var(SEDB_del, Field_area),
  lm_er_var(SEDB_del, LAI),
  lm_er_var(SEDB_del, Plant_cover_C)
)


# Robust regression
# rust_lm <- rlm(SEDB_agg$ER_annual ~ SEDB_agg$Pannual_del, weights = SEDB_agg$Obs)
# summary(rust_lm)
# cooks.distance( rust_lm ) %>% max()

```



```{r plot ER_annual vs num variables, fig.height=10, fig.width=8}
# function for testing ER_annual vs numeric relationship
plot_er_var <- function(sdata, var_column) {
  sdata %>%
    dplyr::select(ER_annual, {{var_column}}) %>% 
    filter(ER_annual > 0) %>% na.omit() -> sdata_agg
  
  colnames (sdata_agg) = c("ER_annual", "var")
  
  sdata_agg %>% 
    group_by(var) %>% 
    summarise(Obs = n(), ER_annual = mean(ER_annual)) %>% 
    mutate(ER_annual = log(ER_annual)) -> sdata_agg
  
  # scatter plot
  sdata_agg %>% 
    ggplot(aes(x = var, y = ER_annual, size = Obs)) +
    geom_point(alpha = 0.75) +
    geom_smooth(method = "lm", show.legend = FALSE) +
    scale_size_continuous(name = "Obs (n)") -> plot_var
  
  return(plot_var)
}

## a. ER_annual vs Precipitation
plot_er_var(SEDB_del, Pannual_del) +
  labs(x = expression(Annual~precipitation~"(mm)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  annotate("text", x = 2700, y = -8, label = "italic(R) ^ 2 == 0.030", hjust = 0, parse = TRUE) +
  ggtitle('(a)') -> plot_precp


# b. vs clay
plot_er_var(SEDB_del, Soil_clay) +
  labs(x = expression(Soil~clay~"(%)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  annotate("text", x = 55, y = -4, label = "italic(R) ^ 2 == 0.119", hjust = 0, parse = TRUE) +
  ggtitle('(b)') -> plot_clay


# c vs bulk density
plot_er_var(SEDB_del, Soil_BD) +
  labs(x = expression(Soil~BD~"("~g~m^3~")"),
       y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  # theme(axis.title.y = element_blank()) +
  annotate("text", x = 1.5, y = -3, label = "italic(R) ^ 2 == 0.180", hjust = 0, parse = TRUE) +
  ggtitle('(c)') -> plot_bd

# d vs pH
plot_er_var(SEDB_del, Soil_pH) +
  labs(x = expression(Soil~pH),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  theme(axis.title.y = element_blank()) +
  annotate("text", x = 7.5, y = -6, label = "italic(R) ^ 2 == 0.108", hjust = 0, parse = TRUE) +
  ggtitle('(d)') -> plot_pH

# e vs SOC
plot_er_var(SEDB_del, Soil_SOC) +
  labs(x = expression(Soil~SOC~"(%)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  theme(axis.title.y = element_blank()) +
  annotate("text", x = 8, y = -6, label = "italic(R) ^ 2 == 0.066", hjust = 0, parse = TRUE) +
  ggtitle('(e)') -> plot_soc

# f. vs plant coverage
plot_er_var(SEDB_del, Plant_cover_C) +
  labs(x = expression(Plant~coverage~"(%)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  theme(axis.title.y = element_blank()) +
  annotate("text", x = 15, y = -5, label = "italic(R) ^ 2 == 0.240", hjust = 0, parse = TRUE) +
  ggtitle('(f)') -> plot_plant


(plot_precp  / plot_clay / plot_bd) | (plot_pH / plot_soc / plot_plant) 
# ggsave("outputs/Figure 4. Erosion v relationship.png", width = 8, height = 8, dpi = 300, units = "in" )
```

```{r plot var not significant}
# function for testing ER_annual vs numeric relationship
plot_er_var2 <- function(sdata, var_column) {
  sdata %>%
    dplyr::select(ER_annual, {{var_column}}) %>% 
    filter(ER_annual > 0) %>% na.omit() -> sdata_agg
  
  colnames (sdata_agg) = c("ER_annual", "var")
  
  sdata_agg %>% 
    group_by(var) %>% 
    summarise(Obs = n(), ER_annual = mean(ER_annual)) %>% 
    mutate(ER_annual = log(ER_annual)) -> sdata_agg
  
  # scatter plot
  sdata_agg %>% 
    ggplot(aes(x = var, y = ER_annual, size = Obs)) +
    geom_point(alpha = 0.75) +
    # geom_smooth(method = "lm", show.legend = FALSE) +
    scale_size_continuous(name = "Obs (n)") -> plot_var
  
  return(plot_var)
}

## a. ER_annual vs MPET
plot_er_var2(SEDB_del, MPET) +
  labs(x = expression(Annual~evaporation~"(mm)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(a)') -> plot_mpet

## b. ER_annual vs sand
plot_er_var2(SEDB_del, Soil_sand) +
  labs(x = expression(Soil~sand~"(%)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(b)') -> plot_sand

## c. ER_annual vs N
plot_er_var2(SEDB_del, Soil_N) +
  labs(x = expression(Soil~N~"(%)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(c)') -> plot_n

## d. ER_annual vs cn
plot_er_var2(SEDB_del %>% filter(Soil_CN < 50), Soil_CN) +
  labs(x = expression(Soil~"C/N"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(d)') -> plot_cn

## e. ER_annual vs ksat
plot_er_var2(SEDB_del %>% filter(Ksat < 50), Ksat) +
  labs(x = expression(Ksat~"(cm/h)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(e)') -> plot_ksat

  
## f. ER_annual vs L
plot_er_var2(SEDB_del %>% filter(L < 20000), L) +
  labs(x = expression(Slope~length~"(m)"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(f)') -> plot_l


# g. vs annual temperature
plot_er_var2(SEDB_del, Tannual_del) +
  labs(x = expression(Annual~temperature~(degree~C)),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  ggtitle('(g)') -> plot_temp

## h vs Field area
SEDB_del %>% 
  dplyr::select(ER_annual, Field_area) %>%
  na.omit() %>% 
  filter(ER_annual > 0 & Field_area > 0) %>% 
  group_by(Field_area) %>% 
  summarise(Obs = n(), ER_annual = mean(ER_annual)) %>% 
  mutate(ER_annual = log(ER_annual), Field_area = log(Field_area)) ->
  SEDB_agg_fieldarea

lm(SEDB_agg_fieldarea$ER_annual ~ SEDB_agg_fieldarea$Field_area, weights = SEDB_agg_fieldarea$Obs) %>% summary()

SEDB_agg_fieldarea %>% 
  ggplot(aes(x = Field_area, y = ER_annual, size = Obs)) +
  geom_point(alpha = 0.75) +
  labs(x = expression(Field~'area, '~log~"("~m^{-2}~")" ),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  theme(axis.title.y = element_blank()) +
  scale_size_continuous(name = "Obs (n)") +
  ggtitle('(h)') -> plot_area

```



```{r}
(plot_mpet / plot_sand / plot_n / plot_cn) | (plot_ksat / plot_l / plot_temp / plot_area) 
ggsave("outputs/Figure SX. Erosion var.png", width = 9, height = 10, dpi = 300, units = "in" )
```



```{r test slope}
plot_er_var(SEDB_del %>% filter(Field_scale == "Runoff plot"), Slope) +
  labs(x = expression(Slope~"("~degree~")"),
      y = expression('Erosion,'~log~"(t ha"^{-1}~"yr"^{-1}~")")) +
  annotate("text", x = 15, y = -10, label = "italic(R) ^ 2 == 0.002", hjust = 0, parse = TRUE)

SEDB_del %>% 
  mutate(IGBP = case_when(IGBP.x %in% c("BDF", "BEF", "DBF", "EBF", "EDF", "ENF", "FOR", "MF", "PLT") ~ "FOR",
                          IGBP.x %in% "BSV" ~ "BSV", 
                          IGBP.x %in% "CRO" ~ "CRO",
                          IGBP.x %in% c("CSH", "OSH", "SAV") ~ "CSH",
                          IGBP.x %in% c("GRA", "PST") ~ "GRA",
                          IGBP.x %in% c("OCD", "ORC") ~ "OCD",
                          # IGBP.x %in% c("", "") ~ "",
                          TRUE ~ "OTH")) -> SEDB_del

SEDB_del %>% filter(IGBP == "OTH") %>% 
  dplyr::select(IGBP.x, IGBP) %>% 
  count(IGBP.x)

SEDB_del %>%
  dplyr::select(ER_annual, Meas_method_group, Slope, IGBP.x) %>% 
  filter(ER_annual > 0) %>% na.omit() %>%  
  filter(Meas_method_group == "Runoff plot") %>% 
  ggplot(aes(x = Slope, y = log(ER_annual))) +
  geom_point(alpha = 0.75) +
  facet_wrap(.~IGBP.x) +
  geom_smooth(method = "lm", show.legend = FALSE) 

SEDB_del %>%
  dplyr::select(ER_annual, Meas_method_group, Slope, IGBP) %>% 
  filter(ER_annual > 0) %>% na.omit() %>%  
  filter(Meas_method_group == "Runoff plot") -> SEDB_veg

SEDB_veg %>% 
  dplyr::select(ER_annual, IGBP, Slope) %>% 
  group_by(IGBP, Slope) %>% 
  summarise(ER_annual = log(mean(ER_annual)),
            obs = n()) -> SEDB_veg

# BSV
SEDB_veg %>% 
  filter(IGBP %in% c("BSV")) ->
  sub1
lm(ER_annual ~ Slope, data = sub1) %>% summary()

# CRO
SEDB_veg %>% 
  filter(IGBP %in% c("CRO")) ->
  sub2
lm(ER_annual ~ Slope, data = sub2) %>% summary()

# CSH
SEDB_veg %>% 
  filter(IGBP %in% c("CSH")) ->
  sub3
lm(ER_annual ~ Slope, data = sub3) %>% summary()

# FOR
SEDB_veg %>% 
  filter(IGBP %in% c("FOR")) ->
  sub4
lm(ER_annual ~ Slope, data = sub4) %>% summary()

# GRA
SEDB_veg %>% 
  filter(IGBP %in% c("GRA")) ->
  sub5
lm(ER_annual ~ Slope, data = sub5) %>% summary()

# OCD
SEDB_veg %>% 
  filter(IGBP %in% c("OCD")) ->
  sub6
lm(ER_annual ~ Slope, data = sub6) %>% summary()

# OTH
SEDB_veg %>% 
  filter(IGBP %in% c("OTH")) ->
  sub7
lm(ER_annual ~ Slope, data = sub7) %>% summary()

```

```{r Er vs slope by vegetation, fig.height=6, fig.width=8}
SEDB_veg %>% 
  ggplot(aes(Slope, ER_annual)) +
  geom_point(aes(size  = obs), shape = 16, alpha = 0.5) +
  geom_smooth(method = "lm", data = SEDB_veg %>% filter(IGBP %in% c("CRO")), fill = "skyblue") +
  facet_wrap(.~ IGBP, scales = "free", nrow = 2) +
  # facet_grid(rows = vars(IGBP), scales = "free") +
  labs(x=expression("Slope ("~degree~")"), y = expression(Annual~soil~erosion~(t~ha^{-1}~yr^{-1}))) +
  labs(size="obs (n)")

# ggsave('outputs/Figure 5-v2.TestSlope.png', width = 8, height = 5)
```


```{r}
## ER_annual vs Elevation
SEDB_del %>% 
  dplyr::select(ER_annual, Elevation) %>%
  na.omit() %>% 
  filter(ER_annual > 0) %>% 
  group_by(Elevation) %>% 
  summarise(Obs = n(), ER_annual = mean(ER_annual)) %>% 
  mutate(ER_annual = log(ER_annual)) ->
  SEDB_agg_elev

lm(SEDB_agg_elev$ER_annual ~ SEDB_agg_elev$Elevation, weights = SEDB_agg_elev$Obs) %>% summary()
```



```{r}
GarciaRuiz_SI <- read.csv('References/259-GarciaRuiz/1-s2.0-S0169555X1500149X-mmc2.csv', sep = ";") %>% 
  filter(!is.na(ref))

write.csv(GarciaRuiz_SI, "References/259-GarciaRuiz/GarciaRuiz_SI.csv", row.names = FALSE)
```