rsofun_rsip_mct.Rmd

---
title: "rsofun at RSIP sites"
author: "Beni Stocker"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(rsofun)
library(ingestr)
library(rbeni)
library(cowplot)
library(ggridges)
```

## Get sites

```{r cars}
# df <- read_csv("~/data/rootingdepth/rsip/RSIP_Analysis_sheet_210409.csv") |>      # previously done with "data/RSIP_Analysis_sheet.csv"
df <- read_csv("~/data/rootingdepth/rsip/RSIP_Analysis_sheet_210721.csv") |>
  rename(lon = Long, lat = Lat) |> 
  rowid_to_column(var = "id") |> 
  
  ## problem: some have a reference error
  dplyr::filter(lon != "#REF!") |> 
  mutate(lon = as.numeric(lon), lat = as.numeric(lat), 
         Dr = as.numeric(Dr),
         wtd = as.numeric(Water_Table_Depth_Fan))
  
print(paste("Total number of entries:", nrow(df)))
```

### Aggregate

1705 distinct sites based on lon and lat info.
```{r}
df_sites <- df |> 
  mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))) |> 
  group_by(sitename) |> 
  summarise(Dr = mean(Dr, na.rm = TRUE), wtd = mean(wtd, na.rm = TRUE))

## add WWF biome info
df_wwf_sites <- ingest(
  df |> 
    distinct(lon, lat) |> 
    mutate(sitename = paste0("i_", as.character(lon), "_", as.character(lat))),
  source = "wwf",
  dir = "~/data/biomes/wwf_ecoregions/official/",
  settings = list(layer = "wwf_terr_ecos")
  )|> 
  mutate(data = purrr::map(data, ~slice(., 1))) |> 
  unnest(data)

df_sites <- df_sites |> 
  left_join(
    df_wwf_sites |> 
      dplyr::select(sitename, biome = BIOME, biome_name = BIOME_NAME),
    by = "sitename"
  )

df_sites <- df_sites %>% 
  separate(sitename, c("prefix", "lon", "lat"), sep = "_") %>% 
  dplyr::select(-prefix) %>% 
  mutate(sitename = paste0("i_", 1:nrow(.))) %>% 
  mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>%
  mutate(
    year_start = 1989,
    year_end = 2018
  )
```

### Add model results.
```{r}
df_sites <- extract_nc(df_sites |> dplyr::select(sitename, lon, lat), "data/zroot_cwd80.nc") |> 
  unnest(data) |> 
  rename(zroot_cwdx80 = V1) |> 
  mutate(zroot_cwdx80 = zroot_cwdx80 / 1000) |> 
  right_join(df_sites, by = c("lon", "lat", "sitename"))
```

<!-- ### Add elevation data -->

<!-- ```{r} -->
<!-- if(grepl('eu-', Sys.info()['nodename'])){ -->
<!--   r <- raster::raster("~/data/etopo/ETOPO1_Bed_g_geotiff.tif") -->
<!-- } -->

<!-- # grab elev data -->
<!-- s <- st_as_sf( -->
<!--   df, -->
<!--   coords = c("lon","lat"), -->
<!--   crs = 4326) -->
<!-- df$elv <- raster::extract(r, s) -->
<!-- df <- df |> -->
<!--   mutate( -->
<!--     elv = ifelse(elv < 0, 0, elv) -->
<!--   ) -->

### Write to file

```{r}
write_csv(df_sites, file = "data/df_sites_rsip.csv")
```


## Get forcing data

Use `submit_forcing_rsip.sh` to run `rscript_forcing_rsip.R`. Thid saves files in subdirectory `"./data/forcing_rsip/"`.

## Test

Use `submit_rsofun_cwdx.sh` to run `rscript_rsofun_cwdx.R` for different S_0 (WHC) and diagnose it following the S_dEF method. This runs `run_rsofun_cwdx_by_chunk.R` and saves files in subdirectory `"./data/out_rsofun_cwdx/"`. 

### Get data

```{r}
# ichunk <- 1
# use_whc <- 100

read_myfile <- function(ichunk, use_whc){
  path <- "./data/out_rsofun_cwdx/"
  filename <- file.path(path, paste0("out_rsofun_cwdx_whc_", as.character(use_whc), "_ichunk_", as.character(ichunk), ".rds"))
  if (file.exists(filename)){
    df <- readRDS(filename) |> 
      dplyr::select(sitename, whc, s_def, s_cwdx80)
    return(df)
  } else {
    return(tibble())
  }
}

df <- purrr::map_dfr(
  as.list(seq(50)),
  ~read_myfile(., 100)
  ) |> 
  mutate(setup = "whc_100") |> 
  bind_rows(
    purrr::map_dfr(
      as.list(seq(50)),
      ~read_myfile(., 200)
      ) |> 
      mutate(setup = "whc_200")
  ) |> 
  bind_rows(
    purrr::map_dfr(
      as.list(seq(50)),
      ~read_myfile(., "NA")
      ) |> 
      mutate(setup = "whc_NA")
  )

saveRDS(df, file = "data/df_test_cwdx_rsofun.rds")
```

### Scatterplots

```{r}
df <- readRDS("data/df_test_cwdx_rsofun.rds")

# out2 <- df |> 
#   dplyr::filter(setup == "whc_NA" & !is.nan(s_cwdx80) & !is.infinite(s_cwdx80) & s_cwdx80 < 1e4) |>
#   mutate(s_cwdx80 = remove_outliers(s_cwdx80)) |> 
#   analyse_modobs2("whc", "s_cwdx80", type = "hex")
# 
# out2$gg <- out2$gg + xlim(0, 2000) + ylim(0, 2000)

n_fun <- function(x){
  return(data.frame(y = 380,
                    label = length(x)))
}

out <- df |> 
  dplyr::filter(setup == "whc_NA") |> 
  analyse_modobs2("whc", "s_def", shortsubtitle = TRUE) 
gg1 <- out$gg +
  labs(x = expression("Prescribed" ~ "S"[0] ~ "(mm)"), y = expression("Diagnosed" ~ "S"[0] ~ "(mm)"))

gg2 <- df |> 
  dplyr::filter(setup != "whc_NA") |> 
  ggplot(aes(as.factor(whc), s_def)) +
  geom_violin(draw_quantiles = 0.5, fill = "grey70") +
  stat_summary(fun.data = n_fun, geom = "text", hjust = 0.5, size = 2.5) +
  labs(x = expression("Prescribed" ~ "S"[0] ~ "(mm)"), y = expression("Diagnosed" ~ "S"[0] ~ "(mm)")) +
  geom_hline(yintercept = c(100, 200), linetype = "dotted") +
  ylim(0, 380) +
  theme_classic()

toprow <- cowplot::plot_grid(gg2, gg1, ncol = 2, rel_widths = c(1, 1), labels = c('a', 'b'))
```

### By biome

```{r}
tmp <- df |> 
  left_join(df_sites, by = "sitename")

q75 <- function(x){quantile(probs = 0.75, ...)}

n_fun <- function(x){
  return(c(y = median(x) + 30, label = length(x))) 
}

## boxplot
gg3 <- tmp |> 
  dplyr::filter(!(biome_name %in% c(NA))) |>
  dplyr::filter(setup != "whc_NA") |> 
  ggplot(aes(x = biome_name, y = s_def, fill = as.factor(setup))) +
  geom_boxplot(outlier.shape = NA, position=position_dodge(0)) +
  stat_summary(fun.data = n_fun, geom = "text", fun.y = median,
               position = position_dodge(width = 0), size = 2.5) +
  scale_fill_manual(values = c("#E69F00", "#56B4E9"), labels = c("100", "200"), name = "Prescribed" ~ "S"[0] ~ "(mm)") +
  labs(x = expression("Biome"), y = expression("Diagnosed" ~ "S"[0] ~ "(mm)")) +
  geom_hline(yintercept = c(100, 200), linetype = "dotted") +
  ylim(0, 300) +
  theme_classic() +
  coord_flip()

cowplot::plot_grid(toprow, gg3, nrow = 2, rel_heights = c(1, 1.5), labels = c('', 'c'))
ggsave("fig/plot_rsofun_rsip.pdf", width = 10, height = 8)
ggsave("fig/plot_rsofun_rsip.png", width = 10, height = 8)
```


### ridge plot

```{r}
gg4 <- df |> 
  dplyr::filter(!(biome_name %in% biome_exclude)) |>
  filter(setup == "whc_NA") |> 
  rename(obs = whc, mod = s_def) |> 
  tidyr::pivot_longer(cols = c(mod, obs), names_to = "source", values_to = "s0") |> 
  dplyr::filter(!is.na(biome_name)) |> 
  dplyr::filter(biome_name!="Mangroves") |>   
  ggplot(aes(x = s0, y = biome_name, color = source, point_color = source, fill = source)) +
  geom_density_ridges(
    jittered_points = TRUE, scale = .95, rel_min_height = .01,
    point_shape = "|", point_size = 1.5, size = 0.25,
    position = position_points_jitter(height = 0)
  ) +
  scale_y_discrete(expand = c(0, 0), name = "") +
  scale_x_continuous(expand = c(0, 0), name = expression("S"[0] ~ "(mm)"), limits = c(-100, 300)) +
  scale_fill_manual(values = c("#D55E0050", "#0072B250"), labels = c("Diagnosed", "Prescribed"), name = "") +
  scale_color_manual(values = c("#D55E00", "#0072B2"), guide = "none") +
  scale_discrete_manual("point_color", values = c("#D55E00", "#0072B2"), guide = "none") +
  coord_cartesian(clip = "off") +
  guides(fill = guide_legend(
    override.aes = list(
      fill = c("#D55E00A0", "#0072B2A0"),
      color = NA, point_color = NA)
    )
  ) +
  theme_ridges(center = TRUE)
```