script for processing GAK1 data

gak1 query.R

@@ -0,0 +1,66 @@
+setwd("/Users/MikeLitzow/Documents/R/NSF-GOA")
+require(mgcv)
+require(zoo)
+require(lubridate)
+require(dplyr)
+
+# gak1 data can be downloaded at:
+# http://www.ims.uaf.edu/gak1/data/TimeSeries/gak1.dat
+
+data <- read.csv("gak1 2-17-17.csv")
+head(data)
+
+# appears to be a bad date in the 2013 data - asked Seth Danielson about this (2.17.17)
+# dropping for now
+## NOW THAT CRUISE DATE HAS BEEN CORRECTED!
+# data[data$cruise == "LO128",]
+
+# data <- filter(data, cruise!="LO128")
+
+
+setwd("/Users/MikeLitzow/Documents/R/FATE/fate-ewi")
+
+data$year <- floor(data$dec.year)
+
+data$day <- yday(date_decimal(data$dec.year))
+
+
+head(data)
+
+sal.set <- as.data.frame(tapply(data$sigma.t, list(data$dec.year, data$depth), mean))
+colnames(sal.set) <- c("d0", "d10", "d20", "d30", "d50", "d75", "d100", "d150", "d200", "d250")
+cor(sigma.set, use="pa")
+cor(sal.set, use="pa")
+
+sal.set$year <- floor(as.numeric(rownames(sal.set))) 
+sal.set$day <- yday(date_decimal(as.numeric(rownames(sal.set))))
+
+# now get FMA salinity
+salFMA <- filter(sal.set, day>31 & day <=120)
+
+sal10mu <- tapply(salFMA$d10, salFMA$year, mean, na.rm=T)
+sal20mu <- tapply(salFMA$d20, salFMA$year, mean, na.rm=T)
+sal30mu <- tapply(salFMA$d30, salFMA$year, mean, na.rm=T)
+
+cor(cbind(sal10mu, sal20mu, sal30mu)) # all the same!
+
+ !is.na(salFMA$d30) # no missing values for 30 m salinity!
+sal30day <- tapply(salFMA$day, salFMA$year, mean, na.rm=T)
+sal30n <- tapply(salFMA$d30, salFMA$year, function(x) sum(!is.na(x)))
+
+plot(names(sal30mu), sal30mu, type="o") # getting fresher!
+
+mod <- gam(sal30mu ~ s(sal30day, k=4))
+summary(mod) # best fit is linear!
+
+adj.sal30 <- residuals(mod, type="response") + summary(mod)$p.coeff
+
+# and check
+plot(names(sal30mu), sal30mu, type="o", col="blue")
+lines(names(adj.sal30), adj.sal30, type="o", col= "dark green")
+# only makes a small difference in a couple years
+
+sal.out <- data.frame(raw.sal30=sal30mu, adj.sal30=adj.sal30, mean.day=sal30day)
+
+write.csv(sal.out, "GAK1 30m salinity.csv")
+

goa ketchikan-seward SLP query.R

@@ -1,6 +1,8 @@
 setwd("/Users/MikeLitzow/Documents/R/FATE/fate-ewi")
 library(ncdf4)
 require(chron)
+require(zoo)
+require(dplyr)
 
 # dowload the latest version of the NCEP/NCAR Reanalysis
 # download.file(url="ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/slp.mon.mean.nc", destfile= "/Users/MikeLitzow/Documents/R/NSF-GOA/slp.mon.mean.nc")
@@ -25,5 +27,86 @@ Kx <- ncvar_get(nc, "lon", start=93, count=1)
 Ky <- ncvar_get(nc, "lat", start=15, count=1)
 Kx;Ky
 
+# and check Seward coordinates
+Sx <- ncvar_get(nc, "lon", start=85, count=1)
+Sy <- ncvar_get(nc, "lat", start=13, count=1)
+Sx;Sy
 
-SLP <- ncvar_get(nc, "slp", start=c(77,11,13), count=c(19,8,length(d)), verbose = T)
+Ketchikan <- ncvar_get(nc, "slp", start=c(93,15,1), count=c(1,1,length(d)), verbose = T)
+Seward <- ncvar_get(nc, "slp", start=c(85,13,1), count=c(1,1,length(d)), verbose = T)
+# and we're interested in the difference between the two
+diff <- Ketchikan-Seward
+names(diff) <- d # label by date
+
+# quick plot to check
+plot(d, diff, type="l")
+lines(d, rollmean(diff, k=37, fill=NA), col="red")
+
+      abline(h=mean(diff))
+
+# now summarize
+yrs <- years(d)
+m <- months(d)
+
+# I reckon we should use monthly anomalies to remove the seasonal signal
+# this should produce a more informative annual mean
+mu <- tapply(diff, m, mean)
+sd <- tapply(diff, m, sd)
+
+
+mu <- rep(mu, length.out=length(d))
+sd <- rep(sd, length.out=length(d))
+diff.anom <- (diff-mu)/sd
+
+# first, annual means
+annual <- tapply(diff, yrs, mean)
+plot(names(annual), annual, type="o")
+
+# limit to complete cases
+f <- function(x) sum(!is.na(x)) # function to count number of observations
+use <- tapply(diff, yrs, f)==12
+annual <- annual[use]
+
+# now winter
+win <- c("Nov", "Dec", "Jan", "Feb", "Mar") # define winter months
+win.yrs <- as.numeric(as.character(yrs)) # as 
+
+# align winter years with year corresponding to January
+early.winter <- c("Nov", "Dec")
+add <- m %in% early.winter
+win.yrs[add] <- win.yrs[add]+1
+
+use <- m %in% win
+
+winter.diff <- diff[use]
+win.yrs <- win.yrs[use]
+
+win.mean <- tapply(winter.diff, win.yrs, mean)
+
+# and limit to complete cases
+use <- tapply(winter.diff, win.yrs, f)==5
+win.mean <- win.mean[use]
+# and add a leading NA to fill in 1948
+win.mean <- c(NA, win.mean)
+names(win.mean)[1] <- 1948
+
+# check with a plot
+plot(names(win.mean), win.mean, type="o")
+
+# now spring
+spr <- c("Apr", "May", "Jun") # define spring months
+
+use <- m %in% spr
+
+spr.diff <- diff[use]
+spr.yrs <- yrs[use]
+
+spr.mean <- tapply(spr.diff, spr.yrs, mean)
+
+# and limit to complete cases
+use <- tapply(spr.diff, spr.yrs, f)==3
+spr.mean <- spr.mean[use]
+
+
+# check with a plot
+plot(names(spr.mean), spr.mean, type="o")