river ice analysis.Rmd

---
title: "001_Import_river_ice_data"
author: "Xiao Yang"
date: "6/3/2019"
output: 
  html_document: 
    df_print: tibble
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)


require(tidyverse)
require(foreach)

calculate_map_bounds_robinson = function(minlat = -55, maxlat = 80) {
  load("outputs/wrs_with_centroid.RData")
  latlon = wrs_lonlat %>% 
    filter(lat >= minlat, lat <= maxlat) %>% 
    st_transform(54030) %>% 
    st_coordinates()
  
  xymin = latlon[, 1:2] %>% apply(2, min)
  xymax = latlon[, 1:2] %>% apply(2, max)
  
  return(list(xymin = xymin, xymax = xymax))
}
```

# Import river ice dataset

```{r}
dat_dirs = c(
  "~/Google_Drive/river_ice_era_temp_dataset/",
  "~/Google_Drive_morphologee1/river_ice_era_temp_dataset/",
  "~/Google_Drive_yxiao/river_ice_era_temp_dataset/",
  "~/Google_Drive_morphologee1/river_ice_era_temp_dataset_cs15-25/"
)

files = dir(path = dat_dirs, full.names = T)

dat = foreach(i = 1:length(files), .combine = "bind_rows") %do% {
  read_csv(files[i]) %>% 
    filter(base_count != -9999)
}

save(dat, file = "outputs/river_ice_era5_temp_raw.RData")

load("outputs/river_ice_era5_temp_raw.RData", verbose = T)

datfil = dat %>% 
  mutate(
    nclearpixels = base_count - cloud_and_shadow_count
    ) %>% 
  filter(
    base_count > 0,
    base_count > cloud_and_shadow_count,
    !is.na(hillShadow_mean),
    nclearpixels > 0
    ) %>% 
  select(-`system:index`, -`.geo`) %>% 
  transmute(
    npixels = base_count,
    nclearpixels = nclearpixels,
    LANDSAT_SCENE_ID = LANDSAT_SCENE_ID,
    CLOUD_COVER = cloudscore,
    topo_shadow = hillShadow_mean,
    temp = pre30T2mMean_mean,
    ice_fraction = snow_ice_count / nclearpixels,
    cloud_fraction = cloud_and_shadow_count / base_count,
    PATH = as.integer(substr(LANDSAT_SCENE_ID, start = 4, stop = 6)),
    ROW = as.integer(substr(LANDSAT_SCENE_ID, start = 7, stop = 9)),
    ice = ice_fraction >= 0.5,
    year = as.integer(substr(LANDSAT_SCENE_ID, start = 10, stop = 13)),
    doy = as.integer(substr(LANDSAT_SCENE_ID, start = 14, stop = 16)),
    date = as.Date(doy - 1, origin = paste0(year, "-01-01")),
    period = factor(doy <= 30 | doy >= 210, levels = c(TRUE, FALSE), labels = c("Freeze-up", "Breakup"))
    )

load("outputs/wrs_with_centroid.RData", verbose = T)

datfil = datfil %>% 
  inner_join(wrs_lonlat %>% as.data.frame() %>% select(-geometry), by = c("PATH", "ROW"))

write_csv(datfil %>% select(date, river_ice_fraction = ice_fraction, cloud_fraction, topo_shadow, LANDSAT_SCENE_ID, PATH, ROW, N_river_pixel = npixels, N_clear_river_pixel = nclearpixels) %>% 
            mutate(PATH = as.integer(PATH),
                   ROW = as.integer(ROW),
                   N_river_pixel = as.integer(N_river_pixel),
                   N_clear_river_pixel = as.integer(N_clear_river_pixel)), path = "outputs/global_river_ice_dataset.csv")

write_csv(datfil %>% select(date, river_ice_fraction = ice_fraction, cloud_fraction, topo_shadow, LANDSAT_SCENE_ID, PATH, ROW, N_river_pixel = npixels, N_clear_river_pixel = nclearpixels) %>% 
            mutate(PATH = as.integer(PATH),
                   ROW = as.integer(ROW),
                   N_river_pixel = as.integer(N_river_pixel),
                   N_clear_river_pixel = as.integer(N_clear_river_pixel)) %>% sample_n(100), path = "outputs/global_river_ice_dataset_sample.csv")

dat_filterEffect = datfil %>% 
  mutate(month = factor(month(date), levels = 1:12, labels = month.abb)) %>% 
  mutate(valid = (nclearpixels >= 333) & (topo_shadow >= 0.95) & (cloud_fraction <= 0.25))

save(dat_filterEffect, file = "outputs/dat_filterEffect.RData")

require(lubridate)
datfil = datfil %>% 
  mutate(month = factor(month(date), levels = 1:12, labels = month.abb)) %>% 
  filter(
    nclearpixels >= 333,
    topo_shadow >= 0.95,
    cloud_fraction <= 0.25
    )

save(datfil, file = "outputs/river_ice_era_temp.RData")

load("outputs/river_ice_era_temp.RData", verbose = T)
```


# Climatology

## calculate map of winter river ice distribution

```{r}
season_table = tibble(
  month = month.abb,
  season = c(rep("Winter", 2), rep("Spring", 3), rep("Summer", 3), rep("Fall", 3), "Winter")
)

load("outputs/world_robinson.RData")
load("outputs/wrs_in_world_robinson.RData")

winter_mean_sf = datfil %>% 
  left_join(season_table) %>% 
  filter((season == "Winter" & lat >= 0) | (season == "Summer" & lat < 0)) %>% 
  group_by(PATH, ROW) %>% 
  summarise(
    mean_ice_fraction = mean(ice_fraction),
    median_ice_fraction = median(ice_fraction)) %>% 
  ungroup() %>% 
  left_join(wrs_in_world_robinson, by = c("PATH", "ROW")) %>% 
  st_as_sf()

xylim_robinson = calculate_map_bounds_robinson()
xymin = xylim_robinson[[1]]
xymax = xylim_robinson[[2]]

winter_mean_map = winter_mean_sf %>% 
  ggplot() +
  geom_sf(data = world_robinson, fill = "black", color = NA) +
  geom_sf(aes(fill = mean_ice_fraction * 100), color = NA) +
  geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
  coord_sf(crs = st_crs(54030), 
           xlim = c(xymin[1], xymax[1]), 
           ylim = c(xymin[2], xymax[2]),
           expand = T) +
  labs(fill = "River ice extent (length percentage)") +
  scale_fill_viridis_c(
    limits = c(0, 100),
    guide = guide_colorbar(
      label.theme = element_text(family = "sans", size = 5),
      title.theme = element_text(family = "sans", size = 5),
      direction = "horizontal",
      nbin = 10,
      draw.ulim = FALSE, 
      draw.llim = FALSE,
      title.position = "top",
      title.hjust = 0,
      barheight = unit(1, units = "mm"),
      barwidth = unit(25, units = "mm")
      )) +
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title = element_blank(),
        line = element_blank(),
        rect = element_blank(),
        panel.grid = element_blank(),
        legend.position = c(0.6, 0.090),
        panel.grid.major = element_line(color = "white", size = 0))

winter_mean_map

winter_mean_map %>% 
  ggsave(filename = "figs/final_figures/river_ice_winter_mean.pdf",
         device = "pdf",
         units = "mm",
         width = 120,
         height = 53,
         colormode = "rgb")

# winter_mean_map %>% 
#   ggsave(filename = "figs/river_ice_winter_mean.png",
#          device = "png",
#          width = 6,
#          height = 3,
#          dpi = "print")
```

## monthly river ice maps

```{r}
require(sf)

monthly_ice = datfil %>% 
  group_by(PATH, ROW, month) %>% 
  summarise(mean_ice_fraction = mean(ice_fraction),
            n = n(),
            npixels = max(npixels)) %>% 
  ungroup()

load("outputs/wrs_in_world_robinson.RData", verbose = T)
load("outputs/world_robinson.RData", verbose = T)

monthly_ice_sf = monthly_ice %>% 
  inner_join(wrs_in_world_robinson %>% select(PATH, ROW)) %>% 
  st_as_sf

for (i in 1:12) {
  monthly_fig = monthly_ice_sf %>%
    filter(month == month.abb[i]) %>%
    ggplot() +
    geom_sf(data = world_robinson, fill = "black", color = "black", lwd = 0.2) +
    geom_sf(aes(fill = mean_ice_fraction * 100), color = NA) +
    coord_sf(crs = st_crs(54030), 
           xlim = c(xymin[1], xymax[1]), 
           ylim = c(xymin[2], xymax[2]),
           expand = T) +
  labs(
    fill = "River ice extent (length percentage)",
    title = (month.name %>% toupper)[i]) +
  scale_fill_gradientn(
    colors = c("brown", "white", "cyan"),
    limits = c(0, 100),
    guide = guide_colorbar(
      direction = "horizontal",
      nbin = 10,
      draw.ulim = FALSE, 
      draw.llim = FALSE,
      title.position = "top",
      title.hjust = 0,
      barheight = unit(2, units = "mm"),
      barwidth = unit(50, units = "mm")
      )) +
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title = element_blank(),
        line = element_blank(),
        rect = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(size = 14, hjust = 0.5, face = "bold"),
        legend.position = c(0.65, 0.090),
        panel.grid.major = element_line(color = "white", size = 0),
        text = element_text(size = 8))

  monthly_fig %>%
    ggsave(
      filename = paste0("figs/monthly_river_ice_mean/", sprintf("%02d", i), ".png"),
      width = 5,
      height = 4,
      dpi = "print"
    )
}
```

## calculate monthly river ice stats

```{r}
## monthly stats for northern hemisphere
load("outputs/grwl_wrs_npixels.RData", verbose = T)
load("outputs/overlap_area_percentage_row.RData", verbose = T)

monthly_ice_global = datfil

all_pr = monthly_ice_global %>% dplyr::select(PATH, ROW) %>% distinct
global_total_rivern_area = all_pr %>% left_join(wrs_in_world_robinson) %>% st_as_sf %>% st_union() %>% st_area

month_stats = monthly_ice_global %>% 
  group_by(month, PATH, ROW) %>% 
  summarise(mrif = mean(ice_fraction),
            ni = n()) %>% 
  ungroup() %>% 
  left_join(grwl_wrs_npixels) %>% 
  left_join(overlap_area_percentage_row) %>% 
  filter(!is.na(n)) %>% 
  mutate(w = n * (1 - overlap_perc)) %>% 
  left_join(wrs_in_world_robinson) %>% 
  st_as_sf %>% 
  group_by(month) %>% 
  do({
    dat = . 
    rif_mean = weighted.mean(dat$mrif, dat$w)
    area_obs = dat %>% st_union() %>% st_area
    area_total = global_total_rivern_area
    tibble(area_obs, area_total, rif_mean)
  }) %>% 
  ungroup() %>% 
  mutate(percent_area_obs = area_obs / area_total * 100)

# ntiles = monthly_ice_nh %>% 
#   select(PATH, ROW) %>% 
#   distinct() %>% 
#   nrow()
# N = 7568214 # calculated from the total pixels in the asset: users/eeProject/GRWL_JRC_occGte90_widthGte90_chn1_lake0

monthly_stats_fig = month_stats %>% 
  ggplot() +
  geom_col(aes(x = month, y = rif_mean * 100), fill = "black") +
  # geom_point(aes(x = month, y = observed_percent * 100), color = "red") +
  geom_text(aes(x = month, y = rif_mean * 100, label = paste0(format(rif_mean * 100, digits = 1, hjust = 0), " (", format(percent_area_obs, digits = 0), "%)"), angle = 90, hjust = 0), nudge_x = 0, nudge_y = 1.5, size = 1.75, color = "black") +
  geom_text(aes(x = month, y = 0, label = month.abb, angle = 90, hjust = 1), nudge_x = 0, nudge_y = -1.5, size = 1.75, color = "black") +
  ylim(-13, 100) +
  labs(
    x = "",
    y = "NH river ice coverage\n(percentage of river network)"
  ) +
  theme_void()
  
monthly_stats_fig

monthly_stats_fig %>% 
  ggsave(filename = "figs/final_figures/global_monthly_river_ice_extent.pdf",
         device = "pdf",
         units = "mm",
         width = 30,
         height = 31, 
         colormode = "rgb")
```


# Changes

## calculate changes in river ice

```{r}
river_stats_history = datfil %>% 
  filter(date <= "1994-03-16" | date >= "2008-12-31") %>% 
  mutate(decade = factor(date <= "2005-01-01", levels = c(TRUE, FALSE), labels = c("1984-1994", "2008-2018")))

## exclude the yellow river tiles from the analysis
# load("outputs/yriver_tiles.RData", verbose = T)
# river_stats_history = river_stats_history %>%
#   left_join(yriver_tiles) %>%
#   filter(is.na(rn))

load("outputs/gridm.RData", verbose = T)
load("outputs/wrs2grid.RData", verbose = T)

common_pr_monthly = river_stats_history %>% 
  inner_join(wrs2grid %>% dplyr::select(grid_index = index, PATH, ROW)) %>% 
  group_by(month, grid_index) %>% 
  do({
    dat = .
    n1 = nrow(dat %>% dplyr::filter(decade == "1984-1994"))
    n2 = nrow(dat %>% dplyr::filter(decade == "2008-2018"))
    tibble(n1 = n1, n2 = n2)
  }) %>% 
  ungroup()

common_pr_monthly %>% 
  gather(key = "decade", value = "nobs", -c(month, grid_index)) %>% 
  ggplot() +
  geom_histogram(aes(nobs, fill = month), position = "stack", binwidth = 5) +
  facet_wrap(~decade)

## number of grids in each months that have at least 5 obs for each decade
common_pr_monthly %>% filter(n1 >= 5, n2 >= 5) %>% ggplot() + geom_bar(aes(x = month))

rif_diff = river_stats_history %>% 
  left_join(wrs2grid %>% dplyr::select(grid_index = index, PATH, ROW)) %>% 
  right_join(common_pr_monthly %>% filter(n1 >= 5, n2 >= 5)) %>% 
  group_by(month, grid_index) %>% 
  do({
    dat = .
    dat1 = dat %>% filter(decade == "1984-1994")
    dat2 = dat %>% filter(decade == "2008-2018")
    rif1 = dat1$ice_fraction
    rif2 = dat2$ice_fraction
    tibble(
      rif_diff = median(rif2) - median(rif1),
      n1 = dat1 %>% dplyr::select(PATH, ROW) %>% distinct() %>% nrow(),
      n2 = dat2 %>% dplyr::select(PATH, ROW) %>% distinct() %>% nrow())
  }) %>% 
  ungroup()
```

## calculate map of river ice change

```{r}
require(sf)
xylim_robinson = calculate_map_bounds_robinson()
xymin = xylim_robinson[[1]]
xymax = xylim_robinson[[2]]

load("outputs/world_robinson.RData", verbose = T)
gridm_rb = gridm %>% 
  st_transform(54030) %>% 
  st_intersection(world_robinson)

temp = rif_diff %>% 
  group_by(grid_index) %>% 
  summarise(mdif = mean(rif_diff),
            ntiles1 = sum(n1),
            ntiles2 = sum(n2),
            nmonths = n()) %>% 
  ungroup() %>% 
  left_join(gridm_rb %>% select(grid_index)) %>% 
  rename(geometry = g) %>% 
  st_as_sf

require(scales)
river_ice_diff = temp %>% 
  ggplot() +
  geom_sf(data = world_robinson, fill = "black", color = NA) +
  geom_sf(aes(fill = mdif * 100), color = NA) +
  geom_sf(data = gridm_rb, fill = NA, color = "grey30", size = 0.1) +
  geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
  coord_sf(crs = st_crs(54030), 
           xlim = c(xymin[1], xymax[1]), 
           ylim = c(xymin[2], xymax[2]),
           expand = T) +
  labs(fill = "River ice extent change \n(percentage point)") +
  scale_fill_gradient2(
    low = "red", high = "blue", mid = "grey", midpoint = 0,
    limits = c(-30, 10),
    guide = guide_colorbar(
      label.theme = element_text(family = "sans", size = 5),
      title.theme = element_text(family = "sans", size = 5),
      direction = "horizontal",
      nbin = 10,
      draw.ulim = FALSE, 
      draw.llim = FALSE,
      title.position = "top",
      title.hjust = 0,
      barheight = unit(1, units = "mm"),
      barwidth = unit(25, units = "mm")
      ), oob = squish) +
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        line = element_blank(),
        rect = element_blank(),
        panel.grid = element_blank(),
        legend.position = c(0.6, 0.080),
        panel.grid.major = element_line(color = "white", size = 0))

river_ice_diff

river_ice_diff %>% 
  ggsave(filename = "figs/final_figures/changes_in_median_historical_river_ice.pdf",
         device = "pdf",
         units = "mm",
         width = 120,
         height = 53, 
         colormode = "rgb")
```

## monthly map of river ice change

```{r}
rif_diff_sf = rif_diff %>% left_join(gridm_rb) %>% st_as_sf
require(scales)

for (i in 1:12) {
  
  month_river_ice_diff = rif_diff_sf %>% 
    filter(month == month.abb[i]) %>% 
    ggplot() +
    geom_sf(data = world_robinson, fill = "black", color = NA) +
    geom_sf(aes(fill = rif_diff * 100), color = NA) +
    geom_sf(data = gridm_rb, fill = NA, color = "grey30", size = 0.1) +
    geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
    coord_sf(crs = st_crs(54030), 
             xlim = c(xymin[1], xymax[1]), 
             ylim = c(xymin[2], xymax[2]),
             expand = T) +
    labs(fill = "River ice extent change \n(percentage point)",
         title = (month.name %>% toupper)[i]) +
    scale_fill_gradient2(
      low = "red", high = "blue", mid = "grey", midpoint = 0,
      limits = c(-20, 20),
      guide = guide_colorbar(
        direction = "horizontal",
        nbin = 10,
        draw.ulim = FALSE, 
        draw.llim = FALSE,
        title.position = "top",
        title.hjust = 0,
        barheight = unit(2, units = "mm"),
        barwidth = unit(50, units = "mm")
        ), oob = squish) +
    theme(axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          line = element_blank(),
          rect = element_blank(),
          plot.title = element_text(size = 14, hjust = 0.5, face = "bold"),
          panel.grid = element_blank(),
          legend.position = c(0.65, 0.080),
          panel.grid.major = element_line(color = "white", size = 0),
          text = element_text(size = 8))
  
  month_river_ice_diff
  
  month_river_ice_diff %>% 
    ggsave(filename = paste0("figs/monthly_river_ice_change/", sprintf("%02d", i), ".png"),
           width = 6,
           height = 3)
}

allmonth_river_ice_diff = rif_diff_sf %>% 
  # filter(month == month.abb[i]) %>% 
  ggplot() +
  geom_sf(data = world_robinson, fill = "black", color = NA) +
  geom_sf(aes(fill = rif_diff * 100), color = NA) +
  geom_sf(data = gridm_rb, fill = NA, color = "grey30", size = 0.05) +
  geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
  coord_sf(crs = st_crs(54030), 
           xlim = c(xymin[1], xymax[1]), 
           ylim = c(xymin[2], xymax[2]),
           expand = T) +
  labs(fill = "River ice extent change \n(percentage point)") +
  scale_fill_gradient2(
    low = "red", high = "blue", mid = "grey", midpoint = 0,
    limits = c(-20, 20),
    guide = guide_colorbar(
      direction = "horizontal",
      nbin = 10,
      draw.ulim = FALSE, 
      draw.llim = FALSE,
      title.position = "top",
      title.hjust = 0,
      barheight = unit(2, units = "mm"),
      barwidth = unit(50, units = "mm")
      ), oob = squish) +
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        line = element_blank(),
        rect = element_blank(),
        plot.title = element_text(size = 14, hjust = 0.5, face = "bold"),
        panel.grid = element_blank(),
        legend.position = "bottom",
        panel.grid.major = element_line(color = "white", size = 0),
        text = element_text(size = 6)) +
  facet_wrap(~month, ncol = 2)
  
  allmonth_river_ice_diff
  
  allmonth_river_ice_diff %>% 
    ggsave(filename = "figs/final_figures/historical_change_by_month.eps",
           units = "mm",
           width = 140,
           height = 240)
```

## monthly difference calculated based on 55 grid

```{r}
load("outputs/gridm.RData", verbose = T)

grid_global = gridm

rif_diff_global = rif_diff %>% 
  inner_join(
    grid_global %>% as.data.frame %>% dplyr::select(grid_index), by = "grid_index"
    )

allgrids_area = rif_diff_global %>% 
  dplyr::select(grid_index) %>% 
  distinct() %>% 
  left_join(gridm) %>% 
  st_as_sf %>% 
  st_union() %>% 
  st_area

monthly_change = rif_diff_global %>% 
  left_join(gridm %>% dplyr::select(grid_index)) %>% 
  st_as_sf() %>% 
  group_by(month) %>% 
  do({
    dat = .
    ice_fraction_diff = mean(dat$rif_diff)
    observed_percent_area = dat %>% st_union %>% st_area / allgrids_area
    tibble(ice_fraction_diff, observed_percent_area)
  }) %>% 
  ungroup()

decadal_change_monthly_obs_bar = monthly_change %>% 
  ggplot() +
  geom_col(aes(x = month, y = ice_fraction_diff * 100), fill = "black") +
  geom_text(aes(x = month, y = ice_fraction_diff * 100, label = paste0(format(ice_fraction_diff * 100, digits = 1, hjust = 0), " (", format(observed_percent_area * 100, digits = 0), "%)"), angle = -90), nudge_x = 0, nudge_y = -3, size = 1.75, color = "black") +
  geom_text(aes(x = month, y = 0, label = month.abb, angle = -90, hjust = 1), nudge_x = 0, nudge_y = 0.2, size = 1.75, color = "black") +
  scale_x_discrete(limits = rev(month.abb)) +
  scale_y_continuous(limits = c(-11, 2)) +
  theme_void()

decadal_change_monthly_obs_bar

decadal_change_monthly_obs_bar %>% 
  ggsave(filename = "figs/final_figures/decadal_change_monthly_obs_bar_global.pdf",
         device = "pdf",
         units = "mm",
         width = 30,
         height = 30, 
         colormode = "rgb")

summary(monthly_change %>% right_join(tibble(month = month.abb)))
```

## Image availability two decades

```{r}
xylim_robinson = calculate_map_bounds_robinson()
xymin = xylim_robinson[[1]]
xymax = xylim_robinson[[2]]

allpr = river_stats_history %>% 
  select(PATH, ROW) %>% 
  distinct()

allpr = allpr %>% 
  mutate(decade = FALSE) %>% 
  bind_rows(allpr %>% 
              mutate(decade = TRUE)) %>% 
  mutate(decade = factor(decade, levels = c(FALSE, TRUE), labels = c("1984-1994", "2008-2018")))

allpr_sf = allpr %>% 
  left_join(wrs_in_world_robinson %>% select(PATH, ROW), by = c("PATH", "ROW")) %>% 
  st_as_sf

image_availability_map = river_stats_history %>% 
  group_by(PATH, ROW, decade) %>% 
  count %>% 
  ungroup %>% 
  right_join(allpr_sf, by = c("PATH", "ROW", "decade")) %>% 
  st_as_sf %>% 
  ggplot() +
  geom_sf(data = world_robinson, fill = "black", color = NA) +
  geom_sf(aes(fill = n), color = NA) +
  geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
  coord_sf(crs = st_crs(54030), 
         xlim = c(xymin[1], xymax[1]), 
         ylim = c(xymin[2], xymax[2]),
         expand = T) +
  facet_wrap(~decade, ncol = 1) +
  scale_fill_viridis_c(
    limits = c(0, 180),
    guide = guide_colorbar(
      direction = "horizontal",
      nbin = 10,
      draw.ulim = FALSE, 
      draw.llim = FALSE,
      title.position = "top",
      title.hjust = 0.5,
      barheight = unit(2, units = "mm"),
      barwidth = unit(60, units = "mm")
      )
  ) +
  labs(fill = "Number of images") +
  theme(axis.text.y = element_blank(),
      axis.text.x = element_blank(),
      line = element_blank(),
      rect = element_blank(),
      panel.grid = element_blank(),
      legend.position = "bottom",
      panel.grid.major = element_line(color = "white", size = 0),
      text = element_text(size = 6))

image_availability_map

image_availability_map %>%
  ggsave(filename = "figs/final_figures/ED1a_image_availability_map.eps",
         units = "mm",
         width = 120,
         height = 132,
         dpi = 300,
         colormode = "rgb")
```

## Image availability bar chart

```{r}
load("outputs/dat_filterEffect.RData", verbose = T)

dat_filterEffect = dat_filterEffect %>% 
  filter(date <= "1994-03-16" | date >= "2008-12-31") %>% 
  mutate(decade = factor(date <= "2005-01-01", levels = c(TRUE, FALSE), labels = c("1984-1994", "2008-2018")))

validObs = dat_filterEffect %>% 
  group_by(month, decade) %>% 
  summarise(
    n = n(),
    nvalid = sum(valid),
    pvalid = sum(valid) / n
  ) %>% 
  ungroup()

valid_obs_fig = validObs %>% 
  ggplot() +
  geom_col(aes(x = month, y = pvalid, fill = decade), position = "dodge") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "",
    y = "Percent of successful acquisitions",
    fill = "Period"
  ) +
  coord_cartesian(expand = c(0, 0), ylim = c(0, 0.61)) +
  theme(text = element_text(size = 6))

valid_obs_fig

valid_obs_fig %>% 
  ggsave(
    filename = "figs/final_figures/ED1b_valid_obs_fig.eps",
    width = 120,
    height = 72,
    units = "mm",
    dpi = 300,
    colormode = "rgb"
  )
```


## Monthly image availability—two decades side-by-side

This map shows the spatial distribution of data that went into the monthly change calculation. Only the tiles that have at least three images in both decades were included in the calculation.

```{r}
xylim_robinson = calculate_map_bounds_robinson()
xymin = xylim_robinson[[1]]
xymax = xylim_robinson[[2]]

allpr = river_stats_history %>% 
  select(PATH, ROW) %>% 
  distinct() %>% 
  group_by(PATH, ROW) %>% 
  do({
    dat = .
    tmp = tibble(PATH = dat$PATH,
           ROW = dat$ROW,
           month = 1:12)
    tmp %>% 
      mutate(decade = FALSE) %>% 
      bind_rows(tmp %>% 
                  mutate(decade = TRUE)) %>% 
      mutate(decade = factor(decade, levels = c(FALSE, TRUE), labels = c("1984-1994", "2008-2018")))
  }) %>% 
  ungroup() %>% 
  mutate(month = factor(month, levels = 1:12, labels = month.abb))

allpr_sf = allpr %>% 
  left_join(wrs_in_world_robinson %>% select(PATH, ROW), by = c("PATH", "ROW")) %>% 
  st_as_sf

monthly_image_availability_map = river_stats_history %>% 
  mutate(month = factor(month, levels = 1:12, labels = month.abb)) %>% 
  group_by(PATH, ROW, decade, month) %>% 
  count %>% 
  ungroup %>% 
  mutate(gte3 = n >= 3) %>% 
  right_join(allpr_sf, by = c("PATH", "ROW", "decade", "month")) %>% 
  mutate(gte3cat = factor(gte3, levels = c(TRUE, FALSE, NA), 
                          labels = c(expression(n >= 3), expression(n < 3), "No data"), exclude = NULL)) %>% 
  # filter(month %in% c("Jan")) %>%
  st_as_sf() %>% 
  ggplot() +
  geom_sf(data = world_robinson, fill = "black", color = NA) +
  geom_sf(aes(fill = gte3cat), color = NA) +
  geom_sf(data = world_robinson, fill = NA, color = "black", size = 0.1) +
  coord_sf(crs = st_crs(54030), 
         xlim = c(xymin[1], xymax[1]), 
         ylim = c(xymin[2], xymax[2]),
         expand = T) +
  scale_fill_manual(
    values = c("green", "yellow", "red"),
    guide = guide_legend(
      direction = "vertical",
      title.position = "top",
      title.hjust = 0
      )
  ) +
  facet_grid(rows = month~decade) +
  labs(fill = "No. of images") +
  theme(axis.text.y = element_blank(),
      axis.text.x = element_blank(),
      axis.title = element_blank(),
      axis.ticks = element_blank(),
      line = element_blank(),
      rect = element_blank(),
      panel.grid = element_blank(),
      legend.position = "right",
      panel.grid.major = element_line(color = "white", size = 0),
      text = element_text(size = 10, face = "bold"))

monthly_image_availability_map

monthly_image_availability_map %>%
  ggsave(filename = "figs/monthly_image_availability_map_comp_two_decades.png",
         width = 7.5,
         height = 11,
         dpi = "print")
```

## Landsat sampling pattern's effect on spatial and temporal aggregation

```{r}
river_stats_history = datfil %>% 
  filter(date <= "1994-03-16" | date >= "2008-12-31") %>% 
  mutate(decade = factor(date <= "2005-01-01", levels = c(TRUE, FALSE), labels = c("1984-1994", "2008-2018")))

## exclude the yellow river tiles from the analysis
# load("outputs/yriver_tiles.RData", verbose = T)
# river_stats_history = river_stats_history %>%
#   left_join(yriver_tiles) %>%
#   filter(is.na(rn))

load("outputs/gridm.RData", verbose = T)
load("outputs/wrs2grid.RData", verbose = T)

common_pr_monthly = river_stats_history %>% 
  inner_join(wrs2grid %>% dplyr::select(grid_index = index, PATH, ROW)) %>% 
  group_by(month, grid_index) %>% 
  do({
    dat = .
    n1 = nrow(dat %>% dplyr::filter(decade == "1984-1994"))
    n2 = nrow(dat %>% dplyr::filter(decade == "2008-2018"))
    tibble(n1 = n1, n2 = n2)
  }) %>% 
  ungroup()

FUN1 = median
FUN2 = mean

monthly_diff = river_stats_history %>% 
  inner_join(wrs2grid %>% dplyr::select(grid_index = index, PATH, ROW)) %>% 
  inner_join(common_pr_monthly, by = c("month", "grid_index")) %>% 
  dplyr::select(lon, lat, doy, month, decade, river_ice_fraction = ice_fraction, grid_index, n1, n2)

sampling_pattern = monthly_diff %>% 
  filter(n1 >= 5,
         n2 >= 5) %>% 
  dplyr::select(lon, lat, doy, month, grid_index, decade, river_ice_fraction, n1, n2) %>% 
  distinct %>% 
  group_by(month, grid_index, decade) %>% 
  summarise(medianDOY = FUN1(doy),
            medianLat = FUN1(lat),
            medianLon = FUN1(lon),
            medianRif = median(river_ice_fraction),
            n1 = first(n1), 
            n2 = first(n2),
            meanDOY = FUN2(doy),
            meanLat = FUN2(lat),
            meanLon = FUN2(lon),
            meanRif = mean(river_ice_fraction)
            ) %>% 
  ungroup

meanDoyDiff = sampling_pattern %>% 
  dplyr::select(decade, meanDOY, grid_index, month) %>% 
  spread(data = ., key = decade, value = meanDOY) %>% 
  mutate(meanDOYDiff = `2008-2018` - `1984-1994`) %>% 
  dplyr::select(-`2008-2018`, -`1984-1994`)

meanLatDiff = sampling_pattern %>% 
  dplyr::select(decade, meanLat, grid_index, month) %>% 
  spread(data = ., key = decade, value = meanLat) %>% 
  mutate(meanLatDiff = `2008-2018` - `1984-1994`) %>% 
  dplyr::select(-`2008-2018`, -`1984-1994`)

meanRifDiff = sampling_pattern %>% 
  dplyr::select(decade, meanRif, grid_index, month) %>% 
  spread(data = ., key = decade, value = meanRif) %>% 
  mutate(meanRifDiff = `2008-2018` - `1984-1994`) %>% 
  dplyr::select(-`2008-2018`, -`1984-1994`)

mean_merged = meanDoyDiff %>% 
  inner_join(meanLatDiff) %>% 
  inner_join(meanRifDiff)

mean_merged = mean_merged %>% 
  gather(key = "sampling_type",
         value = "value",
         -c(meanRifDiff, grid_index, month))

## stats
# 1. mean and sd for difference in temporal and spatial sampling
mean_merged %>% 
  group_by(sampling_type) %>% 
  summarise(mean = mean(value),
            sd = sd(value))

mean_merged %>% 
  # filter(abs(meanRifDiff) >= 0.1) %>%
  group_by(sampling_type) %>% 
  do({
    dat = .
    temp = cor.test(x = dat$meanRifDiff, y = dat$value)
    tibble(r = temp$estimate, p = temp$p.value)
  }) %>% 
  ungroup()

meanDOYDiff_fig = mean_merged %>% 
  filter(sampling_type == "meanDOYDiff") %>% 
  ggplot() +
  geom_histogram(aes(x = value)) +
  xlim(-31, 31) +
  labs(x = "Difference of mean sampling time (day)",
       y = "Count") +
  theme(text = element_text(size = 6))

meanLatDiff_fig = mean_merged %>% 
  filter(sampling_type == "meanLatDiff") %>% 
  ggplot() +
  geom_histogram(aes(x = value)) +
  xlim(-5, 5) +
  labs(x = "Difference of mean latitude (º)",
       y = "Count") +
  theme(text = element_text(size = 6))

mean_doy_rif_diff_fig = mean_merged %>% 
  filter(sampling_type == "meanDOYDiff") %>% 
  ggplot() +
  geom_point(aes(x = value, y = meanRifDiff * 100), size = 0.2) +
  # facet_wrap(~month) +
  labs(x = "Difference of mean sampling time (day)",
       y = "Difference of mean river ice extent\n(percentage point)") +
  theme(text = element_text(size = 6))

mean_lat_rif_diff_fig = mean_merged %>% 
  filter(sampling_type == "meanLatDiff") %>% 
  ggplot() +
  geom_point(aes(x = value, y = meanRifDiff * 100), size = 0.2) +
  # facet_wrap(~month) +
  labs(x = "Difference of mean latitude (º)",
       y = "Difference of mean river ice extent\n(percentage point)") +
  theme(text = element_text(size = 6))

meanDOYDiff_fig %>% ggsave(filename = "figs/final_figures/ED7a_meanDOYDiff_fig.eps", width = 60, height = 60, dpi = 300, colormode = "rgb", units = "mm")
meanLatDiff_fig %>% ggsave(filename = "figs/final_figures/ED7c_meanLatDiff_fig.eps", width = 60, height = 60, dpi = 300, colormode = "rgb", units = "mm")
mean_doy_rif_diff_fig %>% ggsave(filename = "figs/final_figures/ED7b_mean_doy_rif_diff_fig.eps", width = 60, height = 60, dpi = 300, colormode = "rgb", units = "mm")
mean_lat_rif_diff_fig %>% ggsave(filename = "figs/final_figures/ED7d_mean_lat_rif_diff_fig.eps", width = 60, height = 60, dpi = 300, colormode = "rgb", units = "mm")
```


# Models

## global logistic regression model 

```{r}
model_input_global = datfil %>% 
  mutate(
    n_ice = floor(nclearpixels * ice_fraction),
    n_water = nclearpixels - n_ice) %>% 
  filter(lat >= 23.5)

glr1 = glm(cbind(n_ice, n_water) ~ 1 + temp + period, data = model_input_global, family = binomial)
glr2 = glm(cbind(n_ice, n_water) ~ 1 + temp + temp:period, data = model_input_global, family = binomial)
glr3 = glm(cbind(n_ice, n_water) ~ 1 + temp, data = model_input_global, family = binomial)

newData = tibble(temp = seq(-40, 40, by = 1), period = "Breakup") %>% 
  bind_rows(tibble(temp = seq(-40, 40, by = 1), period = "Freeze-up"))

pred1 = newData %>% mutate(ice_fraction = predict(glr1, newdata = newData, type = "response", se.fit = F))
pred2 = newData %>% mutate(ice_fraction = predict(glr2, newdata = newData, type = "response", se.fit = F))
pred3 = newData %>% mutate(ice_fraction = predict(glr3, newdata = newData, type = "response", se.fit = F))

summary(glr2)

model_and_data_fig = model_input_global %>% 
  # sample_frac(0.1) %>% 
  ggplot() +
  geom_hex(aes(x = temp, y = ice_fraction, fill = log10(..count..))) +
  # geom_point(aes(x = temp, y = ice_fraction), size = 0.2, alpha = 0.2) +
  geom_line(data = pred2, aes(x = temp, y = ice_fraction, color = period), lwd = 1) +
  scale_fill_gradientn(colours = rev(grey(level = seq(0, 0.97, length = 100))), breaks = 0:4, labels = as.character(10^(0:4))) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "30-day prior mean SAT (ºC)",
    y = "River ice extent",
    fill = "No. of\nobservations",
    color = "Period"
  ) +
  facet_wrap(~period, nrow = 2) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.25, linetype = 1),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        strip.text = element_text(size = 6, family = "sans"),
        text = element_text(size = 6, family = "sans"),
        strip.background = element_blank())

model_and_data_fig

model_and_data_fig %>% 
  ggsave(filename = "figs/final_figures/fig2a_separate.pdf",
         device = "pdf",
         units = "mm",
         width = 89,
         height = 107,
         colormode = "rgb")

data_fig = model_input_global %>% 
  mutate(odds = n_ice / n_water) %>% 
  # sample_n(1000) %>% 
  ggplot() +
  geom_point(aes(x = temp, y = log(odds)), size = 0.2, alpha = 0.2) +
  # scale_y_continuous(labels = scales::percent) +
  labs(
    x = "30-day prior mean SAT (ºC)",
    y = "River ice extent (log)"
  ) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.25, linetype = 1),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        strip.text = element_text(size = 12, face = "bold"),
        strip.background = element_blank())

data_fig

data_fig %>% 
  ggsave(filename = "figs/fig_data_temp_fig.png",
         width = 5.5,
         height = 3,
         dpi = "print")

model_errors_overall = model_input_global %>% 
  mutate(
    ice_fraction_glr = predict(glr2, newdate = ., type = "response", se.fit = F),
    ice_fraction_0deg = temp <= 0, 
    diff_glr = ice_fraction_glr - ice_fraction,
    diff_0deg = ice_fraction_glr - ice_fraction_0deg) %>% 
  dplyr::select(
    diff_glr,
    diff_0deg
  ) %>% 
  summarise(
    rmse_glr = sqrt(mean(diff_glr^2)),
    rmse_0deg = sqrt(mean(diff_0deg^2)),
    mbs_glr = mean(diff_glr),
    mbs_0deg = mean(diff_0deg)
  )

model_errors_overall

model_errors = model_input_global %>% 
  mutate(
    ice_fraction_glr = predict(glr2, newdate = ., type = "response", se.fit = F),
    ice_fraction_0deg = temp <= 0) %>% 
  dplyr::select(
    ice_fraction_obs = ice_fraction,
    ice_fraction_glr,
    ice_fraction_0deg,
    temp
  ) %>% 
  mutate(temp_grp = cut(temp, breaks = seq(-48, 42, by = 2)),
         model_0deg = ice_fraction_0deg - ice_fraction_obs,
         model_glr = ice_fraction_glr - ice_fraction_obs) %>% 
  dplyr::select(-ice_fraction_obs, -ice_fraction_glr, -ice_fraction_0deg) %>% 
  gather(key = "Model", value = "diff", -c(temp, temp_grp)) %>% 
  group_by(temp_grp, Model) %>% 
    summarise(
      rmse = sqrt(mean(diff^2)),
      mae = mean(abs(diff)),
      mbs = mean(diff),
      temp_mean = mean(temp)) %>% 
  ungroup()

model_comp = model_errors %>% 
  gather(key = "error_type", value = "error_value", -c(temp_mean, Model, temp_grp)) %>% 
  filter(error_type != "mae") %>% 
  mutate(Model = factor(Model, levels = c("model_0deg", "model_glr"), labels = c("0 ºC isotherm", "Logistic regression")),
         error_type = factor(error_type, levels = c("mbs", "rmse"), labels = c("Mean bias", "RMSE"))) %>% 
  ggplot() +
  geom_area(aes(x = temp_mean, y = error_value, fill = Model), stat = "identity", alpha = 0.5, size = 1.5, position = "identity") +
  geom_line(aes(x = temp_mean, y = error_value, group = Model), color = "black", alpha = 1, lwd = 0.5) +
  facet_wrap(~error_type) +
  scale_x_continuous(limits = c(-30, 30), breaks = seq(-20, 20, by = 10)) +
  scale_y_continuous(limits = c(-0.20, 0.8), labels = scales::percent) +
  coord_cartesian(expand = c(0, 0)) +
  labs(
    x = "30-day prior mean SAT (ºC)",
    y = "Ice cover extent error (%)"
  ) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.25, linetype = 1),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        text = element_text(size = 6, family = "sans"),
        strip.background = element_blank())

model_comp

model_comp %>% 
  ggsave(filename = "figs/final_figures/fig2b_ylabel_corrected.pdf",
         device = "pdf",
         units = "mm",
         width = 89,
         height = 45,
         colormode = "rgb")

model_errors %>% 
  group_by(Model) %>% 
  summarise(rmse_max = max(rmse),
            mbs_max = max(mbs)) %>% 
  ungroup()

(0.681 - 0.356) / (0.681)

(0.584 - 0.0698) / (0.584)

## model errors between -10 and 10
model_errors = model_input_global %>% 
  mutate(
    ice_fraction_glr = predict(glr2, newdate = ., type = "response", se.fit = F),
    ice_fraction_0deg = temp <= 0) %>% 
  dplyr::select(
    ice_fraction_obs = ice_fraction,
    ice_fraction_glr,
    ice_fraction_0deg,
    temp
  ) %>% 
  filter(temp >= -10 & temp <= 10) %>% 
  mutate(model_0deg = ice_fraction_0deg - ice_fraction_obs,
         model_glr = ice_fraction_glr - ice_fraction_obs) %>% 
  dplyr::select(-ice_fraction_obs, -ice_fraction_glr, -ice_fraction_0deg) %>% 
  gather(key = "Model", value = "diff", -c(temp)) %>% 
  group_by(Model) %>% 
    summarise(
      rmse = sqrt(mean(diff^2)),
      mae = mean(abs(diff)),
      mbs = mean(diff),
      temp_mean = mean(temp)) %>% 
  ungroup()

model_errors

(0.313 - 0.22) / 0.313

(0.0814 - 0.0106) / 0.0814
  

```

## compare modeled and observed monthly stats

```{r}
## monthly stats for northern hemisphere
load("outputs/grwl_wrs_npixels.RData", verbose = T)
load("outputs/overlap_area_percentage_row.RData", verbose = T)

monthly_ice_nh = datfil %>% 
  filter(lat >= 0) %>% 
  mutate(ice_fraction_model = predict(glr2, newdata = ., type = "response", se.fit = F)) %>% 
  dplyr::select(month, PATH, ROW, ice_fraction_model, ice_fraction_observed = ice_fraction) %>% 
  gather(key = "source", value = "ice_fraction", -c(month, PATH, ROW))

all_pr = monthly_ice_nh %>% select(PATH, ROW) %>% distinct
nh_total_rivern_area = all_pr %>% left_join(wrs_in_world_robinson) %>% st_as_sf %>% st_union() %>% st_area

month_stats = monthly_ice_nh %>% 
  group_by(month, PATH, ROW, source) %>% 
  summarise(mrif = mean(ice_fraction),
            ni = n()) %>% 
  ungroup() %>% 
  left_join(grwl_wrs_npixels) %>% 
  left_join(overlap_area_percentage_row) %>% 
  filter(!is.na(n)) %>% 
  mutate(w = n * (1 - overlap_perc)) %>% 
  left_join(wrs_in_world_robinson) %>% 
  st_as_sf %>% 
  group_by(month, source) %>% 
  do({
    dat = . 
    rif_mean = weighted.mean(dat$mrif, dat$w)
    area_obs = dat %>% st_union() %>% st_area
    area_total = nh_total_rivern_area
    tibble(area_obs, area_total, rif_mean, n = sum(dat$ni))
  }) %>% 
  ungroup() %>% 
  mutate(percent_area_obs = area_obs / area_total * 100)

monthly_stats_comp_fig = month_stats %>% 
  mutate(source = factor(source, 
                         levels = c("ice_fraction_model", "ice_fraction_observed"),
                         labels = c("Modeled", "Observed"))) %>% 
  ggplot() +
  geom_point(aes(x = month, y = rif_mean, color = source, size = n)) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Month",
       y = "River ice fraction",
       color = "",
       size = "Number of images") +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.5, linetype = 1),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        strip.text = element_text(size = 12, face = "bold"),
        strip.background = element_blank())
  
monthly_stats_comp_fig

monthly_stats_comp_fig %>% 
  ggsave(filename = "figs/fig2c.png",
         width = 6,
         height = 3,
         dpi = "print")

## global stats
nh_comp = datfil %>% 
  filter(lat >= 0) %>% 
  mutate(ice_fraction_model = predict(glr2, newdata = ., type = "response", se.fit = F)) %>% 
  dplyr::select(month, PATH, ROW, ice_fraction_model, ice_fraction_observed = ice_fraction, temp)

nh_comp %>% 
  mutate(
    diff = ice_fraction_model - ice_fraction_observed,
    ) %>% 
  summarise(
    rmse = sqrt(mean(diff^2)),
    mae = mean(abs(diff)),
    mbs = mean(diff)
  )
```


## associate data with basins

```{r}
load("outputs/wrs_pr_with_hbID.RData", verbose = T)
load("outputs/hydrobasin_hn_pfaf_lvl_3_simplified_0.1.RData", verbose = T)
require(sf)

dathb = datfil %>% 
  inner_join(wrs_pr_with_hbID %>% mutate(HYBAS_ID = as.character(HYBAS_ID)))
  
ice_affected_basins = dathb %>% 
  group_by(HYBAS_ID) %>% 
  summarise(q95 = quantile(ice_fraction, probs = 0.95),
            nobs = n()) %>% 
  ungroup() %>% 
  filter(q95 >= 0.5) %>% 
  select(HYBAS_ID, q95, nobs) %>% 
  inner_join(hb_sim %>% mutate(HYBAS_ID = as.character(HYBAS_ID)), by = "HYBAS_ID") %>% 
  st_as_sf

remove(hb_sim)

require(mapview)
mapview::mapView(ice_affected_basins["nobs"])
```

## create model input (five fold validation)

```{r}
model_input = dathb %>% 
  right_join(ice_affected_basins %>% as.data.frame() %>% select(HYBAS_ID)) %>%
  mutate(n_ice = floor(nclearpixels * ice_fraction),
         n_water = nclearpixels - n_ice)

save(model_input, file = "output/model_input_era5temp.RData")

load("outputs/model_input_era5temp.RData", verbose = T)
```

## plot basin wise temp-ice relationship

```{r}
basin_ids = model_input %>% 
  select(HYBAS_ID) %>% 
  distinct() %>% 
  pull(HYBAS_ID)

for (i in basin_ids) {
  fig = model_input %>% 
    filter(HYBAS_ID == i) %>% 
    select(ice_fraction, temp, period) %>% 
    ggplot() +
    geom_point(aes(x = temp, y = ice_fraction), size = 0.5, alpha = 0.7) +
    geom_smooth(aes(x = temp, y = ice_fraction), method = glm, method.args = list(family = quasibinomial)) +
    labs(
      title = i
    )
  
  fig %>% 
    ggsave(
      filename = paste("figs/basinwise-ice-temp-relationship/", i, ".png", sep = ""),
      width = 8,
      height = 4
    )
}
```



### 0 deg model

```{r}
model_metric_0deg = model_input %>% 
  transmute(
    nclearpixels = nclearpixels,
    temp = temp,
    observed = ice_fraction,
    predicted = as.integer(temp <= 0),
    basin_id = HYBAS_ID,
    period = period,
    model = "0 degree threshold"
  )
```

### global logreg model

```{r}
gmodel = model_input %>% 
  glm(cbind(n_ice, n_water) ~ 1 + temp, data = ., family = binomial)

summary(gmodel)

model_metric_global = model_input %>% 
  transmute(
    nclearpixels = nclearpixels,
    temp = temp,
    observed = ice_fraction,
    predicted = predict(gmodel, newdata = ., type = "response"),
    basin_id = HYBAS_ID,
    period = period,
    model = "Global logreg model"
  )
```

### per basin model

```{r}
basins = model_input %>% 
  select(HYBAS_ID) %>% 
  distinct() %>% 
  pull()

basin_models = list()
for (i in 1:length(basins)) {
  print(i)
  print(basins[i])
  basin_input = model_input %>% 
    filter(
      HYBAS_ID == basins[i]
    )
  basin_models[[i]] = glm(formula = cbind(n_ice, n_water) ~ 1 + temp, data = basin_input, family = binomial)
  
  temp = basin_input %>% 
      transmute(
        nclearpixels = nclearpixels,
        temp = temp,
        observed = ice_fraction,
        predicted = predict(basin_models[[i]], newdata = ., type = "response"),
        basin_id = HYBAS_ID,
        period = period,
        model = "Basin logreg models"
        )
  
  temp1 = tibble(
    HYBAS_ID = basins[i],
    basin_inter = basin_models[[i]]$coefficients[1],
    basin_slope = basin_models[[i]]$coefficients[2])
  
  if (i == 1) {
    model_metric_basins = temp
    model_para = temp1
  } else if (i != 1) {
    model_metric_basins = bind_rows(model_metric_basins, temp)
    model_para = bind_rows(model_para, temp1)
  }
}
```

### mixed effect model

```{r}
require(lme4)

formula = "cbind(n_ice, n_water) ~ 1 + temp + (1 + temp | HYBAS_ID)"
mem = glmer(formula = formula, data = model_input, family = "binomial", 
            nAGQ = 0)

summary(mem)

model_metric_mem = model_input %>% 
  transmute(
    nclearpixels = nclearpixels,
    temp = temp,
    observed = ice_fraction,
    predicted = predict(mem, newdata = ., type = "response"),
    basin_id = HYBAS_ID,
    period = period,
    model = "Mixed effect model"
  )
```

## merged model function

```{r}
merged_model = function(model_input, model_types) {
  require(lme4)
  ## model 1
  print("building 0 degree threshold model...")
  model_error = model_input %>% 
      transmute(
        nclearpixels = nclearpixels,
        temp = temp,
        observed = ice_fraction,
        predicted = as.integer(temp <= 0),
        basin_id = HYBAS_ID,
        period = period,
        model = "0 degree threshold"
      )
  
  ## model 2
  print("building Global logreg model...")
  gmodel = model_input %>% 
    glm(cbind(n_ice, n_water) ~ 1 + temp + temp:period, data = ., family = binomial)
  
  model_metric_global = model_input %>% 
    transmute(
      nclearpixels = nclearpixels,
      temp = temp,
      observed = ice_fraction,
      predicted = predict(gmodel, newdata = ., type = "response"),
      basin_id = HYBAS_ID,
      period = period,
      model = "Global logreg model"
    )
  
  model_error = model_error %>% 
    bind_rows(model_metric_global)
  
  ## model 3
  print("building Basin logreg models...")
  basins = model_input %>% 
    select(HYBAS_ID) %>% 
    distinct() %>% 
    pull()
  
  basin_models = list()
  for (i in 1:length(basins)) {
    print(i)
    print(basins[i])
    basin_input = model_input %>% 
      filter(
        HYBAS_ID == basins[i]
      )
    basin_models[[i]] = glm(formula = cbind(n_ice, n_water) ~ 1 + temp + temp:period, data = basin_input, family = binomial)
    temp = basin_input %>% 
        transmute(
          nclearpixels = nclearpixels,
          temp = temp,
          observed = ice_fraction,
          predicted = predict(basin_models[[i]], newdata = ., type = "response"),
          basin_id = HYBAS_ID,
          period = period,
          model = "Basin logreg models"
          )
    if (i == 1) {
      model_metric_basins = temp
    } else if (i != 1) {
      model_metric_basins = bind_rows(model_metric_basins, temp)
    }
  }
  
  model_error = model_error %>% 
    bind_rows(model_metric_basins)
  
  ## model 4
  print("building Mixed effect model...")
  formula = "cbind(n_ice, n_water) ~ 1 + temp + temp:period + (1 + temp | HYBAS_ID)"
  mem = glmer(formula = formula, data = model_input, family = "binomial", 
              nAGQ = 0)
  
  model_metric_mem = model_input %>% 
    transmute(
      nclearpixels = nclearpixels,
      temp = temp,
      observed = ice_fraction,
      predicted = predict(mem, newdata = ., type = "response"),
      basin_id = HYBAS_ID,
      period = period,
      model = "Mixed effect model"
    )
  
  model_error = model_error %>% 
    bind_rows(model_metric_mem)
  
  ## output result 
  return(model_error)
}

test = merged_model(model_input)

test %>% 
  mutate(dff = predicted - observed) %>% 
  group_by(model) %>% 
  summarise(
    rmse = sqrt(mean(dff^2)),
    mbs = mean(dff)
  ) %>% 
  ungroup

test %>% 
  mutate(dff = predicted - observed) %>% 
  group_by(model) %>% 
  summarise(
    rmse = sqrt(mean(dff^2)),
    mbs = mean(dff)
  ) %>% 
  ungroup %>% 
  ggplot() +
  geom_col(aes(x = model, y = rmse, fill = model))
```


## comparing model error metrics

```{r}
common_basins = model_metric_mem %>% select(basin_id) %>% distinct() %>% 
  inner_join(model_metric_basins %>% select(basin_id) %>% distinct(), by = "basin_id") %>% 
  inner_join(model_metric_global %>% select(basin_id) %>% distinct(), by = "basin_id") %>% 
  inner_join(model_metric_0deg %>% select(basin_id) %>% distinct(), by = "basin_id") %>% 
  distinct()

model_metric = bind_rows(
  model_metric_0deg,
  model_metric_global,
  model_metric_basins,
  model_metric_mem
) %>% 
  right_join(common_basins)

errors = model_metric %>% 
  mutate(
    temp_grp = cut(temp, breaks = seq(-48, 42, by = 1), include.lowest = T),
    res = predicted - observed
    ) %>% 
  group_by(temp_grp, model) %>% 
  summarise(
    mae = weighted.mean(abs(res), w = nclearpixels),
    mbs = weighted.mean(res, w = nclearpixels),
    rmse = sqrt(weighted.mean(res^2, w = nclearpixels^2)),
    temp_mean = mean(temp)
  ) %>% 
  ungroup()

errors = errors %>% 
  gather(key = "metric", value = "error", -c(temp_mean, model, temp_grp))

errors %>% 
  group_by(metric, model) %>% 
  summarise(max_value = max(error),
            min_value = min(error)) %>% 
  ungroup() %>% 
  arrange(metric, max_value)

errors_fig = errors %>% 
  ggplot() +
  geom_line(aes(x = temp_mean, y = error * 100, color = model), alpha = 1) +
  facet_wrap(~metric) +
  ylim(-40, 100) +
  labs(
    x = "Temperature",
    y = "Error (percentage point)",
    color = "Models"
  )

errors_fig

errors_fig %>% 
  ggsave(
    filename = "figs/errors_fig.png",
    width = 5,
    height = 4
  )


errors_overall = model_metric %>% 
  filter(temp >= -10, temp <= 10) %>% 
  mutate(
    res = predicted - observed
    ) %>% 
  group_by(model) %>% 
  summarise(
    mae = weighted.mean(abs(res), w = nclearpixels),
    mbs = weighted.mean(res, w = nclearpixels),
    rmse = sqrt(weighted.mean(res^2, w = nclearpixels^2))
  ) %>% 
  ungroup()

errors_overall %>% 
  gather(key = "metric", value = "error", -model) %>% 
  group_by(metric) %>% 
  arrange(error) %>% 
  ungroup() %>% 
  ggplot() +
  geom_bar(aes(x = metric, y = error, fill = model), stat = "identity", position = "dodge", color = "black") +
  theme_classic()

```

## export model parameters

```{r}
randomEffect = lme4::ranef(mem, which = "HYBAS_ID", condVar = TRUE)

lattice::dotplot(randomEffect)

randomEffect_tbl = tibble(ranef_intercept = randomEffect[[1]]$`(Intercept)`,
                          ranef_temp = randomEffect[[1]]$temp,
                          HYBAS_ID = row.names(randomEffect[[1]]))

load("output/hydrobasin_hn_pfaf_lvl_3_simplified_0.1.RData", verbose = T)
hb = hb_sim %>% 
  mutate(HYBAS_ID = as.character(HYBAS_ID))

basinwise_random_effect = randomEffect_tbl %>% 
  inner_join(hb %>% select(HYBAS_ID), by = "HYBAS_ID") %>% 
  st_as_sf

# export shp file to be imported into GEE
global_slope = -0.37612
global_intercept = -1.28557

export = basinwise_random_effect %>% 
  transmute(
    intcpt = ranef_intercept + global_intercept,
    temp = ranef_temp + global_slope, 
    HYBAS_ID = HYBAS_ID) %>% 
  st_as_sf() %>% 
  st_simplify(preserveTopology = T, dTolerance = 0.05)

## add parameters for per basin mdoels
export = model_para %>% 
  right_join(export)

st_write(export, dsn = "output/basinwise_random_effect_two_models.shp")
```

## compare modeled month ice extent to Landsat derived value

```{r}
load("outputs/model_input_era5temp.RData", verbose = T)
load("outputs/river_ice_era_temp.RData", verbose = T)
require(lubridate)
require(sf)

datfil_tropical = datfil %>% 
  inner_join(wrs_lonlat %>% filter(lat >= 0, lat <= 23.5) %>% as.data.frame() %>% select(PATH, ROW))

modeled_ice = model_input %>% 
  mutate(ice_fraction_mem = predict(mem, newdata = ., type = "response")) %>% 
  select(-HYBAS_ID, -n_ice, -n_water) %>% 
  bind_rows(datfil_tropical %>% 
              mutate(ice_fraction_mem = 0))

datfil = datfil %>% 
  inner_join(modeled_ice %>% select(LANDSAT_SCENE_ID, ice_fraction_mem), by = "LANDSAT_SCENE_ID")

datfil = datfil %>% 
  select(month, PATH, ROW, npixels, ice_fraction, ice_fraction_mem, date) %>% 
  gather(key = "source", value = "ice_fraction", -c(month, PATH, ROW, npixels, date))

monthly_ice = datfil %>% 
  mutate(month = month(date)) %>% 
  group_by(PATH, ROW, month, source) %>% 
  summarise(
    mean_ice_fraction = mean(ice_fraction),
    rif_quantile10 = quantile(ice_fraction, probs = 0.1),
    rif_quantile90 = quantile(ice_fraction, probs = 0.9),
    n_image = n(),
    ntotalpixels = sum(npixels)) %>% 
  ungroup()

load("outputs/wrs_in_world_robinson.RData", verbose = T)
load("outputs/world_robinson.RData", verbose = T)

## monthly stats for northern hemisphere
load("outputs/wrs_with_centroid.RData", verbose = T)
monthly_ice_nh = monthly_ice %>% 
  inner_join(wrs_lonlat %>% filter(lat >= 0) %>% as.data.frame() %>% select(PATH, ROW))

ntiles = monthly_ice_nh %>% 
  select(PATH, ROW) %>% 
  distinct() %>% 
  nrow()
# N = 7568214 # calculated from the total pixels in the asset: users/eeProject/GRWL_JRC_occGte90_widthGte90_chn1_lake0

monthly_stats = monthly_ice_nh %>% 
  group_by(month, source) %>% 
  summarise(
    mean_ice_fraction = weighted.mean(mean_ice_fraction, w = ntotalpixels),
    rif_quantile10 = weighted.mean(rif_quantile10, w = ntotalpixels),
    rif_quantile90 = weighted.mean(rif_quantile90, w = ntotalpixels),
    n_image = sum(n_image),
    observed_percent = n() / ntiles
    ) %>% 
  ungroup()

# diff_rif = (river_stats_monthly_mean %>% transmute(diff = rif_wmean - rif_wmean_modeled))$diff * 100
# mae = mean(abs(diff_rif))
# mbs = mean(diff_rif)
# rmse = sqrt(mean(diff_rif^2))

monthly_stats %>% 
  ggplot() +
  geom_errorbar(aes(x = month, ymin = rif_quantile10 * 100, ymax = rif_quantile90 * 100)) +
  geom_point(aes(x = month, y = mean_ice_fraction * 100, size = n_image, color = source)) +
  labs(x = "Month",
       y = "River ice fraction (%)",
       size = "Number of images")
```


# Prediction

## decide which climate models to use

```{r}
dat45 = read_csv("data/tempdiffclimatemodel_rcp45.csv") %>% dplyr::select(-aoi)
dat85 = read_csv("data/tempdiffclimatemodel_rcp85.csv") %>% dplyr::select(-aoi)

dat = dat45 %>% 
  bind_rows(dat85) %>% 
  dplyr::select(
    diff_mean = dailyMeanDiff_mean,
    diff_median = dailyMeanDiff_median,
    model,
    scenario)

dat_mean = dat %>% 
  dplyr::select(diff_mean, model, scenario) %>% 
  spread(scenario, diff_mean)

dat_mean %>% 
  ggplot() +
  geom_point(aes(x = rcp45, y = rcp85, color = model))

mode_name_ascend = dat %>% 
  filter(scenario == "rcp85") %>% 
  arrange(diff_mean) %>% 
  pull(model)

comp_fig = dat %>% 
  mutate(model = factor(model, levels = mode_name_ascend, labels = mode_name_ascend, ordered = T)) %>% 
  ggplot() +
  geom_bar(aes(x = model, y = diff_mean, fill = scenario), stat = "identity", position = "dodge") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Model",
       y = "Difference in mean SAT (ºC)",
       title = "Difference in mean surface air tempearture",
       subtitle = "between 2006-2036 and 2069-2099",
       fill = "RCPs")

comp_fig

comp_fig %>% 
  ggsave(
    filename = "figs/final_figures/Difference in mean surface tempearture climate models.tiff",
    units = "mm",
    width = 180,
    height = 90,
    dpi = 300
  )
```

## monthly predicted change map 

### export shp file for clipping the exported raster global maps to only include land

```{r}
latmin = 0
latmax = 90
extent_coords = rbind(c(-180, latmin), c(180, latmin), c(180, latmax), c(-180, latmax), c(-180, latmin))
extentply = st_polygon(list(extent_coords))
extentsfc_nh = st_sfc(extentply, crs = 4326)

latmin = -90
latmax = 0
extent_coords = rbind(c(-180, latmin), c(180, latmin), c(180, latmax), c(-180, latmax), c(-180, latmin))
extentply = st_polygon(list(extent_coords))
extentsfc_sh = st_sfc(extentply, crs = 4326)

world = st_as_sf(maps::map(database = "world", fill = T, plot = F))
world = world %>% st_buffer(0.0)

# world_clipped_nh_land = world %>% 
#   st_intersection(extentsfc_nh)
# 
# world_clipped_sh_land = world %>% 
#   st_intersection(extentsfc_sh)

world_clipped_nh_land = world %>% 
  st_difference(extentsfc_sh)

world_clipped_sh_land = world %>%
  st_difference(extentsfc_nh)

plot(world_clipped_nh_land %>% st_transform(crs = 3995))
plot(world_clipped_sh_land)

st_write(world_clipped_nh_land, dsn = "~/Google_Drive/monthly_map_river_ice_future_change_global/world_clipped_nh.shp")
st_write(world_clipped_sh_land, dsn = "~/Google_Drive/monthly_map_river_ice_future_change_global/world_clipped_sh.shp")
```

### export shp file for clipping the exported raster global maps to only include GRWL covered land

```{r}
latmin = 0
latmax = 90
extent_coords = rbind(c(-180, latmin), c(180, latmin), c(180, latmax), c(-180, latmax), c(-180, latmin))
extentply = st_polygon(list(extent_coords))
extentsfc_nh = st_sfc(extentply, crs = 4326)

latmin = -90
latmax = 0
extent_coords = rbind(c(-180, latmin), c(180, latmin), c(180, latmax), c(-180, latmax), c(-180, latmin))
extentply = st_polygon(list(extent_coords))
extentsfc_sh = st_sfc(extentply, crs = 4326)

# load("outputs/wrs_in_world_robinson.RData", verbose = T)
load("outputs/grwl_by_wrstile.RData", verbose = T)
load("outputs/wrs_clipped_to_land.RData", verbose = T)

grwlPR = grwl_by_wrstile %>% 
  as.data.frame() %>% 
  dplyr::select(PATH, ROW) %>% 
  as_tibble()

grwlWRS_sf = wrs_in_world %>% 
  inner_join(grwlPR, by = c("PATH", "ROW")) %>% 
  st_as_sf() %>%
  st_buffer(0.0)

world_clipped_nh = grwlWRS_sf %>% 
  st_union() %>% 
  st_intersection(extentsfc_nh)
world_clipped_sh = grwlWRS_sf %>% 
  st_union() %>% 
  st_intersection(extentsfc_sh)

plot(world_clipped_nh)
plot(world_clipped_sh)

st_write(world_clipped_nh, dsn = "~/Google_Drive/monthly_map_river_ice_future_change_global/world_clipped_grwl_nh.shp")
st_write(world_clipped_sh, dsn = "~/Google_Drive/monthly_map_river_ice_future_change_global/world_clipped_grwl_sh.shp")
```

### read the geotiff

```{r}
dat_dir = "~/Google_Drive/monthly_map_river_ice_future_change_global/"

filename = dir(dat_dir, pattern = "*.tif")
filepath = dir(dat_dir, pattern = "*.tif", full.names = T)
N = length(strsplit(filename, split = c("_")) %>% unlist())

fileinfo = tibble(filename, filepath) %>% 
  mutate(
    model = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 - 2 == 1],
    rcp = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 - 3 == 1],
    month = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 + 1 == 1]
    ) %>% 
  mutate(month = factor(month, levels = paste0(1:12, ".tif"), labels = month.abb)) %>% 
  mutate(
    cropped_file_name_nh = paste(model, rcp, month, "nh.tif", sep = "_"),
    cropped_file_name_sh = paste(model, rcp, month, "sh.tif", sep = "_"))

img_dat_monthly = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  
  require(raster)
  # https://joeyklee.github.io/broc-cli-geo/guide/XX_raster_cropping_and_clipping.html
  cmdstring = paste("cd ~/Google_Drive/monthly_map_river_ice_future_change_global/; gdalwarp -dstnodata -1 -cutline world_clipped_grwl_nh.shp -crop_to_cutline", fileinfo$filename[i], fileinfo$cropped_file_name_nh[i], sep = " ")
  
  system(cmdstring)
  
  cmdstring = paste("cd ~/Google_Drive/monthly_map_river_ice_future_change_global/; gdalwarp -dstnodata -1 -cutline world_clipped_grwl_sh.shp -crop_to_cutline", fileinfo$filename[i], fileinfo$cropped_file_name_sh[i], sep = " ")
  
  system(cmdstring)
  
  t1 = paste0(dat_dir, fileinfo$cropped_file_name_nh[i]) %>% 
    raster(band = 2) %>% 
    projectRaster(crs = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>% 
    as(Class = "SpatialPixelsDataFrame") %>% 
    as_tibble() %>% 
    mutate(rif = .[, 1] %>% pull) %>% 
    mutate(model = fileinfo$model[i],
           rcp = fileinfo$rcp[i],
           month = fileinfo$month[i]) %>% 
    dplyr::select(x, y, rif, model, rcp, month) %>% 
    mutate(hemisphere = "nh")
  
  t2 = paste0(dat_dir, fileinfo$cropped_file_name_sh[i]) %>% 
    raster(band = 2) %>% 
    projectRaster(crs = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>% 
    as(Class = "SpatialPixelsDataFrame") %>% 
    as_tibble() %>% 
    mutate(rif = .[, 1] %>% pull) %>% 
    mutate(model = fileinfo$model[i],
           rcp = fileinfo$rcp[i],
           month = fileinfo$month[i]) %>% 
    dplyr::select(x, y, rif, model, rcp, month) %>% 
    mutate(hemisphere = "sh")
  
  bind_rows(t1, t2)
  
}

save(img_dat_monthly, file = "outputs/predicted_img_dat_monthly_global.RData")

```

### export the figure

```{r}
load("outputs/predicted_img_dat_monthly_global.RData", verbose = T)

models = img_dat_monthly %>% dplyr::select(model) %>% distinct() %>% pull
rcps = img_dat_monthly %>% dplyr::select(rcp) %>% distinct() %>% pull
ns = img_dat_monthly %>% dplyr::select(hemisphere) %>% distinct() %>% pull

load("predicted_monthly_rif_change_NH.RData", verbose = T)
dat_monthly_rif_nh = dat_monthly_rif
load("predicted_monthly_rif_change_SH.RData", verbose = T)

dat_monthly_rif_nh = dat_monthly_rif_nh %>% 
  mutate(perc_change = diff21 / rif1)

dat_monthly_rif_sh = dat_monthly_rif_sh %>% 
  mutate(perc_change = diff21 / rif1)

xytext_nh = tibble(lon = 180, lat = 35) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_transform(3995) %>% 
  st_coordinates()

xytext_sh = tibble(lon = 0, lat = -35) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_transform(3031) %>% 
  st_coordinates()

dat_monthly_rif_nh = dat_monthly_rif_nh %>% 
  mutate(x = xytext_nh[1], y = xytext_nh[2])

dat_monthly_rif_sh = dat_monthly_rif_sh %>% 
  mutate(x = xytext_sh[1], y = xytext_sh[2])

require(scales)
for (j in 1:length(models)) {
  for (i in 1:length(rcps)) {
    for (k in 1:length(ns)) {
      
      if (ns[k] == "nh") {
        world_clipped = world_clipped_nh_land
        crs = 3995
        xlim = c(-7e06, 8e06)
        ylim = c(-6e06, 7e06)
        stats = dat_monthly_rif_nh %>% 
          filter(
            model == models[j],
            rcp == rcps[i])
        months_included = 1:12
        if (i == 1) {
          ncol = 3
          output_height = 200
          device = "tiff"
        } else if (i == 2) {
            ncol = 4
            output_height = 120
            device = "pdf"
        }
        output_width = 120
        
      } else if (ns[k] == "sh") {
        world_clipped = world_clipped_sh_land
        crs = 3031
        xlim = c(-7.5e6, 7.5e6)
        ylim = c(-7.5e6, 7.5e6)
        stats = dat_monthly_rif_sh %>% 
          filter(
            model == models[j],
            rcp == rcps[i])
        months_included = 6:8
        ncol = 3
        output_height = 55
        output_width = 120
        device = "tiff"
      }
      
      monthly_diffs = img_dat_monthly %>%   
        mutate(month_int = as.integer(month)) %>% 
        filter(
          model == models[j],
          rcp == rcps[i],
          hemisphere == ns[k]) %>% 
        filter(month_int %in% months_included) %>%
        mutate(month_txt = factor(month_int, levels = months_included, labels = month.abb[months_included], ordered = T)) %>% 
        ggplot() +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = "white", color = NA) +
        geom_tile(aes(x, y, fill = rif * 100)) +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = NA, color = "black", lwd = 0.2)  +
        coord_sf(crs = st_crs(crs), xlim = xlim, ylim = ylim) +
        facet_wrap(~month_txt, ncol = ncol) +
        theme(panel.grid.major = element_line(colour = "white"),
          axis.title = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          line = element_blank(),
          rect = element_blank(),
          panel.grid = element_blank(),
          legend.position = "bottom",
          text = element_text(family = "sans", size = 6))
      
      stats_fil = stats %>% 
        mutate(month_int = as.integer(month)) %>% 
        filter(month_int %in% months_included) %>% 
        mutate(month_txt = factor(month_int, levels = months_included, labels = month.abb[months_included], ordered = T))
      
      if (ns[k] == "nh") {
        monthly_diffs = monthly_diffs +
          geom_text(data = stats_fil, aes(x = x, y = y, label = paste0(round(diff21 * 100, digits = 1), " (", round(perc_change * 100, digits = 0), "%)")), size = 5 * 0.3528, color = "red", hjust = "middle") +
          scale_fill_gradientn(
            colors = c("red", "grey"),
            limits = c(-60, 0),
            guide = guide_colorbar(
              title = "River ice extent difference\n(percentage point)",
              title.theme = element_text(size = 6),
              label.theme = element_text(size = 6),
              direction = "horizontal",
              nbin = 20,
              draw.ulim = TRUE,
              draw.llim = TRUE,
              title.position = "top",
              title.hjust = 0.5,
              barheight = unit(2, units = "mm"),
              barwidth = unit(50, units = "mm")
              ),
            oob = squish)
      } else if(ns[k] == "sh") {
        monthly_diffs = monthly_diffs +
          # geom_text(data = stats_fil, aes(x = x, y = y, label = paste0(round(diff21 * 100, digits = 2), "\n(", round(perc_change * 100, digits = 0), "%)")), family = "sans", size = 5 * 0.3528, color = "red", hjust = "middle") +
          scale_fill_gradient2(
            low = "red",
            high = "blue",
            mid = "grey",
            midpoint = 0,
            limits = c(-15, 15),
            guide = guide_colorbar(
              title = "River ice extent difference\n(percentage point)",
              title.theme = element_text(size = 6),
              label.theme = element_text(size = 6),
              direction = "horizontal",
              nbin = 20,
              draw.ulim = TRUE,
              draw.llim = TRUE,
              title.position = "top",
              title.hjust = 0.5,
              barheight = unit(2, units = "mm"),
              barwidth = unit(50, units = "mm")
              ),
            oob = squish)
      }
      
      # monthly_diffs
      if (k == 1 & i == 2) {
        monthly_diffs %>%
          ggsave(filename = paste0("figs/final_figures/corrected_predicted_ice_extent_map_monthly_", models[j], "_", rcps[i], "_", ns[k], ".", device),
                 device = device,
                 units = "mm",
                 width = output_width,
                 height = output_height,
                 dpi = 300,
                 colormode = "rgb")
      } else {
        monthly_diffs %>%
          ggsave(filename = paste0("figs/final_figures/corrected_predicted_ice_extent_map_monthly_", models[j], "_", rcps[i], "_", ns[k], ".", device),
                 device = device,
                 units = "mm",
                 width = output_width,
                 height = output_height,
                 dpi = 300)
      }
      
      
      monthly_diffs %>%
        ggsave(filename = paste0("figs/final_figures/corrected_predicted_ice_extent_map_monthly_", models[j], "_", rcps[i], "_", ns[k], ".png"),
               device = "png",
               units = "mm",
               width = output_width,
               height = output_height,
               dpi = 300)
    }
    
  }
}

## output the data for figures

### fig3 this file is much bigger than the maximum allowed...
# dat_fig3 = img_dat_monthly %>% 
#   filter(hemisphere == "nh",
#          rcp == "rcp85") %>% 
#   select(x, y, rif, month)
# 
# dat_fig3 %>% write_csv(path = "figs/final_figures/fig3.csv")
```


## monthly predicted change stats NH

```{r}
dat_dir = "~/Google_Drive/predicted_monthly_river_ice_change_number"

filenames = dat_dir %>% dir(pattern = "monthlyStats_*")
filepaths = dat_dir %>% dir(pattern = "monthlyStats_*", full.names = T)
N = length(strsplit(filenames, split = c("_")) %>% unlist())

fileinfo = tibble(filenames, filepaths) %>% 
  mutate(
    model = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 1 == 1],
    rcp = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 2 == 1],
    file_index = as.integer(substr((strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 + 1 == 1], start = 1, stop = 1)))

ntropics = 3592377

w = tibble(
  file_index = 0:8, 
  npixels = c(450966, 1391761, 1142861, 139120, 1357906, 2508510, 4161942, 1074039, ntropics))

tropics = tibble(month = 1:12, period1 = 0, period2 = 0, diff = 0, file_index = 8)

w %>% summarise(npixels = sum(npixels))

i = 1
for (modeli in fileinfo %>% dplyr::select(model) %>% distinct() %>% pull) {
  for (rcpi in fileinfo %>% dplyr::select(rcp) %>% distinct() %>% pull) {
    temp = tropics %>% mutate(model = modeli, rcp = rcpi)
    if (i == 1) {
      tropicsAll = temp
      i = i + 1
    } else {
      tropicsAll = bind_rows(tropicsAll, temp)
    }
  }
}

tropicsAll = tropicsAll %>% 
  filter(!(model == "MIROC-ESM" & rcp == "rcp45"))

dat_monthly_rif = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  temp = fileinfo[i, ]
  read_csv(temp$filepaths) %>% 
    mutate(model = temp$model,
           rcp = temp$rcp,
           file_index = temp$file_index) %>% 
    dplyr::select(-`system:index`, -`.geo`)
} %>% 
  bind_rows(tropicsAll) %>% 
  left_join(w, by = "file_index") %>% 
  group_by(model, rcp, month) %>% 
  summarise(
    rif1 = weighted.mean(period1, w = npixels),
    rif2 = weighted.mean(period2, w = npixels),
    dff = rif2 - rif1,
    diff21 = weighted.mean(diff, w = npixels)) %>% 
  ungroup %>% 
  mutate(month = factor(month, levels = 1:12, labels = month.abb))

save(dat_monthly_rif, file = "predicted_monthly_rif_change_NH.RData")

dat_monthly_rif %>% 
  ggplot() +
  geom_bar(aes(x = month, y = dff * 100, fill = rcp), stat = "identity", position = "dodge") +
  facet_wrap(~model) +
  labs(
    x = "Month",
    y = "Changes in river ice extent (percentage points)",
    fill = "RCP"
  )

dat_monthly_rif %>% 
  filter(model == "CESM1-BGC") %>% 
  mutate(
    "Change (% point)" = format(diff21 * 100, digits = 1),
    "Change (%)" = paste0(format(diff21 / rif1 * 100, digits = 0), "%")) %>% 
  rename(
    Model = model,
    Scenario = rcp,
    Month = month) %>% 
  dplyr::select(-rif1, -rif2, -dff, -diff21) %>% 
  filter(Scenario == "rcp85")

dat_monthly_rif %>% 
  filter(model == "CESM1-BGC") %>% 
  filter(month %in% month.abb[-(6:8)]) %>%
  mutate(pchange = dff / rif1 * 100) %>% 
  group_by(rcp) %>% 
  summarise(dff_ns_mean = mean(dff),
            pchange_ns_mean = mean(pchange))
```

## monthly predicted change stats SH

```{r}
dat_dir = "~/Google_Drive/predicted_monthly_river_ice_change_number_SH"

filenames = dat_dir %>% dir(pattern = "monthlyStats_*")
filepaths = dat_dir %>% dir(pattern = "monthlyStats_*", full.names = T)
N = length(strsplit(filenames, split = c("_")) %>% unlist())

fileinfo = tibble(filenames, filepaths) %>% 
  mutate(
    model = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 1 == 1],
    rcp = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 2 == 1],
    file_index = as.integer(substr((strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 + 1 == 1], start = 1, stop = 1)))

w = tibble(
  file_index = 0:7, 
  npixels = c(0, 0, 3034688, 135893, 935265, 138134, 347300, 541229))

w %>% summarise(npixels = sum(npixels))

dat_monthly_rif_sh = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  temp = fileinfo[i, ]
  test = read_csv(temp$filepaths) %>% 
    mutate(model = temp$model,
           rcp = temp$rcp,
           file_index = temp$file_index)
  if (ncol(test) != 6) {
    test
  }
} %>% 
  dplyr::select(-`system:index`, -`.geo`) %>% 
  left_join(w, by = "file_index") %>% 
  group_by(model, rcp, month) %>% 
  summarise(
    rif1 = weighted.mean(period1, w = npixels),
    rif2 = weighted.mean(period2, w = npixels),
    dff = rif2 - rif1,
    diff21 = weighted.mean(diff, w = npixels)) %>% 
  ungroup %>% 
  mutate(month = factor(month, levels = 1:12, labels = month.abb))

save(dat_monthly_rif_sh, file = "predicted_monthly_rif_change_SH.RData")

load("predicted_monthly_rif_change_SH.RData", verbose = T)

dat_monthly_rif_sh %>% 
  ggplot() +
  geom_bar(aes(x = month, y = dff * 100, fill = rcp), stat = "identity", position = "dodge") +
  facet_wrap(~model) +
  labs(
    x = "Month",
    y = "Changes in river ice extent (percentage points)",
    fill = "RCP"
  )

dat_monthly_rif_sh %>% 
  filter(model == "CESM1-BGC") %>% 
  mutate(
    "Change (% point)" = format(diff21 * 100, digits = 1),
    "Change (%)" = paste0(format(diff21 / rif1 * 100, digits = 0), "%")) %>% 
  rename(
    Model = model,
    Scenario = rcp,
    Month = month) %>% 
  dplyr::select(-rif1, -rif2, -dff, -diff21) %>% 
  filter(Scenario == "rcp85")

dat_monthly_rif_sh %>% 
  filter(model == "CESM1-BGC") %>% 
  filter(month %in% month.abb[-(6:8)]) %>%
  mutate(pchange = dff / rif1 * 100) %>% 
  group_by(rcp) %>% 
  summarise(dff_ns_mean = mean(dff),
            pchange_ns_mean = mean(pchange))
```

## global monthly predicted change stats

```{r}
load("predicted_monthly_rif_change_NH.RData", verbose = T)
dat_monthly_rif_nh = dat_monthly_rif %>% 
  filter(model == "CESM1-BGC") %>% 
  mutate(w = 15819482)
load("predicted_monthly_rif_change_SH.RData", verbose = T)
dat_monthly_rif_sh = dat_monthly_rif_sh %>% 
  filter(model == "CESM1-BGC") %>% 
  mutate(w = 5132509)

dat_monthly_rif_global = dat_monthly_rif_nh %>% 
  bind_rows(dat_monthly_rif_sh) %>% 
  group_by(month, rcp) %>% 
  summarise(
    rif1 = weighted.mean(rif1, w = w),
    rif2 = weighted.mean(rif2, w = w),
    diff21 = weighted.mean(diff21, w = w)) %>% 
  ungroup %>% 
  mutate(pchange = diff21 / rif1 * 100)

dat_monthly_rif_global %>% 
  filter(rcp == "rcp85") %>% 
  summary()

dat_monthly_rif_global %>% 
  filter(rcp == "rcp45") %>% 
  summary()
```


## Map comparing durations

### read the geotiff

```{r}

dat_dir = "~/Google_Drive/duration_change_maps/"

filename = dir(dat_dir, pattern = "*.tif")
filepath = dir(dat_dir, pattern = "*.tif", full.names = T)
N = length(strsplit(filename, split = c("_")) %>% unlist())

fileinfo = tibble(filename, filepath) %>% 
  mutate(
    model = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 - 2 == 1],
    rcp = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 - 3 == 1],
    period = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 + 1 == 1]
    ) %>% 
  mutate(period = factor(period, levels = c("08-28.tif", "79-99.tif", "diff.tif"), labels = c("2009-2029", "2080-2100", "Difference"))) %>% 
  mutate(
    cropped_file_name_nh = paste(model, rcp, period, "nh.tif", sep = "_"),
    cropped_file_name_sh = paste(model, rcp, period, "sh.tif", sep = "_")
    )

img_dat = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  
  require(raster)
  # https://joeyklee.github.io/broc-cli-geo/guide/XX_raster_cropping_and_clipping.html
  cmdstring = paste("cd ~/Google_Drive/duration_change_maps/; gdalwarp -dstnodata -1 -cutline world_clipped_grwl_nh.shp -crop_to_cutline", fileinfo$filename[i], fileinfo$cropped_file_name_nh[i], sep = " ")
  
  system(cmdstring)
  
  cmdstring = paste("cd ~/Google_Drive/duration_change_maps/; gdalwarp -dstnodata -1 -cutline world_clipped_grwl_sh.shp -crop_to_cutline", fileinfo$filename[i], fileinfo$cropped_file_name_sh[i], sep = " ")
  
  system(cmdstring)
  
  t1 = paste0(dat_dir, fileinfo$cropped_file_name_nh[i]) %>% 
    raster() %>% 
    projectRaster(crs = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>% 
    as(Class = "SpatialPixelsDataFrame") %>% 
    as_tibble() %>% 
    mutate(duration = .[, 1] %>% pull) %>% 
    mutate(model = fileinfo$model[i],
           rcp = fileinfo$rcp[i],
           period = fileinfo$period[i]) %>% 
    dplyr::select(x, y, duration, model, rcp, period) %>% 
    mutate(hemisphere = "nh")
  
  t2 = paste0(dat_dir, fileinfo$cropped_file_name_sh[i]) %>% 
    raster() %>% 
    projectRaster(crs = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>% 
    as(Class = "SpatialPixelsDataFrame") %>% 
    as_tibble() %>% 
    mutate(duration = .[, 1] %>% pull) %>% 
    mutate(model = fileinfo$model[i],
           rcp = fileinfo$rcp[i],
           period = fileinfo$period[i]) %>% 
    dplyr::select(x, y, duration, model, rcp, period) %>% 
    mutate(hemisphere = "sh")
  
  bind_rows(t1, t2)
}

save(img_dat, file = "outputs/duration_map_global.RData")

```

### export the figure

```{r}
load("outputs/duration_map_global.RData", verbose = T)

img_dat_comp = img_dat %>% 
  filter(period != "Difference") %>% 
  mutate(duration_grp = base::cut(duration, breaks = c(0, 5, 15, 90, 180, 366), labels = c("Ice free", "Intermittent", "0.5-3 months", "3-6 months", ">6 months"), include.lowest = T))

img_dat_comp = img_dat_comp %>% 
  mutate(period = factor(as.character(period), levels = c("2008-2028", "2079-2099"), labels = c("2009-2029", "2080-2100")))

models = img_dat %>% dplyr::select(model) %>% distinct() %>% pull
rcps = img_dat %>% dplyr::select(rcp) %>% distinct() %>% pull
ns = img_dat %>% dplyr::select(hemisphere) %>% distinct() %>% pull

# world_clipped_unioned = world_clipped %>% st_transform(crs = crs) %>% st_union()

for (j in 1:length(models)) {
  for (i in 1:length(rcps)) {
    for (k in 1:length(ns)) {
      
      if (ns[k] == "nh") {
        world_clipped = world_clipped_nh_land %>% st_union
        crs = 3995
        xlim = c(-7e06, 8e06)
        ylim = c(-6e06, 7e06)
        if (i == 1) {
          device = "tiff"
        } else if (i == 2) {
          device = "pdf"
        }
      } else if (ns[k] == "sh") {
        world_clipped = world_clipped_sh_land %>% st_union
        crs = 3031
        xlim = c(-7.5e6, 7.5e6)
        ylim = c(-7.5e6, 7.5e6)
        device = "tiff"
      }
      
      duration_comps = img_dat_comp %>% 
        filter(
          model == models[j],
          rcp == rcps[i],
          hemisphere == ns[k]) %>% 
        ggplot() +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = "white", color = NA) +
        geom_tile(aes(x, y, fill = duration_grp), color = NA) +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = NA, color = "black", lwd = 0.3) +
        scale_fill_viridis_d(
          begin = 0.1,
          guide = guide_legend(
            direction = "horizontal",
            title.position = "top",
            title.hjust = 0.5,
            keywidth = 2.5,
            keyheight = 0.5,
            label.position = "bottom",
            label.hjust = 0.5
          )) +
        coord_sf(crs = crs, x = xlim, y = ylim) +
        labs(fill = "River ice duration") +
        theme(panel.grid.major = element_line(color = "white", size = 0.5),
              axis.text.y = element_blank(),
              axis.text.x = element_blank(),
              line = element_blank(),
              rect = element_blank(),
              text = element_text(size = 6),
              panel.grid = element_blank(),
              legend.position = "bottom",
              axis.title = element_blank()) +
        facet_wrap(~period)
      
      # duration_comps
      
      duration_comps %>% 
        ggsave(filename = paste0("figs/final_figures/corrected_ice_duration_map_comparison_", models[j], "_", rcps[i], "_", ns[k], ".", device),
               device = device,
               units = "mm",
               width = 120,
               height = 70)
      
      duration_comps %>% 
        ggsave(filename = paste0("figs/final_figures/corrected_ice_duration_map_comparison_", models[j], "_", rcps[i], "_", ns[k], ".png"),
               # device = device,
               units = "in",
               width = 6,
               height = 3.5)
    }
    
  }
}

```

## Map of changes in duration

```{r}
img_dat_diff = img_dat %>% 
  filter(period == "Difference")

models = img_dat %>% dplyr::select(model) %>% distinct() %>% pull
rcps = img_dat %>% dplyr::select(rcp) %>% distinct() %>% pull
ns = img_dat %>% dplyr::select(hemisphere) %>% distinct() %>% pull

require(scales)
for (j in 1:length(models)) {
  for (i in 1:length(rcps)) {
    for (k in 1:length(ns)) {
      
      if (ns[k] == "nh") {
        world_clipped = world_clipped_nh_land %>% st_union
        crs = 3995
        xlim = c(-7e06, 8e06)
        ylim = c(-6e06, 7e06)
      } else if (ns[k] == "sh") {
        world_clipped = world_clipped_sh_land %>% st_union
        crs = 3031
        xlim = c(-7.5e6, 7.5e6)
        ylim = c(-7.5e6, 7.5e6)
      }
      
      duration_diff = img_dat_diff %>% 
        filter(
          rcp == rcps[i],
          model == models[j],
          hemisphere == ns[k]) %>% 
        ggplot() +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = "white", color = NA) +
        geom_tile(aes(x, y, fill = duration)) +
        geom_sf(data = world_clipped %>% st_transform(crs = crs), fill = NA, color = "black", lwd = 0.3) +
        scale_fill_gradient2(
          low = "red", 
          mid = "grey", 
          high = "blue", 
          midpoint = 0,
          limits = c(-90, 5),
          # breaks = c(-150, -90, -60, -30, -20, -10),
          # trans = scales::pseudo_log_trans(base = 2, sigma = 3),
          guide = guide_colorbar(
            direction = "horizontal",
            title.position = "top",
            title.hjust = 0,
            label.position = "bottom",
            label.hjust = 0.5,
            barheight = unit(1.5, units = "mm"),
            barwidth = unit(35, units = "mm")
          ),
          oob = squish) +
        coord_sf(crs = crs, x = xlim, y = ylim) +
        labs(fill = "Changes in river ice duration (days)\nbetween 2009-2029 and 2080-2100") +
        theme(panel.grid.major = element_line(color = "white", size = 0.5),
              axis.text.y = element_blank(),
              axis.text.x = element_blank(),
              line = element_blank(),
              rect = element_blank(),
              text = element_text(size = 6),
              panel.grid = element_blank(),
              legend.position = "bottom",
              axis.title = element_blank())
      
      # duration_diff
      
      duration_diff %>% 
        ggsave(filename = paste0("figs/final_figures/corrected_ice_duration_map_difference_", models[j], "_", rcps[i], "_", ns[k], ".pdf"),
               device = "pdf",
               units = "mm",
               width = 60,
               height = 70,
               colormode = "rgb")
    }
    
   }
}

```


```{r}
img_dat_comp = img_dat %>% 
  filter(period != "Difference") %>% 
  mutate(duration_grp = base::cut(duration, breaks = c(0, 5, 15, 90, 180, 366), labels = c("Ice free", "Intermittent", "0.5-3 months", "3-6 months", ">6 months"), include.lowest = T))
dat_dir = "~/Google_Drive/predicted_duration_change"

files = dir(dat_dir, full.names = T)

temp = files[8]

require(raster)

img = raster(temp) %>%
  projectRaster(crs = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") # reproject to EPSG:3995

plot(img)

img_tb = as_tibble(as(img, "SpatialPixelsDataFrame")) %>% 
  mutate(duration_grp = base::cut(test_clipped_output_4, breaks = c(0, 5, 15, 90, 180, 366), labels = c("Ice free", "Intermittent", "0.5-3 months", "3-6 months", ">6 months"), include.lowest = T))

img_df %>% 
ggplot() +
  geom_sf(data = world_clipped %>% st_transform(crs = 3995), fill = "darkgrey", color = NA) +
  geom_tile(aes(x, y, fill = duration_grp)) +
  geom_sf(data = world_clipped %>% st_transform(crs = 3995), fill = NA, color = "black") +
  scale_fill_viridis_d(
    begin = 0.1,
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      title.hjust = 0.5,
      keywidth = 3,
      keyheight = 1,
      label.position = "bottom",
      label.hjust = 0.5
      )) +
  coord_sf(crs = 3995, x = c(-7e06, 8e06), y = c(-6e06, 5e06)) +
  labs(fill = "River ice duration") +
  theme(panel.grid.major = element_line(color = "grey", size = 0.5),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        line = element_blank(),
        rect = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title = element_blank())
```


## duration difference between 0828 and 7999 stats

### global duration difference

```{r}
load("outputs/annual_dat_global_0899.RData", verbose = T)

duration_two_periods = annual_dat_global %>% 
  filter((year >= 2009 & year < 2029) | (year >= 2080 & year < 2100)) %>% 
  mutate(period = base::cut(year, breaks = c(2008, 2050, 2099), labels = c("y0929", "y8000"), include.lowest = T)) %>% 
  group_by(model, rcp, period) %>% 
  summarise(
    mean_ice_duration = mean(mean_ice_duration),
    mean_SAT = mean(gsat)
  ) %>% 
  ungroup()

## duration change
duration_two_periods %>% 
  dplyr::select(-mean_SAT) %>% 
  spread(period, mean_ice_duration) %>% 
  mutate(duration_diff = y8000 - y0929) %>% 
  dplyr::select(-y8000, -y0929)

## SAT change
duration_two_periods %>% 
  dplyr::select(-mean_ice_duration) %>% 
  spread(period, mean_SAT) %>% 
  mutate(meanSAT_diff = y8000 - y0929) %>% 
  dplyr::select(-y8000, -y0929)
```

### NH duration difference

```{r}
dat_dir = "~/Google_Drive/predicted_monthly_river_ice_change_number"

filenames = dat_dir %>% dir(pattern = "monthlyDuration_*")
filepaths = dat_dir %>% dir(pattern = "monthlyDuration_*", full.names = T)
N = length(strsplit(filenames, split = c("_")) %>% unlist())

fileinfo = tibble(filenames, filepaths) %>% 
  mutate(
    model = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 1 == 1],
    rcp = (strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 - 2 == 1],
    file_index = as.integer(substr((strsplit(filenames, split = c("_")) %>% unlist())[(1:N) %% 4 + 1 == 1], start = 1, stop = 1)))

w = tibble(
  file_index = 0:7, 
  npixels = c(450966, 1391761, 1142861, 139120, 1357906, 2508510, 4161942, 1074039))

dat_duration = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  temp = fileinfo[i, ]
  read_csv(temp$filepaths) %>% 
    mutate(model = temp$model,
           rcp = temp$rcp,
           file_index = temp$file_index)
} %>% 
  dplyr::select(-`system:index`, -`.geo`) %>% 
  left_join(w, by = "file_index") %>% 
  group_by(model, rcp) %>% 
  summarise(
    d1 = weighted.mean(duration1, w = npixels),
    d2 = weighted.mean(duration2, w = npixels),
    dff = d2 - d1,
    diff_duration = weighted.mean(diff_duration, w = npixels)) %>% 
  ungroup

dat_duration %>% 
  ggplot() +
  geom_bar(aes(x = model, y = dff, fill = rcp), stat = "identity", position = "dodge") +
  # facet_wrap(~model) +
  labs(
    x = "Model",
    y = "Changes in NH mean river ice duration",
    fill = "RCP"
  )

dat_duration %>% 
  # filter(model == "CESM1-BGC") %>% 
  dplyr::select(model, rcp, dff) %>% 
  mutate(rate = dff / (2089 - 2018) * 100,
         change_by_99 = rate * 90 / 100)
```


## change of % of ice-affected river due to increasing global mean SAT

### read in annual global mean temperature

```{r}
dat_dir = "~/Google_Drive/predicted_global_mean_SAT/"
filename = dir(dat_dir)
filepath = dir(dat_dir, full.names = T)
N = length(strsplit(filename, split = c("_")) %>% unlist())

fileinfo = tibble(filename, filepath) %>% 
  mutate(
    model = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 == 1],
    rcp = (strsplit(filename, split = c("_")) %>% unlist())[(1:N) %% 5 - 1 == 1])

dat_temperature = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  temp = fileinfo[i, ]
  read_csv(temp$filepath) %>% 
    mutate(model = temp$model,
           rcp = temp$rcp)
} %>% 
  dplyr::select(-`system:index`, -`.geo`, -tasmean_count) %>% 
  rename(gsat = tasmean_mean)

dat_temperature = dat_temperature %>% 
  group_by(model, rcp) %>% 
  do({
    dat = .
    base_temp = dat %>% 
      filter(year >= 2009, year < 2029) %>% 
      summarise(base_gsat0929 = mean(gsat)) %>% 
      pull(base_gsat0929)
    dat %>% 
      mutate(delta_gsat = gsat - base_temp)
  }) %>% 
  ungroup()

dat_temperature %>% 
  ggplot() +
  geom_point(aes(x = year, y = delta_gsat, color = rcp, pch = model))

dat_temperature %>% 
  ggplot() +
  geom_point(aes(x = delta_gsat, y = gsat, color = rcp, pch = model))

save(dat_temperature, file = "outputs/annual_global_mean_temperature_0899.RData")

dat_temperature %>% 
  filter(year <= 2028 | year >= 2080) %>% 
  mutate(period = factor(year >= 2050, levels = c(F, T), labels = c("2009-2029", "2079-2099"))) %>% 
  group_by(model, rcp) %>% 
  do({
    dat = .
    tibble(diff_gsat = mean(dat %>% filter(period == "2079-2099") %>% pull(gsat)) - 
             mean(dat %>% filter(period == "2009-2029") %>% pull(gsat)))
  }) %>% 
  ungroup()

```

### changing ice duration and % of ice-affected rivers v.s. global mean temperature change

```{r}
w = read_csv("data/grwlByLonBandsNpoints.csv")
w = w %>% 
  arrange(lonStart) %>% 
  mutate(file_index = 0:(nrow(.) - 1)) %>% 
  dplyr::select(-lonStart)

dat_dir = "~/Google_Drive/predicted_annual_river_ice_global_0899"

files = dir(dat_dir, pattern = "*.csv")
filepaths = dir(dat_dir, pattern = "*.csv", full.names = T)
N = length(strsplit(files, split = c("_")) %>% unlist())

fileinfo = tibble(files, filepaths) %>% 
  mutate(
    model = (strsplit(files, split = c("_")) %>% unlist())[(1:N) %% 3 == 1],
    rcp = (strsplit(files, split = c("_")) %>% unlist())[(1:N) %% 3 - 1 == 1],
    file_index = as.integer(strsplit((strsplit(files, split = c("_")) %>% unlist())[(1:N) %% 3 + 1 == 1], split = c(".csv")) %>% unlist()))

dat = foreach(i = 1:nrow(fileinfo), .combine = "bind_rows") %do% {
  current_file = fileinfo[i, ]
  read_csv(current_file$filepaths) %>% 
    mutate(model = current_file$model,
           rcp = current_file$rcp,
           file_index = current_file$file_index)
} %>% 
  left_join(w, by = "file_index") %>% 
  dplyr::select(-`system:index`, -`.geo`)


annual_dat_global = dat %>% 
  group_by(model, rcp, year) %>% 
  summarise(
    ice_affected = weighted.mean(iceAffected, w = npoints),
    mean_ice_duration = weighted.mean(iceStatus, w = npoints) * 365,
    mean_ice_fraction = weighted.mean(iceFraction, w = npoints),
    mean_tas = weighted.mean(tasmean, w = npoints)) %>% 
  ungroup()

annual_dat_global %>% 
  # filter(model == "CESM1-BGC") %>% 
  ggplot() +
  geom_point(aes(x = mean_tas, y = ice_affected, color = rcp)) +
  geom_smooth(aes(x = mean_tas, y = ice_affected, color = rcp), method = lm) +
  facet_wrap(~model)

## attach gsat
load("outputs/annual_global_mean_temperature_0899.RData", verbose = T)
annual_dat_global = annual_dat_global %>% 
  left_join(dat_temperature, by = c("model", "rcp", "year"))

save(annual_dat_global, file = "outputs/annual_dat_global_0899.RData")

load("outputs/annual_dat_global_0899.RData", verbose = T)

predicted_ice_affected_percent_global = annual_dat_global %>% 
  # filter(model == "CESM1-BGC") %>% 
  ggplot() +
  geom_point(aes(x = year, y = ice_affected, color = rcp)) +
  geom_smooth(aes(x = year, y = ice_affected, color = rcp), method = lm) +
  facet_wrap(~model) +
  labs(
    x = "Year",
    y = "Percent of river affected by ice",
    color = "RCP"
  ) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.5, linetype = 1),
      panel.grid.major.x = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(colour = "black"),
      strip.text = element_text(size = 12, face = "bold"),
      strip.background = element_blank())

predicted_ice_affected_percent_global

predicted_ice_affected_percent_global %>% 
  ggsave(
    filename = "figs/predicted_ice_affected_percent_global.png",
    width = 6,
    height = 5
  )

h = 0.5
w = 0.5

predicted_ice_affected_percent_global_temp = annual_dat_global %>% 
  ggplot() +
  geom_point(aes(x = delta_gsat, y = ice_affected, pch = model, color = rcp), size = 0.7) +
  geom_smooth(aes(x = delta_gsat, y = ice_affected), method = lm, lwd = 0.5, se = F) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), breaks = seq(0.39, 53, by = 0.04)) +
  labs(
    x = "Changes in global mean SAT (ºC)\n compared to 2009-2029",
    y = "Percent of rivers affected by ice",
    color = "RCP",
    pch = "Model"
  ) +
  scale_shape_discrete(guide = guide_legend(title.position = "top", direction = "vertical", keyheight = h, keywidth = w)) +
  scale_color_discrete(guide = guide_legend(title.position = "top", direction = "vertical", keyheight = h, keywidth = w)) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.5, linetype = 1),
      panel.grid.major.x = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      text = element_text(size = 6),
      axis.line = element_line(colour = "black"),
      strip.text = element_text(size = 6),
      strip.background = element_blank(),
      legend.position = "bottom",
      legend.background = element_rect(fill = "white"))

predicted_ice_affected_percent_global_temp

predicted_ice_affected_percent_global_temp %>% 
  ggsave(
    filename = "figs/final_figures/predicted_ice_affected_percent_global_delta_gsat.tiff",
    units = "mm",
    width = 120,
    height = 140,
    dpi = 300
  )



predicted_ice_duration_global_temp = annual_dat_global %>% 
  ggplot() +
  geom_point(aes(x = delta_gsat, y = mean_ice_duration, pch = model, color = rcp, fill = rcp), size = 0.7) +
  geom_smooth(aes(x = delta_gsat, y = mean_ice_duration), method = lm, lwd = 0.5, se = F) +
  labs(
    x = "Changes in global mean SAT (ºC)\n compared to 2009-2029",
    y = "Global mean ice duration (days)",
    color = "RCP",
    pch = "Model"
  ) +
    # scale_shape_manual(values = c("CESM1-BGC" = 15, "GFDL-ESM2M" = 17, "MIROC-ESM" = 21), guide = guide_legend(title.position = "top", direction = "vertical", keyheight = h, keywidth = w)) +
  scale_shape_discrete(guide = guide_legend(title.position = "top", direction = "vertical", keyheight = h, keywidth = w)) +
  scale_color_discrete(guide = guide_legend(title.position = "top", direction = "vertical", keyheight = h, keywidth = w)) +
  guides(fill = F) +
  theme(panel.grid.major.y = element_line(color = "lightgrey", size = 0.5, linetype = 1),
      panel.grid.major.x = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(colour = "black"),
      text = element_text(size = 5),
      strip.background = element_blank(),
      legend.position = "bottom",
      legend.background = element_rect(fill = "white"))

predicted_ice_duration_global_temp

predicted_ice_duration_global_temp %>% 
  ggsave(
    filename = "figs/final_figures/corrected_predicted_ice_duration_global_delta_gsat.pdf",
    units = "mm",
    dpi = 300,
    colormode = "rgb",
    width = 60,
    height = 70,
    useDingbats = FALSE # https://stackoverflow.com/questions/9992275/ggplot2-pdf-import-in-adobe-illustrator-missing-font-adobepistd
  )


summary(lm(mean_ice_duration ~ delta_gsat, data = annual_dat_global))

summary(lm(ice_affected ~ delta_gsat, data = annual_dat_global))

```

```{r}
annual_dat_global %>% 
  ggplot() +
  geom_point(aes(x = ice_affected, y = mean_ice_duration, pch = model, color = rcp)) +
  geom_smooth(aes(x = ice_affected, y = mean_ice_duration), method = lm) +
  labs(
    x = "% of ice-affected rivers",
    y = "Ice duration (day)",
    color = "RCP",
    pch = "Model"
  )
```