In AustralianAntarcticDivision/aceecostats: Ecosystem status and trends

knitr::opts_chunk$set(echo = TRUE)
library(raadtools)
library(tibble)
library(aceecostats)
library(dplyr)
library(tabularaster)
library(ggplot2)
library(viridis)

Ice prepare

library(raadtools)
library(tibble)
library(aceecostats)
library(dplyr)
library(tabularaster)
library(ggplot2)
library(virids)
ncf <- icefiles(hemi = "north")
scf <- icefiles(hemi = "south")



smap <- spTransform(aes_region, projection(readice()))

north <- readice(latest = TRUE, hemi = "north", inputfiles = ncf)
south <- readice(latest = TRUE, hemi = "south", inputfiles = scf)


cn <- cellnumbers(south[[1]], smap)
south_index <- tibble(cell_ = seq(ncell(south))) %>% 
  filter(!is.na(values(south[[1]]))) %>%  
  left_join(cn) %>% filter(!is.na(object_))

read_south <- function(date) {
  south_index %>% mutate(iceconc = extract(readice(date, inputfiles = scf, setNA = FALSE), cell_)[,1]) %>% 
    group_by(object_) %>% summarize(date = date, icep = mean(iceconc), notzero = sum(iceconc > 0))
}


north_index <- tibble(cell_ = seq(ncell(north))) %>% 
  filter(!is.na(values(north[[1]]))) 

read_north <- function(date) {
  obj <- values(readice(date, inputfiles = ncf, hemi= "north",  setNA = FALSE)[[1]])
  north_index %>% 
    summarize(date = date, iceconc = mean(obj[cell_]), notzero = sum(obj>0))
}

system.time({
north_df <- bind_rows(lapply(ncf$date, read_north))
})
## 800s
system.time({
south_df <- bind_rows(lapply(scf$date, read_south))
})
## 3600s
#save.image("ice_image.rdata1")
save(north_df, south_df,scf, ncf, north, south, file = "ice_image.rdata")

load("ice_image.rdata")
south_ <- south_df %>% 
  inner_join(aes_region@data %>% transmute(object_ = as.character(row_number()), SectorName = SectorName, Name = Name) %>%  mutate(object_ = as.character(row_number())))

asint <- function(x, tok) as.integer(format(x, tok))
yrs <- as.character(1978:2016)
colpal <- setNames(c(viridis::viridis(length(yrs)-1), "firebrick"), yrs)

north_df$date <- ncf$date
ggplot(north_df %>% 
         mutate(year = asint(date, "%Y"), doy = asint(date, "%j"), month = asint(date, "%m"))) +
  aes(x = doy, y = iceconc, group = as.character(year), color = as.character(year)) + 
  geom_line() + 
  #facet_wrap(~ob)
  scale_color_manual(values = colpal)


ggplot(south_ %>% ##filter(SectorName == "Atlantic") %>% 
         mutate(year = asint(date, "%Y"), doy = asint(date, "%j"), month = asint(date, "%m"))) +
  aes(x = doy, y = icep, group = year, color = as.character(year)) + 
  geom_line() + 
  facet_wrap(~Name) + 
  scale_color_manual(values = colpal) + 
    guides(color = FALSE) 


globe <- bind_cols(
  south_ %>% #mutate(doy  = asint(date, "%j")) %>% 
  group_by(date) %>% summarize(icepS = mean(icep)) %>% 
  ungroup(), 
  north_df %>% transmute(icepN = iceconc)
)

n_n <- sum(!is.na(values(north)))
s_n <- sum(!is.na(values(south)))
ggplot(globe  %>% 
         mutate(year = asint(date, "%Y"), doy = asint(date, "%j"))) + 
  aes(x = doy, y = icepS*s_n + icepN*n_n, group = year, color = as.character(year)) + geom_line() +  
  scale_color_manual(values = colpal)