data-raw/get-gamma-from-isst.R

In EarthSystemDiagnostics/spectraluncertainty:

# Estimate power at 1/2 years for all hadisst locations
library(tidyverse)
library(lubridate)
#library(ncdf4)

load(file = system.file("extdata/isst.ann.RData", package = "ecusdata"))

# fit spectra and estimate gamma

sub <- isst.ann %>%
  filter(lat == isst.lat.lon$lat[5001], lon == isst.lat.lon$lon[5001])

library(PaleoSpec)

GetSpec <- function(vec){
  ts1 <- ts(vec)
  spec <- SpecMTM(ts1)
  tb <- tibble(freq = spec$freq, spec = spec$spec)
  # remove lowest 2 frequencies
  tb <- tb[3:nrow(tb),]
  return(tb)
}

sp <- GetSpec(sub$sst.mean)

plot(log(spec) ~ log(freq), data = sp)
abline(lm(log(spec) ~ log(freq), data = sp))
abline(v = log(1/2))

# table(isst.ann$n)
#
# dY <- isst.ann %>%
#   #filter(n == 12) %>%
#   arrange(year) %>%
#   group_by(lat, lon) %>%
#   mutate(dY = c(NA, diff(year))) #%>%
#   #summarise(all.1 = all(dY == 1, na.rm = T))

isst.spec <- isst.ann %>%
  #filter(n == 12) %>%
  group_by(lat, lon) %>%
  do({
    GetSpec(.$sst.mean)
  })

save(isst.spec, file = "inst/extdata/isst/isst.spec.RData")
#load(file = "inst/extdata/isst/isst.spec.RData")


# # fit lm and get power at 1/2

# sub <- isst.spec %>%
#   filter(lat == 1.5, lon == -91.5) %>%
#   mutate(lon.lat = paste(lon, lat))
#
# sub %>%
#   ggplot(aes(x = freq, y = spec))+
#   geom_line()+
#   facet_grid(~lon.lat) +
#   scale_x_continuous(trans = "log10") +
#   scale_y_continuous(trans = "log10")


# plot(log(spec)~log(freq), data = sub)
# lm1 <- lm(log(spec)~log(freq), data = sub)
# abline(lm1)
# cfs <- coef(lm1)
# abline(h = cfs[1] + cfs[2] * log(0.5))
# abline(v = log(0.5))
#
# plot((spec)~(freq), data = sub)
# abline(h = exp(cfs[1] + cfs[2] * log(0.5)))
# abline(v = (0.5))


isst.pow <- isst.spec %>%
  do({
    lm1 <- lm(log(.$spec) ~ log(.$freq))
    cfs <- coef(lm1)
    tibble(slope = cfs[2], intercept = cfs[1])
  })

isst.pow <- isst.pow %>%
  mutate(pow.0.5 = exp(intercept + slope * log(0.5)),
         pow.0.1 = exp(intercept + slope * log(1/10)),
         pow.0.05 = exp(intercept + slope * log(1/20)),
         pow.0.02 = exp(intercept + slope * log(1/50)))

save(isst.pow, file = "inst/extdata/isst/isst.pow.RData")
#load("inst/extdata/isst/isst.pow.RData")

hist(isst.pow$pow.0.5, 100)
hist(isst.pow$slope, 100)