# run using R CMD BATCH significance.r 
# outputs to Rplots.pdf and significance.r.Rout
# load custom functions
source("functions.r")

# for generalised least squares
library(nlme)

# options
xunits="year"
textsize=1.4
titlesize=1.8
colfit="red"
pch1=20#points
decimate=1#datasets with >10,000 pts default to 100 for speed
stepsize=5#spacing between starting times of trend/accel fits
# for piecewise
col1=4#blue
col2=2#red

# read basin data
if(2==1){
  indata = read.table("greenland2013/GIS_climate.nasa.txt",header=T)
  title="Greenland mass"
  yunits="gigatons"
  tlims=c(-350,-190)
  alims=c(-60,0)
  indata = indata[which(indata$x>2008.0),]
}
if(2==1){
  indata = read.table("greenland.txt",header=T)
  title="Greenland mass"
  yunits="gigatons"
  tlims=c(-350,-190)
  alims=c(-60,0)
  #indata = indata[which(indata$x>2003.87),]
  #indata = indata[which(indata$x<2009.5),]
}
if(2==1){
  indata = read.table("Greenland_timeseries_SS23x_FIX.txt",header=T)
  title="Greenland mass, SS23x_FIX"
  yunits="gigatons"
  tlims=c(-350,-190)
  alims=c(-60,0)
  #indata = indata[which(indata$x>2003.87),]
  #indata = indata[which(indata$x<2009.5),]
}
if(2==1){
  indata = read.table("antarctica.txt",header=T)
  title="Antarctic mass"
  yunits="gigatons"
  tlims=c(-200,200)
  alims=c(-100,100)
}
if(2==1){
  indata = read.table("wais.txt",header=T)
  title="W. Antarctic mass"
  yunits="gigatons"
  tlims=c(-220,0)
  alims=c(-30,30)
}
if(2==1){
  indata = read.table("eais.txt",header=T)
  title="E. Antarctic mass"
  yunits="gigatons"
  tlims=c(0,300)
  alims=c(-100,100)
}
if(2==1){
  indata = read.table("grnlnd1.txt",header=T)
  title="Greenland01 mass"
  yunits="gigatons"
  tlims=c(-100,0)
  alims=c(-10,20)
}
if(2==1){
  indata = read.table("grnlnd2.txt",header=T)
  title="Greenland02 mass"
  yunits="gigatons"
  tlims=c(-10,0)
  alims=c(-10,15)
}
if(2==1){
  indata = read.table("grnlnd3.txt",header=T)
  title="Greenland03 mass"
  yunits="gigatons"
  tlims=c(-20,-5)
  alims=c(-5,15)
}
if(2==1){
  indata = read.table("grnlnd4.txt",header=T)
  title="Greenland04 mass"
  yunits="gigatons"
  tlims=c(-50,-20)
  alims=c(-10,20)
}
if(2==1){
  indata = read.table("grnlnd5.txt",header=T)
  title="Greenland05 mass"
  yunits="gigatons"
  tlims=c(-100,-40)
  alims=c(-10,5)
}
if(2==1){
  indata = read.table("grnlnd6.txt",header=T)
  title="Greenland06 mass"
  yunits="gigatons"
  tlims=c(-25,-15)
  alims=c(-3,3)
}

# or extract a climate data series
if(2==1){
  indata = read.table("rawdata.txt",header=T)

  if(2==1){
    keep=c("t","giss")
    title="GISS"
    yunits="°C"
    tlims=c(-.02,.04)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","ncdc")
    title="NCDC"
    yunits="°C"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","cru")
    title="CRU"
    yunits="°C"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","rss")
    title="RSS"
    yunits="°C"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","uah")
    title="UAH"
    yunits="°C"
    tlims=c(-.02,.04)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","mei")
    title="Multivariate el Nino index"
    yunits="#"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }
  if(2==1){
    keep=c("t","aod")
    title="Aerosol Optical Depth"
    yunits="#"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }
  if(2==1){
    title="Sunspots"
    keep=c("t","ssn")
    yunits="number"
    tlims=c(-.01,.03)
    alims=c(-.01,.01)
  }

#THINK ABOUT THIS tau = t-1990             # time (zeroed to 1990)

  indata = indata[keep]
  colnames(indata) <- c("x", "y")
  indata = indata[which(!is.na(indata$y)),]

  #To reproduce Tamino's "How Long" analysis.
  #indata = indata[which(indata$x>1975),]
  #indata = indata[which(indata$x<2009.8),]
  
}

# or Church and White reconstructed sea level
if(1==1){
  indata = read.table("Church_White_2011/CSIRO_Recons_gmsl_yr_2011.txt",header=T)
  keep=c("x","y")
  title="GMSL"
  yunits="mm"
  tlims=c(1.4,3.2)
  alims=c(0.0,.1)
  indata = indata[keep]
}

# or Church and White altimetry sea level
if(2==1){
  indata = read.table("Church_White_2011/CSIRO_ALT_gmsl_yr_2011.txt",header=T)
  keep=c("x","y")
  title="GMSL"
  yunits="mm"
  tlims=c(1.4,4.0)
  alims=c(-.1,.25)
  indata = indata[keep]
}

# or Cryosphere Today NH sea ice area
if(2==1){
  indata = read.table("cryosphere-today-NH-area-anomaly-1979-2008.txt",header=T)
  keep=c("x","y")
  title="NH sea ice area anom."
  yunits="10^6 km^2"
  tlims=c(-.3,.3)
  alims=c(-.03,.03)
  indata = indata[keep]
}

# or Cryosphere Today SH sea ice area
if(2==1){
  indata = read.table("cryosphere-today-SH-area-anomaly-1979-2008.txt",header=T)
  keep=c("x","y")
  title="SH sea ice area anom."
  yunits="10^6 km^2"
  tlims=c(-.3,.3)
  alims=c(-.03,.03)
  indata = indata[keep]
}

# or NSIDC NH sea ice extent
if(2==1){
  indata = read.csv("NSIDC/NH_seaice_extent_final.csv")
  keep=c("x","y")
  title="NH sea ice extent."
  colnames(indata)[4] <- "y"
  yunits="10^6 km^2"
  tlims=c(-.3,.3)
  alims=c(-.03,.03)
  indata = indata[keep]
}

# or new RSS data
if(2==1){
  indata = read.table("new_RSS.txt",header=T)
  title="RSS"
  yunits="°C"
  tlims=c(-.02,.04)
  alims=c(-.01,.01)
}

# or new UAH data
if(2==1){
  indata = read.table("new_UAH.txt",header=T)
  title="UAH"
  yunits="°C"
  tlims=c(-.02,.04)
  alims=c(-.01,.01)
}

# if t/alims have length 1, the many_trends/accels plots autosize limits
#tlims=-1
#alims=-1

# Calculate length "n" of timeseries, in case we decimate data for speed.
n = length(indata$x)

# large datasets should be decimated automatically.
#if(n > 10000 && decimate < 100) decimate = 100
n

indata = indata[seq(1,n,by=decimate),]
# remove mean again
indata$y = indata$y - mean(indata$y)

# need to recalculate length if decimated
n = length(indata$x)
n

midpoint=(indata$x[n]+indata$x[1])/2.0

# set center of timeseries to zero
#indata$x = indata$x - midpoint

plot_data_and_fit(indata,fit,n,xunits,yunits,title,colfit,pch1,titlesize,textsize)
fit_many_trends(indata,n,stepsize,tlims,xunits,yunits,title,colfit,pch1,titlesize,textsize)
fit_many_accels(indata,n,stepsize,alims,xunits,yunits,title,colfit,pch1,titlesize,textsize)
#piecewise_fit(indata,n,xunits,yunits,title,col1,col2,pch1,titlesize,textsize)
#plot_deannualized_data(indata,n,xunits,yunits,title,colfit,pch1,titlesize,textsize)

# Lomb-Scargle periodogram
#library(cts)
#spec.ls(indata$x, indata$y,spans=NULL,kernel=NULL,taper=0.1,pad=0,fast=T,type="lomb",detrend=T,plot=T)

# plot fit data, etc.
#summary(fit)
#plot(fit,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize)
#fitacf=acf(fit$res, las=1,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize)
#pacf(fit$res, las=1,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize)
#acf2AR(fitacf)